Relative Motion Guidance, Navigation and Control for Autonomous Orbital Rendezvous

ABSTRACT:

In this paper, the dynamics of the relative motion problem in a perturbed orbital environment are exploited based on Gauss' variational equations. The relative coordinate frame (Hill frame) is studied to describe the relative motion. A linear high fidelity model is developed to describe the relative motion. This model takes into account primary gravitational and atmospheric drag perturbations. In addition, this model is used in the design of a control, guidance, and navigation system of a chaser vehicle to approach towards and to depart from a target vehicle in proximity operations. Relative navigation uses an extended Kalman filter based on this relative model to estimate the relative position and velocity of the chaser vehicle with respect to the target vehicle and the chaser attitude and gyros biases. This filter uses the range and angle measurements of the target relative to the chaser from a simulated Light Detection and Ranging (LIDAR) system, along with the star tracker and gyro measurements of the chaser. The corresponding measurement models, process noise matrix and other filter parameters are provided. Numerical simulations are performed to assess the precision of this model with respect to the full nonlinear model.The analyses include the navigations errors, trajectory dispersions, and attitude dispersions.

KEYWORDS:
Satellite relative motion; Orbital rendezvous

INTRODUCTION

Although significant progress and technical development have been achieved with regards to orbital rendezvous such as International Space Station supply and repair and automated inspection, servicing, and assembly of space systems, there are limitations with the traditional methods that struggle to meet the new demands for orbital rendezvous. Presently, in order to perform such close proximity operations, mission controllers generally require significant cooperation between vehicles and utilize man-in-the-loop to ensure successful maneuvering of both spacecraft. The interest in autonomous rendezvous and proximity operations has increased with the recent demonstration of XSS-11, Demonstration of Autonomous Rendezvous Technology (DART), and Orbital Express. Autonomous rendezvous and proximity operations have also been demonstrated by Japanese EST-VII, and the Russian Progress vehicles. In addition future missions to the ISS will require autonomous rendezvous and proximity operations (Fehse, 2003; Woffinden and Geller, 2007).

Many relative motion modeling and control strategies have been designed using the linearized Clohessy-Wiltshire (CW) equations to describe the relative motion between satellites. The CW equations are valid if two conditions are satisfied:

The distance between the chaser and the target is small compared tothe distance between the target and the center of the attracting planet; and

The target orbit is near circular (Clohessy and Wiltshire, 1960).

The CW equations do not include any disturbance forces, for example, gravitational perturbations and environmental forces (solar radiation pressure and atmospheric drag). Alternative linear equations that have been used in the literature to model the relative motion are the Tschauner-Hempel (TH) equations (Tschauner and Hempel, 1965).These expressions generalize the CW equations and are similar to them in their derivation and types of applications. Tschauner and Hempel derived theses equations from the viewpoint of rendezvous of a spacecraft with an object in an elliptical orbit. They found complete solutions for elliptical orbits in terms of the eccentric anomaly. This advancement was followed by additional papers which present the complete analytical solution explicit in time, expanding the state transition matrix in terms of eccentricity (Yamanaka and Ankersen, 2002; Carter, 1998; Melton, 2000; Broucke, 2003; Inalhan et al., 2002; Sengupta and Vadali, 2007; Cho and Park, 2009). This form of solution is used to analyze the relative motion between the chaser and the target vehicles in the relative frame of motion more efficiently and rapidly than solving the exact nonlinear differential equations in the inertial coordinate system. The TH equations do not take into account any perturbation forces. These perturbations have a significant effect on the satellite relative motion.

Due to the previous limitations of the CW and TH models, this paper proposes an innovative linear model which includes both the perturbation that reflects the Earth's oblateness effect and atmospheric drag perturbation in the Cartesian coordinates orbital frame with little complication. Especially in low Earth orbits (LEOs), these perturbations have a deep influence on the relative dynamics, and their inclusion in the linear model can sensibly increase the performance of the linear filters, allowing greater insight of satellite relative motion, and providing an opportunity to investigate alternative feedback control strategies for the proximity operations.

Unlike the relative translation motion control, the relative rotational control is a traditional feedback control system. During the mission scenarios, the chaser vehicle may need to track the target vehicle to achieve proper docking maneuvers and or visual inspection tasks. The paper uses an extended Kalman filter formulation to estimate the relative motion and chaser attitude using range and angle measurements from a LIDAR system coupled with gyro and star tracker measurements of the chaser (Woffinder and Geller, 2007; Jenkins and Geller, 2007; Junkins et al., 2005; Woffinden, 2004). The Kalman filter basically consists of two main stages. The first stage is the propagation stage, where the states are propagated numerically and it is based on the proposed linear model. The second stage comes when the measurements from the sensors are available and it is used to update the states of the first stage. The corresponding measurement models, process noise matrix, and other filter parameters are provided. Momentum wheels are assumed for attitude control and thrusters are assumed for translation control. The effects of the navigation filter, pointing algorithms, and control algorithms are included in the analysis.

The objective of this paper is as follows:

• To develop linearized high fidelity models for relative motion in a perturbed orbit;

• To design a navigation filter that can determine the relative position and velocity between target and chaser vehicles as well as orientations and angular rates of the chaserthat support closed-loop proximity attitude control operations and maneuvers; and

• To design a control system for the chaser vehicle to either approach or depart fromthe target vehicle in proximity operations in a general perturbed orbit for coupled translation and rotation relative motion.

The analysis in the current paper is summarized as follows. First, we present the relative dynamics equation of motion for the chaser with respect to the target in a general perturbed orbit, along with attitude dynamics models.Next, a linear high fidelity relative motion model is developed to describe relative motion in proximity operations based on Gauss' variational method. Then, the relative navigation and an extended Kalman filter are presented for the relative motion and attitude estimations, along with the relative translational and rotational controller. In the simulation section, the accuracy and performanceof the relative navigation and controller, based on the high fidelity model, are illustrated through different numerical examples and comparisons are made with the truth nonlinear model. Finally, conclusion of the work is presented and suggestions are made for future work.

RELATIVE MOTION MODELS

Consider an Earth-centered inertia (ECI) frame, with orthonormal basis {i_X, i_Y, i_Z,}. The vectors i_X and i_Y lie in the equatorial plane, with i_X coinciding with the line of equinoxes, and i_Z passing through the North Pole. Relative motion is conveniently described in a Local-Vertical-Local-Horizontal (LVLH) frame,which is attached to the target spacecraft, as shown in Fig.1. This frame has basis {i_X, i_Y, i_Z,} with i_X lying along the radius vector form the Earth's center to the spacecraft, i_Z coinciding with the normal to the plane defined by the position and velocity vectors of the spacecraft, and i_Y = i_Z x i_X. The LVLH frame rotates with angular velocityvector ω, and its current orientation with respect to the ECI frame is given by the 3-1-3 direction cosine matrix, comprising right ascension of ascending node Ω, inclination i, perigee argument ω plus true anomaly f, respectively (Fig.2). The angular velocity can also be expressed in terms of orbital elements and their rates.

Figure 1
Relative Motion Coordinates.

Figure 2
Orbital Elements.

Let the position of the chaser vehicle in the target's LVLH frame be denoted by ρ=xi_x+yi_y+zi_z, where x, y and z denote the components of the position vector along the radial, transverse, and out-of-plane directions, respectively. ρis determined from ρ=R_c -R_t, where R_c and R_t are the chaser and target absolute position vectors. Then, the most general equations modeling relative motion are given by the following:

(1) $\ddot{\mathit{\rho }}=[{\mathit{f}}_{\mathit{c}}{]}^{\mathit{LVLH}}-[{\mathit{f}}_{\mathit{c}}{]}^{\mathit{LVLH}}-2\mathit{\omega }\times \dot{\mathit{\rho }}-\mathit{\omega }\times \mathit{\omega }\times \mathit{\rho }-\dot{\mathit{\omega }}\times \mathit{\rho }$

where [f_c]^LVLH and [f_t]^LVLH are the external accelerations acting on the chaser and the target, respectively in the LVLH frame of the target vehicle. In Eq. (1), (..) and (.) denote the first and second derivatives with respect to time.

It is assumed, in this paper, that the externalaccelerations arise due to two basic groups of accelerations, defined by the following equation:

(2) $f=f_g+f_a+f_c+f_w$

The first group of accelerationsis due to gravitational effects, f_g, atmospheric drag, f_a, and control, f_c. Since Earth isn't perfectly spherical, more accurate gravity models exist, taking into account Earth's irregular shape. One irregularity that has a significant influence on space missions is the Earth's bulge at the equator. This phenomenon is captured in the J₂ gravity model (Vallado, 2001; Schaub and Junkins, 2003). The second group of accelerations, f_w, is considered to be small accelerations, due to the gravity fields of other planets, solar pressure, or venting, which also perturbs the spacecraft's motion. These small accelerationsare grouped together and modeled as zero mean normally distributed random variables (Woffinden and Geller, 2007).

In the literature, the most popular methods to model the spacecraft's orbit are known as Cowell's method and Gauss's method (Vallado, 2001; Schaub and Junkins, 2003). The Cowell's method is basically defined by specifying the position (R) and velocity (V) vectors of the spacecraft in the inertial coordinate frame, while Gauss' method is defined by an equivalent set of elements called orbital elements (a,e,i,Ω,ω,f) which correspond to the semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of periapsis, and true anomaly, respectively, as shown in Fig. 2.

Table 1 summarizes the dynamic equations that are used in order to describe all of these methods. In this table, [.]^I and [.]^LVLH denote that the forces are defined in inertial and LVLH coordinate frames, respectively; a_x, a_y and a_z are the components of disturbance accelerations acting on the target in the LVLH reference coordinate frame; s_(.)=sin(.) and c_(.)=cos(.); and R⊕ are the Earth gravitational constant and the radius of the Earth; the terms R and V refer to the magnitude of the position and velocity vectors, respectively; the quantity H denotes the magnitude of the specific angular momentum vector defined by H=R×V; X, Y and Z are the components of the spacecraft position vector; C_D is the atmospheric drag coefficient; A denotes the cross sectional area; m is the spacecraft mass; and finally, ρ is the atmospheric density. Exponential atmospheric behavior is used to model Earth's atmospheric density. This model and its corresponding parameters are defined inVallado (2001).

Table 1
Orbit Model Methods Summary.

Method	Dynamic equations
Cowell's	$\begin{array}{c}\ddot{R}=f\hfill \\ f=f_g+f_a+f_c+f_w\hfill \\ f_g=a_{t\mathrm{wo}-\mathrm{body}}+a_{J_2}\hfill \\ a_{\mathrm{two}-\mathrm{body}}=-(\frac{\mu }{R^3})R\hfill \\ [a_{J_2}{]}^I=-\frac{3}{2}{J}_2\frac{\mu R_{\oplus }^2}{R^4}\mathrm{[(}1-5(Z/R{)}^2)\frac{X}{R}\left(1-5{\left(Z/R\right)}^2\right)\frac{Y}{R}\left(3-5{\left(Z/R\right)}^2\right)\frac{Z}{R}{]}^T\hfill \\ [f_a{]}^I=-\frac{1}{2}\frac{C_{D}A}{m}\rho VV\hfill \end{array}$
Gauss's	$\begin{array}{c}\frac{\mathrm{da}}{dr}=\frac{2a^2}{H}[e{s}_f{a}_{\chi }+\frac{p}{R}{a}_y]\hfill \\ \frac{\mathrm{de}}{\mathrm{dr}}=\frac{1}{H}\{p{s}_f{a}_{\chi }+\mathrm{[(}p+R)c_f+\mathrm{Re}]a_y\}\hfill \\ \frac{\mathrm{de}}{d\iota }=\frac{1}{H}\{p{s}_f{a}_{\chi }+\mathrm{[(}p+R)c_f+\mathrm{Re}]a_y\}\hfill \\ \frac{\mathrm{di}}{\mathrm{dr}}=\frac{R{c}_{\omega +f}}{H}{a}_Z\hfill \\ \frac{d\omega }{\mathrm{dl}}=\frac{1}{\mathrm{He}}\{-p{c}_f{a}_{\chi }+(p+R)s_f{a}_y\}-\frac{R{s}_{\omega +f^{C}i}}{H{s}_i}{a}_Z\hfill \\ \frac{\mathrm{df}}{\mathrm{dl}}=\frac{H}{R^2}+\frac{1}{\mathrm{He}}[p{c}_f{a}_{\chi }-(p+R)s_f{a}_y]\hfill \\ [a{]}^{\mathrm{LVLH}}=[a_{J_2}+f_a+f_c+f_w{]}^{\mathrm{LVLH}}\hfill \\ [a_{j_2}{]}^{\mathrm{LVLH}}=-\frac{3}{2}{J}_2\frac{\mu R_{\oplus }^2}{R^4}[1-3s_i^2{s}_{\omega +f}^2{s}_i^2{s}_{2(\omega +f)}{s}_{2i}{s}_{\omega +f}{]}^T\hfill \\ [f_a{]}^{\mathrm{LVLH}}=-\frac{1}{2}\frac{C_{D}A}{m}\mathit{\rho }V{V}^{\mathrm{LVLH}}\hfill \\ V^{\mathrm{LVLH}}=\frac{H}{p}[{\mathrm{es}}_{f}p/R0{]}^T\hfill \end{array}$

In order to use the generalized relative dynamic model defined by Eq. (1), the angular velocity vector, ω, and the angular acceleration vector, ώ, ofthe LVLH frame with respect to the ECI frame, needs to be determined. Table 2 summarizes the equations that can be used to compute these vectors. These equations are derived based on using either Cowell's method (position and velocity vectors) or Gauss's method (orbital elements). In this table, the matrix T_I^LVLH denotes the direction cosine matrix of the LVLH coordinate frame with respect to the ECI coordinate frame.

Table 2
LVLH Coordinate Frame Orientation.

	Dynamic equations
Given Inertial Position and Velocity	$\begin{array}{c}{T}_I^{\mathrm{LVLH}}=[i_x{i}_y{i}_z{]}^T\hfill \\ i_x=\frac{R}{R}/i_Z=\frac{H}{H},i_y=i_Z\times i_x\hfill \\ \mathbf{\omega }=\left(\begin{array}{c}R[f_Z{]}^{\mathrm{LVLH}}/H\\ 0\\ H/R^2\end{array}\right),\dot{\mathbf{\omega }}=\left(\begin{array}{c}\frac{(\dot{R}\left[f_Z{]}^{\mathrm{LVLH}}+R\right.\left[{{\displaystyle \dot{f}}}_Z{]}^{\mathrm{LVLH}})H-R\right.\left[f_Z{]}^{\mathrm{LVLH}}\dot{H}\right.}{H^2}\\ 0\\ \frac{\dot{\mathrm{HR}}-2H\dot{R}}{R^3}\end{array}\right)\hfill \\ [\dot{f}{]}^{\mathrm{LVLH}}=T_I^{\mathrm{LVLH}}[\dot{f}{]}^I-\mathbf{\omega }\times T_I^{\mathrm{LVLH}}[f{]}^I\hfill \\ \dot{H}=R{f}_{y\prime }\dot{R}=\frac{R·V}{R},H\mathbf{=}R\mathbf{\times }V\hfill \end{array}$
Given Orbital Elements	$\begin{array}{c}{T}_I^{\mathrm{LVLH}}=T^3(\omega +f)T^1(i)T^3(\Omega )=\left\{\begin{array}{ccc}{c}_{\Omega }{c}_{\omega +f}-s_{\Omega }{s}_{\omega +f^{C}i}& S_{\Omega }{C}_{\omega +f}+c_{\Omega }{s}_{\omega +f^{C}i}& s_{\omega +f^{S}i}\\ -c_{\Omega }{s}_{\omega +f}-s_{\Omega }{c}_{\omega +f^{C}i}& -S_{\Omega }{S}_{\omega +f}+c_{\Omega }{c}_{\omega +f^{C}i}& c_{\omega +f^{S}i}\\ S_{\Omega }{S}_i& -C_{\Omega }{S}_i& c_i\end{array}\right\}\hfill \\ \mathbf{\omega }=T^3(\omega +f)T^1(i)\left\{\begin{array}{c}0\\ 0\\ \dot{\Omega }\end{array}\right\}+T^3(\omega +f)+\left\{\begin{array}{c}\raisebox{1ex}{$\mathrm{di}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.\\ 0\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ 0\\ \dot{\mathit{\omega }}+f\end{array}\right\}=\left(\begin{array}{c}{s}_{\mathit{\omega }+f}{s}_j\dot{\Omega }+c_{\mathit{\omega }+f}\raisebox{1ex}{$\mathrm{di}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.\\ 0\\ c_i\dot{\Omega }+(\dot{\omega }+f)\end{array}\right)\hfill \\ T^1(\theta )=\left\{\begin{array}{ccc}1& 0& 0\\ 0& c_{\theta }& s_{\theta }\\ 0& -S_{\theta }& C_{\theta }\end{array}\right\},T^2(\theta )=\left\{\begin{array}{ccc}{C}_{\theta }& 0& -S_{\theta }\\ 0& 1& 0\\ \mathrm{Sg}& 0& C_{\theta }\end{array}\right\},T^3(\theta )=\left\{\begin{array}{ccc}{c}_{\theta }& s_{\theta }& 0\\ -s_{\theta }& c_{\theta }& 0\\ 0& 0& 1\end{array}\right\}\hfill \end{array}$

The Euler's equation of motion is used to describe the attitude dynamics for both target and chaser vehicles, and a quaternion formulation is used for attitude kinematics. The dynamics for both vehicles are given below as (Woffinden and Geller, 2007):

(3) ${\dot{\mathit{q}}}_{\mathit{l}}^{\mathit{t}}=\frac{1}{2}{\omega }^{\mathit{t}}\otimes {\mathit{q}}_{\mathit{l}}^t$

(3b) ${{\displaystyle \dot{\omega }}}^t={\mathit{I}}_{\mathit{t}}^{-1}[{\tau }_{\mathit{T}}-{\omega }^{\tau }\times {\mathit{I}}_{\mathit{T}}{\omega }^t]$

(3c) ${\mathit{\tau }}_t={\tau }_{{\tau }_{\mathit{g}}}+{\mathit{T}}_{{\mathit{t}}_{\mathit{d}}}$

(4a) ${\dot{\mathit{q}}}_{\mathit{l}}^{\mathit{c}}=\frac{1}{2}{\omega }^{\mathit{c}}\otimes {\mathit{q}}_{\mathit{l}}^{\mathit{c}}$

(4b) ${{\displaystyle \dot{\omega }}}^{\mathit{c}}={\mathit{I}}_{\mathit{c}}^{-1}[{\tau }_{\mathit{c}}-{\omega }^{\mathit{c}}\times {\mathit{I}}_{\mathit{c}}{\omega }^{\mathit{c}}]$

(4c) ${\tau }_{\mathit{C}}={\tau }_{{\mathit{c}}_{\mathit{C}}}+{\mathit{\tau }}_{{\mathit{c}}_{\mathit{g}}}+{\mathit{\tau }}_{{\mathit{c}}_{\mathit{d}}}$

where ⊕ is the quaternion multiplication operator defi ned by Lear (1985).

(5) $\mathit{p}\otimes \mathit{q}=\left[\begin{array}{cccc}{\mathit{p}}_4& {\mathit{p}}_3& -{\mathit{p}}_2& {\mathit{p}}_1\\ -{\mathit{p}}_3& {\mathit{p}}_4& {\mathit{p}}_1& {\mathit{p}}_2\\ {\mathit{p}}_2& -{\mathit{p}}_1& {\mathit{p}}_4& {\mathit{p}}_3\\ -{\mathit{p}}_1& -{\mathit{p}}_2& -{\mathit{p}}_3& {\mathit{p}}_4\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_1\hfill \\ {\mathit{q}}_2\hfill \\ {\mathit{q}}_3\hfill \\ {\mathit{q}}_4\hfill \end{array}\right]$

and the i vehicle gravity gradient torque is defined by

(6) ${\mathit{T}}_{\mathit{ig}}=3\frac{\mathit{\mu }}{{\left|{\mathit{R}}_{\mathit{i}}\right|}^5}({\mathit{R}}_{\mathit{i}}\times {\mathit{I}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}})$

In Eqs. 3 and 4, the target states include the quaternion, q_I^t, that defines the orientation of the target with respect to the inertial frame, and the target's angular rate, ^_ωt. Similarly the chaser states are q_I^c and I_c. andare the target and chaser inertia matrices, respectively. The gravity gradient torque, τ_(ig), for both vehicles for the chaser and τ_(tg) for the target) is derived from the ^(_τ(cg) Earth'spoint mass gravity models. The random disturbances, τ_(td) and, τ_(tg) are included in the models to account for disturbance torques acting on each vehicle such as drag, solar radiation and other unmodeled disturbances. These unmodeled disturbances are represented as uncorrelated white noise, with mean and variance defined by a trial and error technique outlined by Lear (1985). The control input, τ_(cc), is the torque executed by the actuators (momentum wheels) on the chaser spacecraft.

It is assumed that the available sensors are the LIDAR for tracking the target and an assembly of a star tracker and gyros for attitude determination. The parameter states for these sensors are modeled as first-order Markov processes with large time constants, causing them to behave like biases. The parameter states include the gyros bias b_ω^c, star camera misalignments Є_s^s, and LIDAR misalignments Є_l^l. The dynamics model associated with these states is given by:

(7a) ${\dot{\mathit{b}}}_{\mathit{\omega }}^{\mathit{c}}=-\frac{{\mathit{b}}_{\omega }^{\mathit{c}}}{{\tau }_{\mathit{b}}^{\omega }}+{\mathit{w}}_{\mathit{b}}^{\omega }$

(7b) ${\dot{\varepsilon }}_{\mathit{s}}^{\mathit{s}}=-\frac{{\varepsilon }_{\mathit{s}}^{\mathit{s}}}{{\mathit{\tau }}^{\mathit{S}}}+{\mathit{w}}^{\mathit{s}}$

(7c) ${\dot{\varepsilon }}_{\mathit{l}}^{\mathit{l}}=-\frac{{\varepsilon }_{\mathit{l}}^{\mathit{l}}}{{\tau }^{\mathit{l}}}+{\mathit{w}}^{\mathit{l}}$

where, w_b^ω, and w^s and ^_wl are white noise terms, driving the firstorder Markov processes and τ_b^ω, τ^s and τ^l are the corresponding time constants.

The actuator models used in the simulation include momentum wheels for orientation control and thrusters for translational control. The mathematical model for the actual control torque, generated by the wheels, and the impulsive thrust, by the thrusters, are:

(8) ${\mathit{\tau }}_{\mathit{CC}}=\mathit{\delta }\mathit{T}({\varepsilon }_{\tau }^{\mathit{c}})\left[\left\{{\mathit{I}}_{3\times 3}+\mathit{Diag}\left({\mathit{f}}_{\mathit{\tau }}^{\mathit{c}}\right)\right\}{{\displaystyle \hat{\mathit{\tau }}}}_{\mathit{cc}}+{\mathit{b}}_{\mathit{\tau }}^{\mathit{c}}+{\mathit{v}}_{\mathit{\tau }}^{\mathit{c}}\right]$

(9) ${\mathit{u}}_{{\mathit{C}}_{\Delta \mathit{V}}}=\delta \mathit{T}({\varepsilon }_{\Delta \mathit{V}}^{\mathit{c}}\mathrm{)[\{}{\mathit{I}}_{3\times 3}+\mathit{Diag}({\mathit{f}}_{\Delta \mathit{V}}^{\mathit{c}}\mathrm{)\}}\Delta {\mathit{V}}_{\mathit{c}}+{\mathit{b}}_{\Delta \mathit{V}}^{\mathit{c}}+{\mathit{v}}_{\Delta \mathit{V}}^{\mathit{c}}]$

The generated torque and impulsive include errors such as noises ν^c, biases b^c, scale factor biases f^c, and misalignments Є^c. These errors can be modeled also as white noises.

The simulation contains gyros, star tracker, and LIDAR sensor models. The models for these measurements are given by:

Gyro Model:

(10) ${{\displaystyle \tilde{\mathit{\omega }}}}^{\mathit{c}}=\mathit{\delta }\mathit{T}({\mathit{\varepsilon }}_{\omega }^{\mathit{c}}\mathrm{)[\{}{\mathit{I}}_{3\times 3}+\mathit{Diag}({\mathit{f}}_{\mathit{\omega }}^{\mathit{c}}\mathrm{)\}}{\mathit{\omega }}^{\mathit{c}}+{\mathit{b}}_{\omega }^{\mathit{c}}+{\mathit{v}}_{\omega }^{\mathit{c}}]$

Star Tracker Model:

(11) ${\tilde{\mathit{q}}}_{\mathit{I}}^{\mathit{S}}=\delta \mathit{q}({\mathit{v}}_{\mathit{s}}^{\mathit{s}})\otimes \mathit{q}({\varepsilon }_{\mathit{s}}^{\mathit{s}})\otimes {\mathit{q}}_{\mathit{c}}^{\mathit{s}}\otimes {\mathit{q}}_{\mathit{J}}^{\mathit{c}}$

LIDAR Model:

(12) $\left(\begin{array}{c}\tilde{\mathit{\alpha }}\\ \tilde{\mathit{\beta }}\\ \tilde{\mathit{\rho }}\end{array}\right)=\left(\begin{array}{c}{\tan }^{-1}(\frac{{\mathit{i}}_{\mathit{Y}}}{{\mathit{i}}_{\chi }})\\ {\sin }^{-1}({\mathit{i}}_{\mathit{Z}})\\ \mathit{\rho }\end{array}\right)+\left(\begin{array}{c}{\mathit{v}}_{\mathit{\alpha }}\\ {\mathit{v}}_{\mathit{\beta }}\\ {\mathit{v}}_{\mathit{\rho }}\end{array}\right)$

where:

(13) ${\mathit{i}}_{\mathit{los}}^{\mathit{l}}=\left(\begin{array}{c}{\mathit{i}}_{\chi }\\ {\mathit{i}}_{\mathit{y}}\\ {\mathit{i}}_{\mathit{Z}}\end{array}\right)=\left(\begin{array}{c}{\mathit{c}}_{\mathit{\alpha }}{\mathit{c}}_{\mathit{\beta }}\hfill \\ {\mathit{s}}_{\mathit{\alpha }}{\mathit{c}}_{\mathit{\beta }}\hfill \\ {\mathit{s}}_{\mathit{\beta }}\hfill \end{array}\right)={\mathit{T}}^{\mathit{l}\hat{\mathit{l}}}{\mathit{T}}^{\hat{\mathit{l}}\mathit{s}}{\mathit{T}}^{\mathit{s}\hat{\mathit{s}}}{\mathit{T}}^{\hat{\mathit{s}}\mathit{l}}{\mathit{T}}^{\mathit{lT}}(-\frac{\rho }{\rho })$

The gyro models include bias b_ω^c, scale factor bias f_ω^c, and angular random walk noise ν_ω^c. The starcamera model accounts for the uncertainty in the alignment of the star camera frame ∈_s^s with respect to the chaser frame and sensor noise ν_s^s. The q_c^s refers to the fixed orientation of the star camera coordinate frame with respect to the chaser body coordinate frame. The LIDAR model includes angle measurements (azimuth, α, and elevation, β) noises v_α, v_β and range (ρ) noise, v_ρ. The transformation matrix denoted by ^Tab is the transformation matrix used to transform any vector from coordinate b to coordinate a. The term i_los^l is the line of sight vector in the LIDAR coordinate frame (Fig. 3). The transformations T^lI, T^IS, T^ŜS, T^ŜI, and T^IT are a series of transformation matrices to transform the line of sight vector from target LVLH coordinate frame to the LIDAR coordinate frame. These transformations include errors from sensor misalignments, noises, and attitude determination errors.

Figure 3
Line of Sight Vector.

The small angle rotations can be written in terms of quaternions as

(14) $\mathbf{\delta }\mathit{q}(\mathbf{\theta })\approx \left[\frac{\mathit{\theta }/2}{1}\right]$

or attitude matrices as

(15) $\mathbf{\delta }\mathit{T}(\mathbf{\theta })\approx \mathit{I}-{\mathbf{\theta }}^{\times }$

where θ=θu is a small rotation vector, and θ^× operating on vector ω is a cross product matrix defined by the ordinary cross product θ^× ω = ω × θ.

LINEAR HIGH FIDELITY RELATIVE MODEL

In this section, a linear time varying high fidelity model is obtained to describe the relative motion dynamics. This model is derived based on two main assumptions. The first assumption is that the relative distant between the chaser and the target vehicles is much less than the target orbital radius. The second one assumes that the main disturbance accelerations, that affect both vehicles are the gravitational acceleration and the atmospheric drag acceleration. Based on these assumptions, all terms mentioned in the general relative dynamic expression, Eq. (1), are expanded considering only first order terms to obtain the new proposed model. Table 3 summarizes the procedures that have been followed to obtain this model. In this table the linear time varying model reduces to the following form

(16) $\dot{\mathit{x}}=\mathit{Ax}$

Table 3
Relative Orbit Model Summary.

Model	Equations
Nonlinear	$\ddot{\mathit{\rho }}=[f_c{]}^{\mathrm{LVLH}}-[f_t{]}^{\mathrm{LVLH}}-2\mathbf{\omega }\times \ddot{\mathit{\rho }}-\mathbf{\omega }\times \mathbf{\omega }\times \mathit{\rho }-\dot{\mathbf{\omega }}\times \mathit{\rho }$
Linear Time Varying	$\begin{array}{c}{T}_I^{\mathrm{LVLH}}=T^3(\omega +f)T^1(i)T^3(\Omega )=\left\{\begin{array}{ccc}{c}_{\Omega }{c}_{\omega +f}-s_{\Omega }{s}_{\omega +f^{C}i}& s_{\Omega }{c}_{\omega +f}+c_{\Omega }{s}_{\omega +f^{C}i}& s_{\omega +f^{S}i}\\ -c_{12}{s}_{\omega +f}-s_{\Omega }{c}_{\omega +f}{c}_i& -s_{\Omega }{s}_{\omega +f}+c_{\Omega }{c}_{\omega +f}{c}_i& c_{\omega +f}{s}_i\\ s_{\Omega }{s}_i& -c_{\Omega }{s}_i& c_i\end{array}\right\}\hfill \\ \mathbf{\omega }=T^3(\omega +f)T^1(i)\left\{\begin{array}{c}0\\ 0\\ \dot{\Omega }\end{array}\right\}+T^3(\omega +f)\left[\begin{array}{c}\raisebox{1ex}{$\mathrm{di}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.\\ 0\\ 0\end{array}\right]+\left\{\begin{array}{c}0\\ 0\\ \dot{\omega }+\dot{f}\end{array}\right\}=\left(\begin{array}{c}{s}_{\mathit{\omega }+f}{s}_i\dot{\Omega }+c_{\mathit{\omega }+f}\raisebox{1ex}{$\mathrm{di}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.\\ 0\\ c_i\dot{\Omega }+(\dot{\omega }+\dot{f})\end{array}\right)\hfill \\ T^1(\theta )=\left\{\begin{array}{ccc}1& 0& 0\\ 0& c_{\theta }& s_{\theta }\\ 0& -s_{\theta }& c_{\theta }\end{array}\right\},T^2(\theta )=\left\{\begin{array}{ccc}{c}_{\theta }& 0& -s_{\theta }\\ 0& 1& 0\\ s_{\theta }& 0& c_{\theta }\end{array}\right\},T^3(\theta )=\left\{\begin{array}{ccc}{c}_{\theta }& s_{\theta }& O\\ -s_{\theta }& c_{\theta }& O\\ 0& 0& 1\end{array}\right\}\hfill \\ \ddot{\mathbf{\rho }}=[f_c{]}^{\mathrm{LVLH}}-[f_t]-2\mathbf{\omega }\times \dot{\mathbf{\rho }}-\mathbf{\omega }\times \mathbf{\omega }\times \mathbf{\rho }-\dot{\mathbf{\omega }}\times \mathit{\rho }-\hfill \\ 2\mathbf{\omega }\times \dot{\mathbf{\rho }}-\mathbf{\omega }\times \mathbf{\omega }\times \mathbf{\rho }-\dot{\mathbf{\omega }}\times \mathbf{\rho }=\left\{\begin{array}{cccccc}{\omega }_Z^2& {{\displaystyle \dot{\omega }}}_Z& -{\omega }_x{\omega }_Z& 0& 2{\omega }_Z& 0\\ -{{\displaystyle \dot{\omega }}}_Z& {\omega }_x^2+{\omega }_Z^2& {{\displaystyle \dot{\omega }}}_{\chi }& -2{\omega }_Z& 0& 2{\omega }_x\\ -{\omega }_{\chi }{\omega }_Z& -{{\displaystyle \dot{\omega }}}_{\chi }& {\omega }_{\chi }^2& 0& -2{\omega }_{\chi }& 0\end{array}\right\}\left(\begin{array}{c}x\\ y\\ Z\\ \dot{x}\\ \dot{y}\\ \dot{Z}\end{array}\right)\hfill \\ [f_c]-[f_t]\approx [\Delta f_g+\Delta f_a{]}^{\mathrm{LVLH}}\hfill \\ [\Delta a_{\mathrm{lwo}-\mathrm{body}}{]}^{\mathrm{LVLH}}\approx \frac{\mu }{R^3}\left(\begin{array}{ccc}2& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right)\left(\begin{array}{c}\chi \\ y\\ z\end{array}\right)\hfill \\ [\Delta a_{J_2}{]}^{\mathrm{LVLH}}\approx -\frac{3}{2}{J}_2\frac{\mu R_{\oplus }^2}{R^5}\left\{\begin{array}{ccc}-4(1-3i_{x_3}^2)& -8i_{x_3}{i}_{y_3}& -8i_{x_3}{i}_{z_3}\\ -8i_{x_3}{i}_{y_3}& +2(1-i_{z_3}^2)1-7i_{x_3}^2& 2i_{y_3}{i}_{z_3}\\ -8i_{x_3}{i}_{z_3}& 2i_{y_3}{i}_{z_3}& -2(1-i_{z_3}^2)3-5i_{x_3}^2\end{array}\right\}\left(\begin{array}{c}\chi \\ \gamma \\ z\end{array}\right)\hfill \\ [\Delta f_a{]}^{\mathrm{LVLH}}\approx \left\{\begin{array}{ccc}{C}_c{\omega }_{Z}R& 0& 0\\ 0& C_c{\omega }_{Z}R& 0\\ 0& 0& C_c{\omega }_{Z}R\end{array}\right\}\left(\begin{array}{c}\dot{\chi }\\ \dot{y}\\ \dot{Z}\end{array}\right)\hfill \\ \dot{x}=\mathrm{Ax},x=[xyz\dot{x}\dot{y}\dot{z}{]}^T\hfill \\ A\mathbf{=}\left\{\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right\}{A}_{11}\mathbf{=}{0}_{3\mathbf{\times }3},A_{12}\mathbf{=}{I}_{3\mathbf{\times }3}\hfill \\ A_{21}=\left\{\begin{array}{ccc}{\mathit{\omega }}_z^2+2c_1-4C_2\left(1-3i_{x_3}^2\right)& {\dot{\mathit{\omega }}}_z-8C_2{i}_{x_3}{i}_{y_3}-C_c{\mathit{\omega }}_z^{2}R& -{\mathit{\omega }}_x{\mathit{\omega }}_z-8C_2{i}_{x_3}{i}_{z_3}\\ -{\dot{\mathit{\omega }}}_z-8C_2{i}_{x_3}{i}_{x_{y3}}+C_c{\mathit{\omega }}_z^{2}R& \left({\mathit{\omega }}_x^2+{\mathit{\omega }}_z^2\right)-C_1+C_2\left(1-7i_{x_3}^2+2\left(1-i_{z_3}^2\right)\right)& {\dot{\mathit{\omega }}}_x+2C_2{i}_{y_3}{i}_{z_3}-C_c{\mathit{\omega }}_x{\mathit{\omega }}_{z}R\\ -{\mathit{\omega }}_x{\mathit{\omega }}_z-8C_2{i}_{x_3}{i}_{z_3}& -{\dot{\mathit{\omega }}}_x+2C_2{i}_{y3}{i}_{z_3}+C_c{\mathit{\omega }}_x{\mathit{\omega }}_{z}R& {\mathit{\omega }}_x^2-C_1+C_2\left(3-5i_{x_3}^2-2\left(1-i_{z_3}^2\right)\right)\end{array}\right\}\hfill \\ A_{22}=\left\{\begin{array}{ccc}{C}_c{\omega }_{Z}R& 2{\omega }_Z& 0\\ -2{\omega }_Z& C_c{\omega }_{Z}R& 2{\omega }_{\chi }\\ 0& -2{\omega }_{\chi }& C_c{\omega }_{Z}R\end{array}\right\}\hfill \\ C_1=\frac{\mu }{R^3},C_2=-\frac{3}{2}{J}_2\frac{\mu R_{\oplus }^2}{R^5},C_c=-\frac{1}{2}\frac{C_{D}A}{m}\rho \hfill \end{array}$

where x is the state vector. This model can be used to approximate the time varying state transition matrix by expanding the time invariant exponential matrix solution in a Taylor series to fourth order, as follows:

(17) ${\mathit{\phi }}_{\mathit{LTV}}={\mathit{e}}^{\mathit{A}\Delta \mathit{t}}\approx \mathit{I}+\mathit{A}\Delta \mathit{t}+\frac{{\mathit{A}}^2\Delta {\mathit{T}}^2}{2!}+\frac{{\mathit{A}}^3\Delta \mathit{T}3}{3!}+\frac{{\mathit{A}}^4\Delta {\mathit{T}}^4}{4!}$

This matrix is used in the next section as a part of the extended Kaman filter, to propagate the states forward in time and to compute the filter parameters.

NAVIGATION CONTROL MODEL ALGORITHMS

The main objective of the navigation system is to estimate the target's relative position, relative velocity and orientation given noisy sensor measurements, imperfect dynamic models, and uncertain initial conditions. The logic behind the navigation filter is to process information collected from sensors and various mathematical models to generate the best possible estimation of the states. Space navigation application of the Kalman filter is presented in this section. The dynamic models for a closed loop GN&C system are shown in Fig. 4.

Figure 4
Closed Loop GN&C System.

The navigation model uses an extended Kalman filter to estimate the relative position and velocity of the chaser vehicle with respect to the target vehicle, and the approximated analytical state transition matrix solution. Orbital elements of the target are numerically propagated with respect to time using Gauss's variational equations, with J₂ and drag perturbations. These orbital elements are used to compute the transformation matrix of the target vehicle with respect to the inertial frame, as well as to assist in estimating LIDAR measurements. The dynamic models used to propagate the navigation states are:

(18a) $\frac{\mathit{d}\hat{\mathit{a}}}{\mathit{dT}}=\frac{2{\hat{\mathit{a}}}^2}{\hat{\mathit{H}}}\left[{\widehat{\mathit{es}}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\mathit{x}}+\frac{\hat{\mathit{p}}}{\hat{\mathit{R}}}{\hat{\mathit{a}}}_{\mathit{y}}\right]$

(18b) $\frac{\mathit{d}\hat{\mathit{e}}}{\mathit{dT}}=\frac{1}{\hat{\mathit{H}}}\{{\widehat{\mathit{ps}}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\mathit{x}}+\mathrm{[(}\hat{\mathit{p}}+\hat{\mathit{R}}){\mathit{c}}_{\hat{\mathit{f}}}+\hat{\mathit{R}}\hat{\mathit{e}}]{\hat{\mathit{a}}}_{\mathit{y}}\}$

(18c) $\frac{\mathit{d}\widehat{i}}{\mathit{dt}}=\frac{\widehat{R}{c}_{\widehat{\omega }+\widehat{f}}}{\widehat{H}}{\widehat{\mathit{a}}}_{\mathit{z}}$

(18d) $\frac{\mathit{d}\widehat{\Omega }}{\mathit{dt}}=\frac{\widehat{R}{s}_{\widehat{\omega }+\widehat{f}}}{\widehat{H}\mathit{s}\widehat{i}}{\widehat{\mathit{a}}}_{\mathit{z}}$

(18e) $\frac{\mathit{d}\hat{\mathit{\omega }}}{\mathit{dt}}=\frac{1}{\hat{\mathit{H}}\hat{\mathit{e}}}\{-{\widehat{\mathit{pc}}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\varkappa }+(\hat{\mathit{p}}+\hat{\mathit{R}}){\mathit{s}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\mathit{y}}\}-\frac{{\mathit{R}}_{\mathit{s}\hat{\mathit{\omega }}}{+}_{\hat{\mathit{f}}}{\mathit{c}}_{\hat{\mathit{i}}}}{{\widehat{\mathit{Hs}}}_{\hat{\mathit{i}}}}{\hat{\mathit{a}}}_{\mathit{z}}$

(18f) $\frac{\mathit{d}\hat{\mathit{f}}}{\mathit{dt}}=\frac{\hat{\mathit{H}}}{\mathit{R}}+\frac{1}{\hat{\mathit{H}}\hat{\mathit{e}}}[{\widehat{\mathit{pc}}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\mathit{\chi }}-(\hat{\mathit{p}}+\hat{\mathit{R}}){\mathit{s}}_{\hat{\mathit{f}}}{\hat{\mathit{a}}}_{\mathit{y}}]$

where

(19) $[\hat{\mathit{a}}{]}^{\mathrm{LVLH}}={\left[{{\displaystyle \hat{\mathit{a}}}}_{\mathit{J}2}+{{\displaystyle \hat{\mathbf{f}}}}_{\mathit{a}}\right]}^{\mathrm{LVLH}}$

The orbit perturbed acceleration term, â, is different form the term used in the truth model in which it does not contain the unmodeled disturbance acceleration term f. This navigation target model is used only to assist in the process of estimation. The dynamic modelfor the relative navigation states are:

(20) $\left[\begin{array}{c}\hat{\mathit{\rho }}(\mathit{t})\hfill \\ \widehat{\dot{\mathit{\rho }}}(\mathit{t})\hfill \end{array}\right]={\mathit{\phi }}_{\mathrm{LTV}}(\mathit{t},{\mathit{t}}_0)\left[\begin{array}{c}\hat{\mathit{\rho }}({\mathit{t}}_0)\hfill \\ \widehat{\dot{\mathit{\rho }}}({\mathit{t}}_0)\hfill \end{array}\right]$

where φ_LTV is the state transition matrix, and it is defi ned by Eq. (17) for the relative linear time varying model.

The navigation model for the target angular motion is used only to produce a reference attitude trajectory. This trajectory will be tracked by the chaser attitude control system.

(21a) ${\widehat{\dot{\mathit{q}}}}_{\mathit{l}}^{\mathit{t}}=\frac{1}{2}{\hat{\mathit{\omega }}}^{\mathit{t}}\oplus {\hat{\mathit{q}}}_{\mathit{l}}^{\mathit{t}}$

(21b) ${\widehat{\dot{\mathit{\omega }}}}^{\mathit{t}}={\hat{\mathit{I}}}_{\mathit{t}}^{-1}\left[{\hat{\mathit{\tau }}}_{\mathit{t}}-{\hat{\mathit{\omega }}}^{\mathit{t}}\times {\hat{\mathit{I}}}_{\mathit{t}}{\hat{\mathit{\omega }}}^{\mathit{t}}\right]$

(21c) ${\hat{\mathit{\tau }}}_{\mathit{t}}={\hat{\mathit{\tau }}}_{\mathit{tg}}$

For the chaser vehicle, the propagation of the state can be accomplished by using numerical integration techniques. However, in general, the gyros observations are sampled Here, the state vector x can represent relative position, at a high rate (usually higher than or at least equal to the same rate as the vector attitude observations). A discrete propagation is usually sufficient. Discrete propagation can be derived using a power series approach (Crassidis and Junkins, 2004).

(22) ${\hat{\mathit{q}}}_{\mathit{l}}^{\mathit{c}}\left(\mathit{t}\right)=\Omega \left({\hat{\mathit{\omega }}}^{\mathit{c}}\right){\hat{\mathit{q}}}_{\mathit{l}}^{\mathit{c}}\left({\mathit{t}}_0\right)$

where

(23a) ${\hat{\mathit{\omega }}}^{\mathit{c}}={\tilde{\mathit{\omega }}}^{\mathit{c}}-{\hat{\mathit{b}}}_{\mathit{\omega }}^{\mathit{c}}$

(23b) $\begin{array}{c}{{\displaystyle \widehat{\omega }}}^{\mathit{c}}={{\displaystyle \tilde{\omega }}}^{\mathit{c}}-{{\displaystyle \widehat{\mathit{b}}}}_{\omega }^{\mathit{c}}\hfill \\ \Omega \left({\hat{\mathit{\omega }}}^{\mathit{c}}\right)=\left[\begin{array}{cc}\cos \left(\frac{1}{2}‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}\right){\mathit{I}}_{3\times 3}-\left[{\hat{\mathit{\psi }}}^{\mathit{c}}\times \right]& {\hat{\mathit{\psi }}}^{\mathit{c}}\\ -{\hat{\mathit{\psi }}}^{{\mathit{c}}^{\mathit{T}}}& \cos \left(\frac{1}{2}‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}\right)\end{array}\right],\hfill \\ {\hat{\mathit{\psi }}}^{\mathit{c}}=\frac{\sin \left(\frac{1}{2}‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}\right){\hat{\mathit{\omega }}}^{\mathit{c}}}{‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖}\hfill \end{array}$

The propagation dynamic model for the error parameters is given by

(24) $\left\{\begin{array}{c}{\hat{\mathit{b}}}_{\mathit{\omega }}^{\mathit{c}}(\mathit{t})\\ {\hat{\mathit{\varepsilon }}}_{\mathit{S}}^{\mathit{S}}(\mathit{t})\\ {\hat{\mathit{\varepsilon }}}_{\mathit{l}}^{\mathit{l}}(\mathit{t})\end{array}\right\}={\mathit{\phi }}_{\mathrm{Markov}}(\mathit{t},{\mathit{t}}_0)\left\{\begin{array}{c}{\hat{\mathit{b}}}_{\mathit{\omega }}^{\mathit{c}}({\mathit{t}}_0)\\ {\hat{\mathit{\varepsilon }}}_{\mathit{S}}^{\mathit{S}}({\mathit{t}}_0)\\ {\hat{\mathit{\varepsilon }}}_{\mathit{l}}^{\mathit{l}}({\mathit{t}}_0)\end{array}\right\}$

where ϕ_Markovis defined as follows:

(25) ${\mathit{\phi }}_{\mathit{Marko}t\mathit{v}}=\left[\begin{array}{ccc}\mathit{e}-\frac{\Delta t}{{\tau }_{\mathit{b}}^{\omega }}{\mathit{I}}_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& \mathit{e}-\frac{\Delta t}{{\tau }^{\mathit{S}}}{\mathit{I}}_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& \mathit{e}\frac{\Delta t}{{\tau }^{\mathit{l}}}{\mathit{I}}_{3\times 3}\end{array}\right]$

An extended Kalman filter is derived from the nonlinear models as illustrated in the equations below (Brown and Hawag, 1997).

(26a) $\dot{\mathit{x}}=\mathit{f}(\mathit{x},\mathit{u},\mathit{t})+\mathit{w}(\mathit{t}\mathrm{),}\mathit{w}\left(\mathit{t}\right)\backsim \mathit{N}\left(0,\mathit{Q}\right)$

(26b) ${\tilde{\mathit{z}}}_{\mathit{k}}=\mathit{h}({\mathit{x}}_{\mathit{k}},\mathit{t})+{\mathit{v}}_{\mathit{k}},{\mathit{v}}_{\mathit{k}}\sim \mathit{N}(0,{\mathit{R}}_{\mathit{k}})$

Here, the state vector ẋ can represent relative position, velocity, and orientations of the chaser as well as other parameters that need to be estimated for the use by other flight algorithms. The time derivatives of the states x are a function of the states, inputs, time, and additive process noise w. This process noise is used to approximate unmodeled disturbances and other random disturbances to the dynamics. The measurements _~zk are modeled as a function of the states, time, and measurement noise v_k. The process noise and measurement noise are normally distributed with zero mean and covariance Q and R_k respectively.

The following steps summarize the Kalman filter equations, that are used to estimate the relative motion states and it is based on minimizing mean square of the error.

• Enter prior estimate of x-_k and its error covariance P-_k and compute the Kalman gain

  (27a) ${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}}^{-}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}{\left({\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}}^{-}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}\right)}^{-1}$

• Update estimate by measurement z~_k

  (27b) ${\hat{\mathit{z}}}_{\mathit{k}}^{-}=\mathit{h}\left({\hat{\mathit{x}}}_{\mathit{k}}^{-}\right)$
  (27c) ${\hat{\mathit{x}}}_{\mathit{k}}={\hat{\mathit{x}}}_{\mathit{k}}^{-}+{\mathit{K}}_{\mathit{k}}\left({\tilde{\mathit{z}}}_{\mathit{k}}-{\hat{\mathit{z}}}_{\mathit{k}}^{-}\right)$

• Compute error covariance for updated estimate

  (27d) ${\mathit{P}}_{\mathit{k}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}\right){\mathit{P}}_{\mathit{k}}^{-}$

• Project ahead

  (27e) ${\hat{\mathit{\varkappa }}}_{\mathit{k}+1}^{-}={\mathit{\phi }}_{\mathit{k}}{\hat{\mathit{\varkappa }}}_{\mathit{k}}$
  (27f) ${\mathit{P}}_{\mathit{k}+1}^{-}={\mathit{\phi }}_{\mathit{k}}{\mathit{P}}_{\mathit{k}}{\mathit{\phi }}_{\mathit{k}}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}}$

The term ϕ_k is the state transition matrix, and H_k is the measurement partial matrix that represents the sensitivity of the measurements to changes in the states. The state vector of the Kalman filter is defined to be:

(28) $\mathit{x}={\left[\mathit{\rho }\dot{\mathit{\rho }}{\mathit{\theta }}_{\mathit{c}}{\mathit{b}}_{\mathit{\omega }}^{\mathit{c}}{\mathbf{\varepsilon }}_{\mathit{S}}^{\mathit{S}}{\mathbf{\varepsilon }}_{\mathit{l}}^{\mathit{l}}\right]}^T$

and Kalman filter matrices are given by:

(29a) ${\mathit{\phi }}_{\mathit{k}}=\left[\begin{array}{ccc}{\mathit{\phi }}_{\mathit{LTV}}& {\mathbf{0}}_{6\times 6}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathit{\phi }}_{\mathit{Attitude}}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathbf{0}}_{6\times 6}& {\mathit{\phi }}_{\mathit{Markov}}\end{array}\right]$

(29b) ${\mathit{P}}_{\mathit{k}}=\left\{\begin{array}{ccc}{\mathit{P}}_{\mathrm{LTV}}& 0_{6\times 6}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathit{P}}_{\mathit{A}\mathrm{ttitude}}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathbf{0}}_{6\times 6}& {\mathit{P}}_{\mathrm{Markov}}\end{array}\right\}$

(29c) ${\mathit{Q}}_{\mathit{k}}=\left\{\begin{array}{ccc}{\mathit{Q}}_{\mathrm{LTV}}& {\mathbf{0}}_{6\times 6}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathit{Q}}_{\mathit{A}\mathrm{ttitude}}& {\mathbf{0}}_{6\times 6}\\ {\mathbf{0}}_{6\times 6}& {\mathbf{0}}_{6\times 6}& {\mathit{Q}}_{\mathrm{Markov}}\end{array}\right\}$

The state vector contains θ_c instead of q^c_I because the quaternion must obey a normalization constraint, which can be violated by the linear measurement updates associated with the filter. The most common approach to overcome this shortfall involves using a multiplicative error quaternion, where after neglecting higher order terms, the four component quaternion can effectively be replaced by a three component error vector θ_c (Crassidis and Junkins, 2004).Therefore, within first order, the quaternion update is given by:

(30) ${\mathit{q}}_{\mathit{T}}^{\mathit{c}+}=\mathbf{\delta }\mathit{q}({\mathit{\theta }}_{\mathit{c}})\otimes {\mathit{q}}_{\mathit{l}}^{C^{-}}$

and the discrete attitude error state transition matrix can also be derived using a power series approach to be:

(31) ${\mathit{\phi }}_{\mathrm{Attitude}}=\left\{\begin{array}{cc}{\mathit{\phi }}_{11}^{\mathit{c}}& {\mathit{\phi }}_{12}^{\mathit{c}}\\ {\mathit{\phi }}_{2l}^{\mathit{c}}& {\mathit{\phi }}_{22}^{\mathit{c}}\end{array}\right\}$

where

(32a) ${\mathit{P}}_{\mathit{LTV}}^{-}=\left[\begin{array}{cccccc}{\mathit{\sigma }}_{\mathit{x}}^2& 0& 0& 0& {\epsilon \mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\dot{\mathit{y}}}& 0\\ 0& {\mathit{\sigma }}_{\mathit{y}}^2& 0& {\epsilon \mathit{\sigma }}_{\dot{\mathit{x}}}{\mathit{\sigma }}_{\mathit{y}}& 0& 0\\ 0& 0& {\mathit{\sigma }}_{\mathit{z}}^2& 0& 0& 0\\ 0& {\epsilon \mathit{\sigma }}_{\dot{\mathit{x}}}{\mathit{\sigma }}_{\mathit{y}}& 0& {\mathit{\sigma }}_{\dot{\mathit{x}}}^2& 0& 0\\ {\epsilon \mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\dot{\mathit{y}}}& 0& 0& 0& {\mathit{\sigma }}_{\dot{\mathit{y}}}^2& 0\\ 0& 0& 0& 0& 0& {\mathit{\sigma }}_{\dot{\mathit{z}}}^2\end{array}\right]$

(32b) $\begin{array}{c}{\mathit{\phi }}_{12}^{\mathit{c}}=\left[{\hat{\mathit{\omega }}}^{\mathit{c}}\times \right]\frac{\left\{1-\cos \left(‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}\right)\right\}}{{‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖}^2}-{\mathit{I}}_{3\times 3}\Delta \mathit{t}-\hfill \\ {\left[{\hat{\mathit{\omega }}}^{\mathit{c}}\times \right]}^2\frac{\left\{‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}-\sin \left(‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖\Delta \mathit{t}\right)\right\}}{{‖{\hat{\mathit{\omega }}}^{\mathit{c}}‖}^3}\hfill \end{array}$

(32c) ${\mathit{\phi }}_{21}^{\mathit{c}}={\mathbf{0}}_{3\times 3}$

(32d) ${\mathit{\phi }}_{22}^{\mathit{c}}={\mathit{I}}_{3\times 3}$

By following the line steps of Woffinden and Geller (2007), Woffinden (2004) and Lear (1985), the initial error covariance matrix P_o, which represents how accurate the initial states are known, is given below for the proposed linear relative model, attitude, and error parameters.

(33a) ${\mathit{P}}_{\mathit{LTV}}^{-}=\left[\begin{array}{cccccc}{\sigma }_{\mathit{x}}^2& 0& 0& 0& {\epsilon \mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\dot{\mathit{y}}}& 0\\ 0& {\sigma }_{\mathit{y}}^2& 0& {\epsilon \mathit{\sigma }}_{\dot{\mathit{x}}}{\mathit{\sigma }}_{\mathit{y}}& 0& 0\\ 0& 0& {\sigma }_{\mathit{z}}^2& 0& 0& 0\\ 0& {\epsilon \mathit{\sigma }}_{\dot{\mathit{x}}}{\mathit{\sigma }}_{\mathit{y}}& 0& {\sigma }_{\dot{\mathit{x}}}^2& 0& 0\\ {\epsilon \mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\dot{\mathit{y}}}& 0& 0& 0& {\sigma }_{\dot{\mathit{y}}}^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{\dot{\mathit{z}}}^2\end{array}\right]$

(33b) ${\mathit{P}}_{\mathrm{Attitude}}^{-}=\left\{\begin{array}{cc}{\sigma }_{{\mathit{w}}_{\mathit{b}}^{\theta }}^2& 0_{3\times 3}\\ 0_{3\times 3}& {\sigma }_{{\mathit{w}}_{\mathit{b}}^{\omega }}^2\end{array}\right\}$

(33c) ${\mathit{P}}_{\mathrm{Markov}}^{-}=\left\{\begin{array}{cc}\sigma _{{\mathit{w}}^{\mathit{s}}}^2{}_{{\mathit{I}}_{3\times 3}}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& \sigma _{{\mathit{w}}^{\mathit{l}}}^2{}_{{\mathit{I}}_{\mathbf{3}\times \mathbf{3}}}\end{array}\right\}$

Parameters σ_x, σ_y and σ_z denote the standard deviation uncertainties of the relative position components, and σ_x, σy and σ_z are for the relative velocity components. Th e coeffi cient ε refers to the uncertainty correlation coupling between relative position and velocity components in the LVLH coordinate frame, and it ranges between a positive and a negative one. Th e standard deviations σ_wb^θ, σ_wb^ω, σ_w^S and σ_w^l are referring to the uncertainties of initial attitude, gyro biases, star tracker misalignments, and LIDAR misalignments, respectively. Th e discrete process noise matrix components of the relative motion canbe approximated by:

(34a) ${\mathit{Q}}_{\mathit{LTV}}=\left[\begin{array}{cccccc}{\mathit{\sigma }}_{{\mathit{w}}_{\mathit{x}}}^2\left(\frac{\Delta {\mathit{t}}^3}{3}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{x}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)& 0& 0\\ 0& {\mathit{\sigma }}_{\mathit{wy}}^2\left(\frac{\Delta {\mathit{t}}^3}{3}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{y}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)& 0\\ 0& 0& {\mathit{\sigma }}_{\mathit{wz}}^2\left(\frac{\Delta {\mathit{t}}^3}{3}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{z}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)\\ {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{x}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{x}}}^2\left(\Delta \mathit{t}\right)& 0& 0\\ 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{y}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{y}}}^2\left(\Delta \mathit{t}\right)& 0\\ 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{z}}}^2\left(\frac{\Delta {\mathit{t}}^2}{2}\right)& 0& 0& {\mathit{\sigma }}_{{\mathit{w}}_{\mathit{z}}}^2\left(\Delta \mathit{t}\right)\end{array}\right]$

(34b) ${\mathit{Q}}_{\mathit{Attitude}}=\left\{\begin{array}{cc}\left({\mathit{\sigma }}_{{\mathit{v}}_{\mathit{\omega }}^{\mathit{c}}}^2\Delta \mathit{t}+\frac{1}{3}{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{b}}^{\mathit{\omega }}}^2\Delta {\mathit{t}}^3\right){\mathit{I}}_{3\times 3}& -\left(\frac{1}{2}{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{b}}^{\mathit{\omega }}}^2\Delta {\mathit{T}}^2\right){\mathit{I}}_{3\times 3}\\ -\left(\frac{1}{2}{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{b}}^{\mathit{\omega }}}^2\Delta {\mathit{t}}^2\right){\mathit{I}}_{3\times 3}& -\left({\mathit{\sigma }}_{{\mathit{v}}_{\mathit{b}}^{\mathit{\omega }}}^2\Delta \mathit{t}\right){\mathit{I}}_{3\times 3}\end{array}\right\}$

(34c) ${\mathit{Q}}_{\mathit{Markov}}=\left[\begin{array}{cc}{\mathit{\sigma }}_{{\mathit{v}}^{\mathit{s}}}^2\Delta {\mathit{tI}}_{3\times 3}& 0_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathit{\sigma }}_{{\mathit{v}}^{\mathit{l}}}^2\Delta {\mathit{TI}}_{3\times 3}\end{array}\right]$

Here, σ_Wx, σ_Wy and σ_Wx are the standard deviations for the random unmodeled acceleration disturbances that act on the relative motion during the sample time period ∆t and σ_Vc, σ_Vbω, σ_Vs and σ_Vl are the random process uncertainty noises for gyros, gyro biases, star tracker misalignments, and LIDAR misalignments, respectively.

The measurements sensitivity matrices H_k and sensor measurements noise matrices R_k are defined for both star sensor and LIDAR as:

(35) ${\mathit{H}}_{\mathit{k}}=\left[\frac{{\mathit{H}}_{\mathit{k}}^{\mathit{lidar}}}{{\mathit{H}}_{\mathit{k}}^{\mathit{s}}}\right]$

(36) ${\mathit{R}}_{\mathit{k}}=\left[\begin{array}{cc}{\mathit{R}}_{\mathit{k}}^{\mathit{lidar}}\hfill & 0_{3\times 3}\hfill \\ 0_{3\times 3}\hfill & {\mathit{R}}_{\mathit{k}}^{\mathit{s}}\hfill \end{array}\right]$

The measurement partials for the azimuth, elevation and range measurements are computed with the help of the LIDAR measurement range vector. Utilizing Eq. (13) and small angle approximations leads to the following equation for the relative range in terms of the navigation states:

(37) ${\mathit{\rho }}^{\mathit{l}}=\rho \left(\frac{\begin{array}{cc}{\mathit{c}}_{\alpha }& {\mathit{c}}_{\beta }\\ {\mathit{s}}_{\alpha }& {\mathit{c}}_{\beta }\end{array}}{{\mathit{s}}_{\beta }}\right)=[{\mathit{I}}_{\mathit{3}\times \mathit{3}}-{\mathbf{\varepsilon }}_{\mathit{l}}^{{\mathit{l}}^X}\mathrm{][}{\mathit{I}}_{3\times 3}-{\mathbf{\varepsilon }}_{\mathit{s}}^{{\mathit{l}}^{\times }}]\lceil {\mathit{I}}_{3\times 3}-\mathbf{\theta }_{\mathit{c}}^{\times }\rceil {\mathit{T}}^{\widehat{\mathit{sI}}}({\mathbf{q}}_{\mathit{l}}^{{\mathit{c}}^{-}}){\mathit{T}}^{\mathit{IT}}(-\mathit{\rho })$

Using the chain rule, the partial of the range vector with respect to the navigation states can be expressed as (Woffinden and Geller, 2007)^:

(38) $\begin{array}{c}{\left.\frac{\partial {\rho }^{\mathit{l}}\left(\mathit{x}\right)}{\partial \mathit{x}}\right|}_{\widehat{\mathit{x}}}={\left[\frac{\partial {\rho }^{\mathit{l}}\left(\mathit{x}\right)}{\partial \alpha }\frac{\partial \alpha \left(\mathit{x}\right)}{\partial \mathit{x}}+\frac{\partial {\rho }^{\mathit{l}}\left(\mathit{x}\right)}{\partial \beta }\frac{\partial \beta \left(\mathit{x}\right)}{\partial \mathit{x}}+\frac{\partial {\rho }^{\mathit{l}}\left(\mathit{x}\right)}{\partial \rho }\frac{\partial \rho \left(\mathit{x}\right)}{\partial \mathit{x}}\right]}_{\widehat{\mathit{x}}}\hfill \\ ={\left[\left(\mathit{\rho }{\mathit{p}}_{\alpha }^{\mathit{l}}\right){\mathit{h}}_{\alpha }+\left(\mathit{\rho }{\mathit{p}}_{\beta }^{\mathit{l}}\right){\mathit{h}}_{\beta }+{\mathit{p}}_{\rho }^{\mathit{l}}{\mathit{h}}_{\rho }\right]}_{\widehat{\mathit{x}}}\hfill \end{array}$

(39) ${\mathit{p}}_{\alpha }^{\mathit{l}}=\left[\frac{\begin{array}{c}-{\mathit{c}}_{\mathit{\beta }}{\mathit{s}}_{\mathit{\alpha }}\\ {\mathit{c}}_{\mathit{\beta }}{\mathit{c}}_{\mathit{\alpha }}\end{array}}{0}\right],{\mathit{p}}_{\mathit{\beta }}^{\mathit{l}}=\left[\frac{\begin{array}{c}-{\mathit{s}}_{\mathit{\beta }}{\mathit{c}}_{\mathit{\alpha }}\\ -{\mathit{s}}_{\mathit{\beta }}{\mathit{s}}_{\mathit{\alpha }}\end{array}}{{\mathit{c}}_{\mathit{\beta }}}\right],{\mathit{p}}_{\mathit{\rho }}^{\mathit{l}}=\left[\frac{\begin{array}{c}{\mathit{c}}_{\beta }{\mathit{c}}_{\alpha }\\ {\mathit{c}}_{\beta }{\mathit{s}}_{\alpha }\end{array}}{{\mathit{s}}_{\beta }}\right]$

The measurement geometry can now be computed by taking the advantages of the property that p_α^l, p_β^l and p_p^l are orthogonal to each other and taking the dot product with respect to each of them.

(40) ${\mathit{h}}_{\mathit{a}}=\frac{({\mathit{p}}_{\chi }^{\mathit{l}}{)}^T}{\mathit{\rho }{\mathit{c}}_{\beta }^2}\left[\frac{\partial {\mathit{\rho }}^{\mathit{l}}(\mathit{x})}{\partial \mathit{x}}\right],{\mathit{h}}_{\mathit{\beta }}=\frac{({\mathit{p}}_{\mathit{\beta }}^{\mathit{l}}{)}^T}{\mathit{\rho }}\left[\frac{\partial {\mathit{\rho }}^{\mathit{l}}(\mathit{x})}{\partial \mathit{x}}\right],{\mathit{h}}_{\rho }=({\mathit{p}}_{\rho }^{\mathit{l}}{)}^T\left[\frac{\partial {\mathit{\rho }}^{\mathit{l}}(\mathit{x})}{\partial \mathit{x}}\right]$

The evaluation of the relative range vector with respect to the navigation states yields

(41) $\left[\frac{\partial {\mathit{\rho }}^{\mathit{l}}(\mathit{x})}{\partial \mathit{x}}\right]=[-{\mathit{T}}^{\mathit{lT}}{\mathbf{0}}_{\mathit{3}\times 3}-{\mathit{T}}^{\mathit{l}\hat{\mathit{l}}}{\mathit{T}}^{\widehat{\mathit{ls}}}[{\mathit{\rho }}^{\hat{\mathit{s}}}\times ]{\mathbf{0}}_{3\times \mathit{3}}-{\mathit{T}}^{\mathit{l}\hat{\mathit{l}}}[{\mathit{\rho }}^{\mathit{s}}\times ]-[{\mathit{\rho }}^{\hat{\mathit{l}}}\times \mathrm{]]}$

Now, the LIDAR measurement sensitivity matrix and covariance matrix can be written as:

(42) ${\mathit{H}}_{\mathit{k}}^{\mathit{lidar}}=[{\hat{\mathit{h}}}_{\mathit{\alpha }}{\hat{\mathit{h}}}_{\mathit{\beta }}{\hat{\mathit{h}}}_{\mathit{\rho }}{]}^T$

and

(43) ${\mathit{R}}_{\mathit{k}}^{\mathit{li}d\mathit{ar}}=\left[\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{\alpha }}^2\hfill & 0\hfill & 0\hfill \\ 0\hfill & {\mathit{\sigma }}_{\mathit{\beta }}^2\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\mathit{\sigma }}_{\mathit{\rho }}^2\hfill \end{array}\right]$

When processing star tracker data, a derived measurement is calculated (Woffinden and Geller, 2007). This quantity is effectively the residual to be processed by the filter.

(44) $\left[\frac{\frac{1}{\mathit{z}}{\tilde{\mathit{Z}}}_{\mathit{S}}}{1}\right]={\tilde{\mathit{q}}}_{\mathit{l}}^{\mathit{s}}\otimes \left[{\hat{\mathit{q}}}_{\mathit{c}}^{\mathit{I}}\otimes {\hat{\mathit{q}}}_{\mathit{s}}^{\mathit{c}}\otimes \mathit{\delta }\hat{\mathit{q}}\left({\hat{\mathit{\varepsilon }}}_{\mathit{s}}^{\mathit{s}}\right)\right]$

The derived star tracker measurement can be written as a function of the navigation states as:

(45) ${\tilde{\mathit{Z}}}_{\mathit{S}}={\mathit{h}}_{\mathit{s}}(\mathit{x})+{\mathit{v}}_{\mathit{S}}^{\mathit{S}}={\mathit{\theta }}_{\mathit{c}}^{\mathit{s}}+{\mathit{\varepsilon }}_{\mathit{s}}^{\mathit{s}}+{\mathit{v}}_{\mathit{S}}^{\mathit{S}}$

Therefore, the measurement sensitivity matrix for the star tracker can be derived to be

(46) ${{\displaystyle \left.{\mathit{H}}_{\mathit{k}}^{\mathit{s}}=\frac{\partial {\mathit{h}}_{\mathit{s}}(\mathit{x})}{\partial \mathit{x}}\right|}}_{\hat{\mathit{x}}}=[{\mathbf{0}}_{6\times 6}{\mathit{I}}_{3\times 3}{\mathbf{0}}_{3\times 3}{\mathit{I}}_{3\times 3}{\mathbf{0}}_{3\times 3}]$

and the star tracker measurement covariance is

(47) ${\mathit{R}}_{\mathit{k}}^{\mathit{s}}={\sigma }_{\mathit{s}}^2{\mathit{I}}_{3\times 3}$

For close proximity operations, a propositional-derivative (PD) controller is employed for both the rotational and translational controls. The commanded torques for the chaser spacecraft to match its orientation with the target vehicle are computed as

(48) ${\hat{\mathit{\tau }}}_{\mathit{cc}}={\mathit{K}}_{\mathit{q}}\mathbf{\delta }\hat{\mathit{q}}+{\mathit{K}}_{\omega }\mathbf{\delta }\hat{\omega }$

where

(49) $\delta \hat{\mathit{q}}=\delta {\hat{\mathit{q}}}_{{\mathit{e}}_{\mathit{v}}}\delta {\hat{\mathit{q}}}_{{\mathit{e}}_{\mathit{s}}}$

(49b) $\delta \hat{\omega }={\hat{\omega }}_{\mathit{des}}^{\mathit{c}}-{\hat{\omega }}^{\mathit{c}}$

and

(50a) ${{\displaystyle \widehat{\mathit{q}}}}_{{\mathit{I}}_{\mathit{des}}}^{\mathit{c}}={{\displaystyle \widehat{\mathit{q}}}}_{\mathit{I}}^{\mathit{t}}$

(50b) ${{\displaystyle \widehat{\omega }}}_{\mathit{des}}^{\mathit{c}}={\mathit{T}}^{\mathit{ct}}{{\displaystyle \widehat{\omega }}}^{\mathit{t}}$

(50c) ${\hat{\mathit{q}}}_{\mathit{e}}={\hat{\mathit{q}}}_{{\mathit{I}}_{\mathit{des}}}^{\mathit{c}}\otimes {\hat{\mathit{q}}}_{\mathit{I}}^{\mathit{c}}=\left(\frac{\partial {\hat{\mathit{q}}}_{{\mathit{e}}_{\mathit{v}}}}{\partial {\hat{\mathit{q}}}_{{\mathit{e}}_{\mathit{s}}}}\right)$

̂q^c_Idesc and ̂ω^c_desc are the desired orientation and angular velocity, respectively, to be tracked by the chaser vehicle. The angular offset and angular rate offset between target and chaser are denoted by δ̂q_e and δ̂ω, respectively. The proportional and derivative control gains K_q and K_ω are determined based on the desired natural frequency ω_θ, damping ratio ζ_θ of the attitude control system, and the moment of inertia of the chaser spacecraft I_c (Wie, 1998).

(51) ${\mathit{K}}_{\mathit{q}}={\omega }_{\theta }^2{\mathit{I}}_{\mathit{c}}{\mathit{K}}_{\omega }=2{\zeta }_{\theta }{\omega }_{\theta }{\mathit{I}}_{\mathit{c}}$

On the other hand, The translation control algorithm computes the required continuous thrust, f_c, based on the previous linear model, in order to track the desired trajectory specified by the following guidance algorithm:

(52a) ${\mathit{U}}_{\Delta \mathit{V}}={\mathbf{f}}_c={\mathit{K}}_{\mathit{p}}\mathit{\delta }\mathit{\rho }+{\mathit{K}}_{\dot{\mathit{\rho }}}\mathit{\delta }\dot{\mathit{\rho }}$

(52b) $\mathit{\delta }\mathit{\rho }={\mathit{\rho }}_{\mathit{des}}-\hat{\mathit{\rho }}$

(52c) $\mathbf{\delta }\dot{\mathbf{\rho }}=\dot{\mathit{\rho }}\mathit{des}-\widehat{\dot{\mathbf{\rho }}}$

The proportional and derivative control gains K_ρ and K_ρ̇. are determined based on the desired natural frequency ω_ρ and damping ratio ζ_ρ of the translational control system.

(53) ${\mathit{K}}_{\rho }={\omega }_{\mathit{\rho }}^2{\mathit{I}}_{3\times 3}{\mathit{K}}_{\dot{\rho }}=2{\zeta }_{\mathit{\rho }}{\omega }_{\mathit{\rho }}{\mathit{I}}_{3\times 3}$

Variables ρ_des and ρ̇_des are, respectively, the desired relative position and relative velocity to be tracked by the chaser vehicle, and it is defined by the guidance algorithms. It is worth noting that the equivalent continuous velocity increment ∆V, based on the continuous thrust, can be approximated for small to be

(54) $\Delta \mathit{V}\approx {\mathit{U}}_{\Delta \mathit{V}}\Delta \mathit{t}$

SIMULATION EXAMPLES

The key metrics of the analysis fall into three main categories. The first is navigation performance, which is how well the states are estimated by the filter. This metric is measured by the navigation error, the difference between the true states and the filter states. The second is trajectory control performance, which is a measure of how closely the chaser vehicle is able to follow the guidance algorithms. The third is fuel performance, or ∆V fuel usage, and it is computed based on the linear model developed in the previous section.

The preceding guidance and navigation algorithms are illustrated now through different examples. Initial conditions for simulation are listed in Tables 4 to 6.

Table 4
Navigation Filter Parameters.

Parameter	Value
Initial Relative Position and Velocity Uncertainties	σx = σy = σz = 33.3 m σẋ = σẏ = σż = 0.01 m/s
Process Noise	σwx = σwy = 5 × 10-4m/s3/2 σwz = 5 × 10-6m/s3/2
Measurements Noise	σα = σβ = 0.06 deg, σρ = 0.5m
Simulation Step	0.1 s
Measurements Update	1 Hz

Table 5
Vehicles Orbital Elements.

Paramete	Target	Chaser
a,km	6723.2576	6723.2576
e	0.1	0.1
i, deg	51.6467	51.6467
Ω,deg	188. 0147	188. 0147
ω,deg	174.3022	174.3022
f,deg	270.0882	270.0832

Table 6
Simulation Initial Conditions.

Initial Conditions		Parameter	Value
Inertia	Chaser	Ic	$\left\{\begin{array}{ccc}1.31& 0& 0\\ 0& 10.220& 0\\ 0& 0& 10.65\end{array}\right\}\times {10}^6\mathrm{kg}-m^2$
	Target	It	$\left\{\begin{array}{ccc}128& 0& 0\\ 0& 107& 0\\ 0& 0& 201\end{array}\right\}\times {10}^6\mathrm{kg}-m^2$
Initial Relative Attitude Errors		δϕ,δθ,δψ	[(7.5&-7.5&7.5)]deg
Control Parameters	Rotationalnatural frequency	ωθ	1/30 s-1
	Rotational damping ratio	ζθ	0.7
	Translational natural frequency	ωρ	1/50 s-1
	Translational damping ratio	ζρ	0.7
Unmodeled Disturbances	Rotational disturbances	τId	kg-km2/s2
	Translational disturbances	fw	km/s2
Sensors Errors	Gyro error (3)	Drift rate	3 deg/hr/axis
		Random walk	0.05 mrad/s1/2
	Star Tracker error (3)	Misalignment	1 mrad/axis
		Noise	1 mrad/axis
		Measurements	1 Hz
	LIDAR error (3)	Misalignment	1 mrad/axis
		Noise	[1 mrad 1 mrad 0.5 m]
		Measurements	1 Hz

A Simulink model is built using the MATLAB software to demonstrate the closed-loop guidance transfer of the chaser in order to approach and to depart from the target vehicle in any orbit, either circular or elliptic, given uncertain initial conditions, noisy measurements, and limited dynamics. This model consists of three main parts, guidance, navigation, and control. The proposed linear time varying model is used in designing the navigation filter and in maneuver targeting of the guidance system. The target is assumed to be in a passive nadir pointing mode andnot in maneuvering. The chaser uses star tracker data and gyro data to determine attitude and attitude rate. Momentum wheels and PD controller are used to point the chaser LIDAR at the target. The chaser uses LIDAR data to determine the relative position and velocity of the target. Maneuver targeting algorithms, based on PD controller, are used to compute commands in the chaser body frame as to track the desired trajectory.

The performance of the navigation system is shown in Figs. 5 to 7. In this case, the thrusters are off and both target and chaser vehicles are initially in the same neighborhood (Table 5). Figure 5 shows the relative position and relative velocity between the vehicles during simulation. Figure 6 depicts how accurately the navigation system can estimate the chaser's relative position and velocity. Form this figure, the filter is able to converge within few minutes and the relative position and velocity can be accurately estimatedwithin the accuracy of the sensors. The attitude navigation errors and the PD control tracking performance are shown in Fig. 7. As indicated by this figure, the chaser attitude navigation system is able to converge quickly and the chaser attitude PD controller can track the target attitude and angular velocity trajectories.

Figure 5
Relative Motion Without.

Figure 6
Relative Motion Navigation Performance Without.

Figure 7
Chaser Attitude Navigation and Control Performance.

The basic glidelope rendezvous and close proximity operations scenario used to evaluate the performance of the entire closedloop relative position and attitude control system with the navigation filter consists of two main segments: the inbound and theoutbound segments. Each segment of the glideslope is followed by 3 minutes of station keeping. First, the inbound segment: the chaser starts to approach the target form [58-580 0] m behind the target and ends at [0-100 0] m. After 3 minutes of station keeping at -100 m behind the target, the chaser starts to depart away from the target and leading to a new location -1000 m behind the target. The chase then stays at rest at that location for another 3 minutes. The results of this scenario are shown in Figs. 8 and 9. In all of these figures, different segments of the glideslope are shown, and the variations of in-plane relative motion of the chaser with respect to target vehicle are presented. Figure 8 shows the relative position and velocity plots of relative motion along with the required in order to achieve this trajectory maneuver, while Fig. 9 shows the error in relative position and velocity between the truth model and the navigation model. In all of the above glideslopes, the overall performance of the rendezvous and proximity operations are satisfactory.

Figure 8
Relative Motion Inbound-Outbound Glideslope (Summary Scenario).

Figure 9
Relative Motion Navigation and Control Performance (Summary Scenario).

The continuous thrust is calculated using the estimated relative position and velocity, either from the Kalman filter or from the knowledge of initial conditions, not from the true relative position and velocity of the chaser. As such, the chaser is not expected to reach its intended place exactly, but in the neighborhood thereof. Aided by the sensors, the initial estimation errors subside to an optimal level, determined by the ratio of the process noise matrix , and the measurement noise matrix ,earlier defined. Because of the active range and the angle measurements from the LIDAR system, and relatively small measurement errors, the true and the estimated relative position and velocity states are almost indistinguishable, as seen in previous figures during the steady state.

CONCLUSION

The results of this study indicate that the proposed linear model is clearly effective at estimating the relative position and velocity and controlling the relative trajectory. In addition, this model is not restricted to a circular orbit but it can be used as well for an eccentric orbit. Furthermore, by using this model, simple guidance algorithms for glideslope are developed to autonomously approach and depart form a target vehicle. The relative navigation in this study is utilizing range, azimuth, and elevation measurements of the target relative to the chaser froma simulated LIDAR system, along with the star tracker and gyro measurements of the chaser and an extended Kalman filter. The vehicle attitude dynamics, attitude tracking control, attitude determination, and uncertainties like measurement biases and sensor misalignments are considered in this study to fire the thrusters in the right direction and spin the momentum wheels at the proper rate in the chaser coordinate frame. The analyst must consider, in addition, off nominal situations, limitations and operational range of the sensors, and limitations of the actuators. These topics and others will be addressed in the future.

REFERENCES

• Brown, R.G. and Hawag, P., 1997, "Introduction to Random Signals and Applied Kalman Filtering", 3rd Edition, John Wiley & Son Inc., United States.
• Broucke, R.A., 2003, "Solution of the Elliptic Rendezvous Problem with the Time as Independent Variable", Journal of Guidance, Control, and Dynamics,Vol. 26, No. 4, pp. 615-621. doi: 10.2514/2.5089.
  » https://doi.org/10.2514/2.5089
• Carter,T.E., 1998, "State Transition Matrices for Terminal Rendezvous Studies: Brief Survey and New Example", Journal of Guidance, Control, and Dynamics, Vol. 21, No. 1, pp. 148-155.doi: 10.2514/2.4211.
  » https://doi.org/10.2514/2.4211
• Cho, H.C. and Park, S.Y., 2009, "Analytical Solution for Fuel-Optimal Reconfiguration in Relative Motion",Journal of Optimization Theory and Applications, Vol. 141, No.3, pp. 495-512.doi: 10.1007/s10957008-9482-3.
  » https://doi.org/10.1007/s10957008-9482-3
• Clohessy, W.H. and Wiltshire, R.S., 1960, "Terminal Guidance System for Satellite Rendezvous", Journal of the Aerospace Sciences, Vol. 27, No. 9, pp. 653-678.
• Crassidis, J.L. and Junkins, J.L., 2004, "Optimal Estimation of Dynamic System", 1st Edition, CRC Press LLC, United States.
• Fehse, W., 2003, "Automated Rendezvous and Docking of Spacecraft", 1st Edition, Cambridge University Press, United Kingdom.
• Inalhan, G., Tillerson, M., and How, J. P., 2002, "Relative Dynamics and Control of Spacecraft Formation in Eccentric Orbits", Journal of Guidance, Control, and Dynamics,Vol. 25, No. 1, pp. 48-58.
• Jenkins, S.C. and Geller, D.K., 2007, "State Estimation and Targeting For Autonomous Rendezvous and Proximity Operations", AAS 07-316, Proceedings of theAIAA/AASAstrodynamics Specialists Conference, Mackinac Island, MI.
• Junkins, J.L., Kim, S., Crassidis, L., Cheng, Y., and Fosbury, A.M., 2005,"Kalman Filtering for Relative Spacecraft Attitude and Position Estimation", AIAA 2005-6087, Proceedings of the AIAA Guidance, Navigation, and Control Conference, San Francisco, California.
• Lear, W.M., 1985, "Kalman Filtering Techniques", NASA Johnson Space Center: Mission Planning andAnalysis Division, Houston, TX, JSC-20688.
• Melton, R.G., 2000, "Time-Explicit Representation of Relative Motion Between Elliptica lOrbits" Journal of Guidance, Control and Dynamics, Vol. 23, No. 4, pp. 604-610.doi: 10.2514/2.4605.
  » https://doi.org/10.2514/2.4605
• Sengupta, P. and Vadali, S.R., 2007, "Relative Motion and the Geometry of Formations in KeplerianElliptic Orbits with Arbitrary Eccentricity", Journal of Guidance, Control, and Dynamics,Vol. 30, No. 4, pp. 953-964.doi: 10.2514/1.25941.
  » https://doi.org/10.2514/1.25941
• Schaub, H. and Junkins, J.L., 2003, "Analytical Mechanics of Space Systems", American Institute of Aeronautics and Astronautics, Reston, Virginia, United States.
• Tschauner, J. and Hempel, P., 1965, "Rendezvous zueinem in Elliptischer Bahn Umlaufenden Ziel", ActaAstronautica, Vol. 11, pp. 104-109.
• Vallado, D.A., 2001, "Fundamentals of Astrodynamics and Applications", 2nd Edition, Microcosm Press, El Segundo, California.
• Woffinden, D.C. and Geller, D.K., 2007, "Navigating the Road to Autonomous Orbital Rendezvous", Journal of Spacecraft and Rockets, Vol. 44, No. 4, 2007, pp. 898-909.doi: 10.2514/1.30734.
  » https://doi.org/10.2514/1.30734
• Woffinden, D.C. and Geller, D.K., 2007, "Relative Angles-Only Navigation and Pose Estimation for Autonomous Orbital Rendezvous", Journal of Guidance, Control and Dynamics, Vol. 30, No.5, pp. 1455-1469.doi: 10.2514/1.28216.
  » https://doi.org/10.2514/1.28216
• Woffinden, D.C., 2004, "On-Orbit Satellite Inspection Navigation and Analysis", M.S. Thesis, Aeronautics and Astronautics Department,Massachusetts Institute of Technology, Cambridge, MA, United States.
• Wie, B., 1998,"Space Vehicle Dynamics and Control", American Institute of Aeronautics and Astronautics, Reston, Virginia, United States.
• Yamanaka,K. and Ankersen,F., 2002, "New State Transition Matrix for Relative Motion on an Arbitrary Elliptical Orbit"Journal of Guidance, Control, and Dynamics, Vol. 25, No. 1, pp. 60-66.

Publication Dates

• Publication in this collection
  Jul-Sep 2014

History

• Received
  28 Jan 2014
• Accepted
  29 July 2014