3.1. System Dynamics

During the Mars orbit phase, the probe dynamics should be established in the J2000·0 Mars-centre Mars centre equator inertial coordinate system. The dynamic model can be written as

(1)$$ \ddot{\bi r} = {\bi a}_0 + {\bi a}_{ns} + {\bi a}_{nbody} + {\bi a}_{srp} + {\bi a}_{drag} + {\bi a}_\varepsilon $$

For Mars probe, besides the Mars body centre gravitational force a₀, other perturbation forces include the Mars non-spherical gravitational force a_ns, the third body gravitational perturbations from the Sun and other major planets a_nbody, the solar radiation pressure a_srp. Atmosphere drag a_drag should also be considered for very low orbit, which is given in the Appendix. In addition, a_ε is other much smaller force acting on the spacecraft that can be ignored, such as gravitational forces of Martian moons. The RKF7(8) algorithm is applied for the orbit propagation of this dynamic model.

Because the on board computer has limited capacity, the on board dynamic navigation model must be simplified. In order to analyse how the model errors affect navigation precision, a precise dynamic model (model-1) is set as a criterion and a simplified model (model-2) is used for state estimation for on board navigation, as shown in Table 2.

Table 2.

Different options of precise dynamics model and simplified model.

In order to analyse the propagating errors of the simplified model, there are four different orbits being propagated for two days, which are LMO with low-inclination (a = 3996 km, e = 0·001, i = 0°), LMO with high-inclination (a = 3696 km, e = 0·001, i = 75°), HMO with low-inclination (a = 20430 km, e = 0·001, i = 0°) and HMO with high-inclination (a = 20000 km, e = 0·001, i = 75°).

The position difference between model-1 and model-2 is compared in Figure 4 and Figure 5. The figures show that the position errors due to dynamics model error of HMO are several tens of metres for two days' propagation. As to LMO, the position error for two days' propagation is several tens of kilometres.

Figure 4.

HMO propagation difference between simplified model and precise model.

Figure 5.

LMO propagation difference between simplified model and precise model.

3.2. Measurement models.

3.2.1. Optical measurement

An image photo of the probe on board optical sensor is sketched in Figure 6. The angle θ = arcos(n_s·n_p) is the angle between n_s and n_p, where, n_s is the line-of-sight vector from the probe to the background star in the probe centre inertial coordinate system, and n_p is that of the Mars moon. These vectors are in the probe centre inertial coordinate system.

Figure 6.

Optical measurement of Martian moon and background star.

If there are three background stars and a Martian moon in image, then the image should have three angles (θ ₁,θ ₂,θ ₃) between the Martian moon and the three background stars. The position unit vectors of the three background stars (n_s1, n_s2, n _s3) can be obtained from the Tycho-2 star catalogue. According to the angles (θ ₁,θ ₂,θ ₃) and the unit vectors (n_s1, n_s2, n _s3), the unit vector n_p is obtained as

(2)$$\left\{\matrix{\theta _1 = \arccos ({\bi n}_{\rm p} \cdot {\bi n}_{{\rm s1}}) \cr \theta_2 = \arccos ({\bi n}_{\rm p} \cdot {\bi n}_{{\rm s2}}) \cr \theta _3 = \arccos ({\bi n}_{\rm p} \cdot {\bi n}_{{\rm s3}} )} \right. $$

If there are more than three background stars in image, the least square algorithm can be used. In the probe centre inertial coordinate system, n_p(x _np, y _np, z _np) is the unit vector of the Martian moon, the measurement data right ascension α and declination δ can be obtained.

(3)$$\alpha = \arctan \left(\displaystyle{{y_{{\rm np}}} \over {x_{{\rm np}}}} \right),\quad\delta = \arcsin \left(\displaystyle{{z_{{\rm np}}} \over {r_{{\rm np}}}} \right)$$

In the Mars centre inertial coordinate system, r(x, y, z) is the position vector of probe, and r_m(x _m, y _m, z _m) is the Martian moon position vector. Then, r_p(x _p, y _p, z _p) = r_m – r is the position vector of the Martian moon in the probe centre inertial coordinate system, the measurement model right ascension α and declination δ are obtained:

(4)$$\alpha = \arctan \left(\displaystyle{{y_{\rm p}} \over {x_{\rm p}}} \right),\quad\delta = \arcsin \left(\displaystyle{{z_{\rm p}} \over {r_{\rm p}}} \right)$$

The partial derivative of right ascension α and declination δ for state variables are as follows,

(5)$$\displaystyle{\partial \alpha \over \partial {\bi r}} = \displaystyle{1 \over x_{\rm p}^2 + y_{\rm p}^2} \left(\matrix{y \cr - x \cr 0} \right), \quad \displaystyle{\partial \alpha \over \partial \dot{\bi r}} = 0$$

(6)$$\displaystyle{{\partial \delta} \over {\partial {\bi r}}} = \displaystyle{1 \over {r_{\rm p}^2 \sqrt {x_{\rm p}^2 + y_{\rm p}^2}}} \left( \matrix{ - x_{\rm p} z \cr - y_{\rm p} z \cr x_{\rm p}^2 + y_{\rm p}^2 } \right),\quad\displaystyle{{\partial \delta} \over {\partial \dot{\bi r}}} = 0$$

3.2.2. Range measurement and Doppler measurement

There are two probes orbiting Mars, and the first probe position and velocity vectors are ${\bi r}_1 (x_1, y_1, z_1 )$, and $ \dot{\bi r}_1 (\dot x_1, \dot y_1, \dot z_1 )$, and those of the second probe are r₂(x ₂, y ₂, z ₂), and $ \dot{\bi r}_2 (\dot x_2, \dot y_2, \dot z_2 )$.

Then, the models of range measurement ρ and Doppler measurement $\dot \rho $ are expressed

(7)$$\rho = \Vert {{\bi r}_1 - {\bi r}_2} \Vert = \sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2 + (z_1 - z_2 )^2} $$

(8)$$\dot \rho = \displaystyle{{({\bi r}_1 - {\bi r}_2 ) \cdot ( \dot{\bi r}_1 - \dot{\bi r}_2 )} \over \rho} = \displaystyle{{(x_1 - x_2 )(\dot x_1 - \dot x_2 ) + (y_1 - y_2 )(\dot y_1 - \dot y_2 ) + (z_1 - z_2 )(\dot z_1 - \dot z_2 )} \over {\sqrt {(x_1 - x_2 )^2 + (y_1 - y_2 )^2 + (z_1 - z_2 )^2}}} $$

The partial derivative of range measurement ρ and Doppler measurement $\dot \rho $ for state variables are computed as follows,

(9)$$\displaystyle{{\partial \rho} \over {\partial {\bi r}_1}} = \displaystyle{1 \over \rho} \left( \matrix{x_1 - x_2 \cr y_1 - y_2 \cr z_1 - z_2 } \right),\quad\displaystyle{{\partial \rho} \over {\partial \dot{\bi r}_1}} = \left( \matrix{0 \cr 0 \cr 0} \right)$$

(10)$$\displaystyle{{\partial \rho} \over {\partial {\bi r}_2}} = \displaystyle{1 \over \rho} \left( \matrix{x_2 - x_1 \cr y_2 - y_1 \cr z_2 - z_1 } \right),\quad\displaystyle{{\partial \rho} \over {\partial \dot r_2}} = \left( \matrix{0 \cr 0 \cr 0 } \right)$$

(11)$$\displaystyle{{\partial \dot \rho} \over {\partial {\bi r}_1}} = \displaystyle{1 \over \rho} \left( \matrix{(\dot x_1 - \dot x_2 ) - \displaystyle{{\dot \rho} \over \rho} (x_1 - x_2 ) \cr (\dot y_1 - \dot y_2 ) - \displaystyle{{\dot \rho} \over \rho} (y_1 - y_2 ) \cr (\dot z_1 - \dot z_2 ) - \displaystyle{{\dot \rho} \over \rho} (z_1 - z_2 )} \right),\quad\displaystyle{{\partial \dot \rho} \over {\partial \dot{\bi r}_1}} = \displaystyle{1 \over \rho} \left( \matrix{x_1 - x_2 \cr y_1 - y_2 \cr z_1 - z_2} \right)$$

(12)$$\displaystyle{{\partial \dot \rho} \over {\partial {\bi r}_2}} = \displaystyle{1 \over \rho} \left( \matrix{(\dot x_2 - \dot x_1 ) - \displaystyle{{\dot \rho} \over \rho} (x_2 - x_1 ) \cr (\dot y_2 - \dot y_1 ) - \displaystyle{{\dot \rho} \over \rho} (y_2 - y_1 ) \cr (\dot z_2 - \dot z_1 ) - \displaystyle{{\dot \rho} \over \rho} (z_2 - z_1 )} \right),\quad\displaystyle{{\partial \dot \rho} \over {\partial \dot{\bi r}_2}} = \displaystyle{1 \over \rho} \left( \matrix{x_2 - x_1 \cr y_2 - y_1 \cr z_2 - z_1} \right)$$

3.3. Estimation methods

Generally, there are two types of estimation methods used in space navigation. One is the batch method, typically weighted least squares algorithms; the other is sequential processing algorithms including the Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF) and so on. Usually, the batch estimation method is not used for on board estimation because it requires that all measurements to be considered are available before a priori information can be improved. So, the sequential processing algorithm including EKF and UKF is selected to be applied and compared in this paper.

The positions and velocities of the two probes must be estimated altogether to fulfil the requirement of SST radio range and range-rate measurement.

So the state variables include the first probe's position vector r₁, velocity vector $\dot {\bi r}_1 $ and the second probe's position vector r₂ and velocity vector $\dot {\bi r}_2 $. The state equation and measurement equation are

(13)$$\left\{ \matrix{ \dot{\bi x}(t) = {\bi f}[{\bi x}(t),t] + {\bi w}(t) \cr {\bi y}(t) = {\bi h}[{\bi x}(t),t] + {\bi R}(t)} \right.$$

where ${\bi x}(r_1, \dot r_1, r_2, \dot r_2 )$ is the state parameter, y is measurement consisting of (α, δ) or $(\rho, \dot \rho )$, w is the model noise, and R is the measurement noise.

At time t _k, x_k is the a priori state value, and P_k is the a priori covariance matrix. If at time t _k+1, new measurement data y_k+1 are obtained, we calculate the estimated value x_k+1 at t _k+1. After linearization, Equation (13) is written as

(14)$$\left\{ \matrix{{\bi x}_{k + 1} = {\bf \Phi} \left( {t_{k + 1}, t_k} \right){\bi x}_k + {\bi w}_k \hfill \cr {\bi y}_k {\bi = H}_k {\bi x}_k + {\bi R}_k } \right.$$

Φ(t _k+1, t _k) is the state transition matrix, and can be found by

(15)$$\left\{ \matrix{{\dot{\bf \Phi}} \left( {t,t_0} \right) = \displaystyle{{\partial \dot{\bi x}(t_0 )} \over {\partial {\bi x}(t_0 )}}{\bf \Phi} \left( {t,t_0} \right) \hfill \cr {\bf \Phi} \left( {t_0, t_0} \right){\bi \;= I}} \right.$$

The approximate method to calculate Φ(t _k+1, t_k) is

(16)$${\bf \Phi} \left( {t_{k + 1}, t_k} \right)= {\bi I} + {\bi F} \vert_{x_{{k + 1} \vert k}} (t_{k + 1} - t_k ) + \displaystyle{1 \over 2}{\bi F}^2 \vert_{x_{ {k + 1} \vert k}} (t_{k + 1} - t_k )^2 $$

where I is a 12 × 12 unit matrix in this system, and F is the Jacobi matrix of dynamics equation. F can be written as

(17)$$ {\bi F} \vert_{x_{ {k + 1} \Vert k}} = {\displaystyle{{\partial \dot{\bi x}} \over {\partial {\bi x}}}} \vert_{x_{ {k + 1} \vert k}} = {\displaystyle{{\partial ( \dot{\bi r}_1, \ddot{\bi r}_1, \dot{\bi r}_2, \ddot{\bi r}_2 )} \over {\partial ({\bi r}_1, \dot{\bi r}_1, {\bi r}_2, \dot{\bi r}_2 )}}} \vert_{x_{ {k + 1} \vert k}} = \left[ {\matrix{ {{\bf 0}_{3 \times 3}} & {{\bi I}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} \cr {\displaystyle{{\partial \ddot{\bi r}_1} \over {\partial {\bi r}_1}}} & {\displaystyle{{\partial \ddot{\bi r}_1} \over {\partial \dot{\bi r}_1}}} & {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} \cr {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {{\bi I}_{3 \times 3}} \cr {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {\displaystyle{{\partial \ddot{\bi r}_2} \over {\partial {\bi r}_2}}} & {\displaystyle{{\partial \ddot{\bi r}_2} \over {\partial \dot{\bi r}_2}}} \cr}} \right]$$

One can refer to the pioneering work of Kalman (Reference Kalman1960) and Kalman and Bucy (Reference Kalman and Bucy1961) for EKF.

A brief process for using the EKF is given here:

1. 1. The propagation of the estimated state from t _k to t _k+1, is subject to the dynamics model Equation (1), and ${\bar{\bi x}}_{k + 1} $ is the predicted state.

2. 2. ${\bar{\bi P}}_{k + 1}$ is the predicted covariance matrix, updated by

   (18)$${\bar{\bi P}}_{k + 1} = {\bf \Phi} \left( {t_{k + 1}, t_k} \right){\bi P}_{\bi k} {\bf \Phi} ^{\rm T} \left( {t_{k + 1}, t_k} \right) + {\bi w}_k $$

3. 3. Calculate the Kalman gain matrix through

   (19)$${\bi K}_{k + 1} = {\bar{\bi P}}_{k + 1} {\bi H}_{k + 1}^{\rm T} [{\bi H}_{k + 1}^{} {\bar{\bi P}}_{k + 1} {\bi H}_{k + 1}^{\rm T} + {\bi R}]^{ - 1} $$

4. 4. Update the estimated state through

   (20)$${\bi x}_{k + 1} = {\bar{\bi x}}_{k + 1} + {\bi K}_{k + 1} [{\bi y}_{k + 1} - {\bi H}_{k + 1} {\bar{\bi x}}_{k + 1} ]$$

5. 5. Update the covariance matrix by

   (21)$${\bi P}_{k + 1} = [{\bi I} - {\bi K}_{k + 1} {\bi H}_{k + 1} ]{\bar{\bi P}}_{k + 1} $$

The UKF is a new filtering algorithm for non-linear systems invented by Julier and Uhlmann (Reference Julier and Uhlmann1997). The state transition matrix and the Jacobi matrices of the dynamic equations and observation equations, necessary in EKF, are not required in UKF. The key point of the UKF is based on an Unscented Transformation (UT) that uses a set of appropriately chosen weighted points to parameterise the mean and covariance of the probability distributions. A detailed description of the UKF theory can be found in Julier and Uhlmann (Reference Julier and Uhlmann2004).

A brief process for using the UKF is given here:

1. 1. Create the set of sigma points by UT as follows

   (22)$$\left\{ \matrix{{\bi \chi} _0 = {\bi x}_k, w_0 = \displaystyle{\kappa \over {n + \kappa}} \hfill \cr {\bf \chi} _i = {\bi x}_k + \sqrt {n + \kappa} \left( {\sqrt {{\bi P}_{xx}}} \right)_i, w_i = \displaystyle{\kappa \over {2(n + \kappa )}} \hfill \cr {\bi \chi} _{i \;+\; n} = {\bi x}_k - \sqrt {n + \kappa} \left( {\sqrt {{\bi P}_{xx}}} \right)_i, w_{i + n} = \displaystyle{\kappa \over {2(n + \kappa )}} \hfill} \right.$$

   where, κ ∈ R, and usually n + κ = 3 when x is assumed to be drawn from Gaussian distribution. In this system, n = 12 and κ = − 9.

2. $\left( {\sqrt {{\bi P}_{xx}}} \right)_i $ is the i-th column of the matrix square root of P_xx

   w _i is the weighted factor associated with the i-th point.

3. 2. The propagation of the estimated state and the covariance matrix

   (23)$${\bar{\bi x}}_{k + 1\vert k} = \sum\limits_{i = 0}^{2n} {w_i} {\bi \chi} _{i,k + 1\vert k} $$
   (24)$${\bi P}_{xx,k + 1\vert k} = \sum\limits_{i = 0}^{2n} {w_i} \left[ {{\bi \chi} _{i,k + 1\vert k} - {\bar{\bi x}}_{k + 1\vert k}} \right]\left[ {{\bi \chi} _{i,k + 1\vert k} - {\bar{\bi x}}_{k + 1\vert k}} \right]^{\rm T} +\; {\bi w}$$

4. 3. The prediction of the observation and its covariance matrix are expressed by

   (25)$${\bi \gamma} _{i,k + 1\vert k} = {\bi h}({\bi \chi} _{i,k + 1\vert k} )$$
   (26)$$\tilde{\bi y}_{k + 1\vert k} = \sum\limits_{i = 0}^{2n} {w_i} {\bi \gamma} _{i,k + 1\vert k} $$
   (27)$${\bi P}_{yy,k + 1\vert k} = \sum\limits_{i = 0}^{2n} {w_i} \left[ {{\bi \gamma} _{i,k + 1\vert k} - \tilde{\bi y}_{k + 1\vert k}} \right]\left[ {{\bi \gamma} _{i,k + 1\vert k} - \tilde{\bi y}_{k + 1\vert k}} \right]^{\rm T} +\; {\bi R}_{k +1} $$

5. 4. Calculate the Kalman gain matrix as follows

   (28)$${\bi P}_{xy,k + 1\vert k} = \sum\limits_{i = 0}^{2n} {w_i} \left[ {{\bi \chi} _{i,k + 1\vert k} - {\bar{\bi x}}_{k + 1\vert k}} \right]\left[ {{\bi \gamma} _{i,k + 1\vert k} - \tilde{\bi y}_{k + 1\vert k}} \right]^{\rm T} $$
   (29)$${\bi K}_{k + 1} = {\bi P}_{xy,k + 1\vert k} \cdot {\bi P}_{yy,k + 1\vert k}^{ - 1} $$

   5. Update the estimated state and the covariance matrix as follows

   (30)$${\bi x}_{k + 1} = {\bar{\bi x}}_{k + 1\vert k} + {\bi K}_{k + 1} \left[ {{\bi y}_{k + 1} - \tilde{\bi y}_{k + 1\vert k}} \right]$$
   (31)$${\bi P}_{xx,k + 1} = {\bi P}_{xx,k + 1\vert k} - {\bi K}_{k + 1} \cdot {\bi P}_{yy,k + 1\vert k} \cdot {\bi K}_{k + 1}^{\rm T} $$

4.1. Scenario Setting

We set two scenarios as samples to substantiate our methods. The first scenario consists of a group of low-inclination orbits of LMO (a = 3996 km, e = 0·001, i = 0°) and HMO (a = 20430 km, e = 0·001, i = 0°). The second scenario consists of high-inclination orbit including LMO (a = 3696 km, e = 0·001, i = 75°) and HMO (a = 20000 km, e = 0·001, i = 75°). Table 3 lists the initial positions and velocities of the low-inclination scenario and Table 4 lists those of the high-inclination scenario in the J2000·0 Mars-equator coordinate System. The effect of probe attitude on visibility is not taken into account in this paper, because the probe can adjust attitude or sensor direction to be adapted to suit the observation.

Table 3.

Orbit parameters of low-inclination scenario.

Table 4.

Orbit parameters of high-inclination scenario.

According to the visibility condition in Section 2, Figures 7 and 8 give the visibility of Scenario-1 and Scenario-2, respectively. In these figures, labels 1, 2, 3, 4, and 5 on vertical axis represent SST, optical measurement from LMO to Phobos, optical measurement from LMO to Deimos, optical measurement from HMO to Phobos, and optical measurement from HMO to Deimos, respectively.

Figure 7.

The visibility of Scenario-1.

Figure 8.

The visibility of Scenario-2.

4.2. Simulation conditions

True data simulation is carried out by propagating the orbit from the initial state under a precise dynamic model (model-1). The true ephemeris is regarded as the criterion for comparison. The simplified dynamics model (model-2) is used for navigation. The difference between the precise and simplified models is taken as model error. The ephemeris of Mars and Martian moons are obtained through JPL's website http://ssd.jpl.nasa.gov/horizons.cgi.

The simulating range measurement error of SST is 10 m (1σ) and the range-rate error of SST is 0·1 m/s (1σ). These simulating measurement errors are in accord with current measurement level and the accuracy can be easily achieved for SST.

To make the measurement error approximate to true measurement condition, the optical measurement error is set according to the following the Equation (32).

(32)$$\sigma = \sigma _1 + \arctan \Big(\displaystyle{{\sigma _{Eph}} \over {\Delta R}}\Big) + \arctan \Big(\displaystyle{{R_{moon}} \over {\Delta R}}\Big)\sin \displaystyle{\psi \over 2}$$

where σ ₁ is the random error with Gaussian distribution of 0·01°, which includes measurement noise and the errors of the stellar location of background stars. σ _Eph is Martian moon ephemeris error of 10 km, ΔR is the range from the probe to Martian moon, R _moon is the approximate radius of Martian moon of 5 km and ψ is the phase angle from Martian moon to the Sun and the probe.

As the stellar location accuracy of the present stellar catalogue is far higher than 1″ (Hog et al., Reference Hog, Fabricius, Makarov, Urban, Corbin, Wycoff, Bastian, Schwekendiek and Wicenec2000), the effect of stellar location accuracy on navigation accuracy is included in the random error (σ ₁) of 0·01°.

In the filter process, the initial position and velocity uncertainty is set as 10 km and 5 m/s respectively for all simulation samples. The initial state covariance and the measurement noise covariance used for navigation filter were as follows:

$$\displaylines{{\bi P}_0 = diag[(10{\rm \;km)}^2, (10{\rm \;km)}^2, (10{\rm \;km)}^2, (5{\rm \;m/s)}^2, (5{\rm \;m/s)}^2, (5{\rm \;m/s)}^2, \cr (10{\rm \;km)}^2, (10{\rm \;km)}^2, (10{\rm \;km)}^2, (5{\rm \;m/s)}^2, (5{\rm \;m/s)}^2, (5{\rm \;m/s)}^2 ]} $$

${\bi R} = diag[(0.1^{\circ} {\rm )}^2, (0.1^{\circ} {\rm )}^2 ]$ for optical data and ${\bi R} = diag[(10{\rm \;m)}^2, (0.1{\rm \;m/s)}^2 ]$ for SST data.

4.3. Optical navigation by Viewing to Martian Moons Alone

Theoretically, by using optical measurement alone, LMO and HMO navigation can be performed respectively. In order to analyse the accuracy of the navigation using optical data alone, numerical simulations of autonomous navigation by using only optical data of Martian moons are carried out. The observation data is updated once per 1 minute. The predicted value will be given for navigation during the observation interval.

No matter what the inclination of the orbits, EKF for LMO cannot achieve results. The reason is that the non-linearity of LMO is stronger than that of HMO in the same time span. The results of EKF for LMO by optical measurement alone are divergent due to the inability to acceptably linearize a strongly non-linear system. The reason that the results of UKF for LMO are convergent is due to its properties without linearization. However, the navigation accuracy of UKF is still low, 10 km position error and 8 m/s velocity error in Scenario-1, and 5 km position error and 6 m/s velocity error in Scenario-2.

The simulated navigation results of HMO in Scenario-1 are shown in Figure 9 and those in Scenario-2 are shown in Figure 10. The results indicate that autonomous navigation by using optical data alone for HMO, EKF and UKF lead to similar results once converged, and the accuracy of UKF is a little better than that of EKF. The maximum position error is about 10 km and the maximum velocity error is about 0·7 m/s for HMO in Scenario-1, and that in Scenario-2 is 5 km and 0·3 m/s respectively.

Figure 9.

Estimate errors of HMO of low-inclination by optical data alone.

Figure 10.

Estimate errors of HMO of high inclination by optical data alone.

4.4. Navigation by combining Optical Data and SST Data

By combining SST data and optical measurement, autonomous navigation of the two probes must be conducted simultaneously. In terms of the given simulation conditions mentioned above, a series of numerical simulations of autonomous navigation by combining optical data of viewing Martian moons and Satellite-to-Satellite radio tracking data are carried out in this section. Similarly, the predicted value will be given for navigation during the observation interval. In order to analyse the navigation performance with different measurement update rate, the measurement date update rate is set to three cases: 10 minutes, 1 minute, and 10 seconds, respectively.

Tables 5 and 6 give the statistical information on the navigation accuracy of Scenario-1 after the filter has converged from the 24th to 48th hour. The simulation navigation results in Scenario-1 with data update rate of 1 minute are shown in Figures 11 and 12.

Figure 11.

Estimate errors of LMO in Scenario-1 with 1 minute update rate.

Table 5.

Accuracy comparison of LMO with low inclination.

Table 6.

Accuracy comparison of HMO with low inclination.

In Scenario-1, when the update rate of measurement data is set to be 10 minutes, no matter LMO or HMO, the navigation accuracies by UKF are both better than that of EKF. For LMO, the position and velocity errors of EKF are about 2·2 km of root mean square (RMS) and 1·9 m/s (RMS), respectively. By contrast, those of UKF are 1·8 km and 1·5 m/s, respectively. For HMO, the position and velocity errors of EKF are about 3 km (RMS) and 0·2 m/s (RMS), respectively. By contrast, those of UKF are 1 km and 0·08 m/s, respectively. Especially for HMO, the performance in accuracy of UKF is about three times better than that of EKF.

However, when the update rate of measurement data is set to be either ten seconds or one minute, no matter LMO or HMO, EKF and UKF have similar excellent performance. For LMO, the position errors of both EKF and UKF are about 1·6 km (RMS) and 2·8 km (Maximum), and the velocity errors are about 1·3 m/s (RMS) and 2·2 m/s (Maximum). Meanwhile, for HMO, in the same scenario, the performance in accuracy of UKF is a little better than that of EKF. The position and velocity errors of EKF are about 1·2 km (RMS) and 0·07 m/s (RMS), respectively. By contrast, those of UKF are 0·8 km and 0·05 m/s, respectively.

Tables 7 and 8 give statistical information of navigation results of LMO and HMO in Scenario-2 after the filter has converged from the 24th to 48th hour. The simulated navigation results for the cases of LMO and HMO in Scenario-2 of high inclination orbits are shown in Figures 13 and 14.

Figure 12.

Estimate errors of HMO in Scenario-1 with 1 minute update rate.

Figure 13.

Estimate errors of LMO in Scenario-2 with 1 minute update rate.

Figure 14.

Estimate errors of HMO in Scenario-2 with 1 minute update rate.

Table 7.

Accuracy comparison of LMO with high inclination.

Table 8.

Accuracy comparison of HMO with high inclination.

In Scenario-2, with the update rate of measurement data set to be ten minutes, no matter LMO or HMO, the UKF navigation accuracy is also better than that of EKF. For LMO, the position and velocity errors of EKF are about 1·1 km Root Mean Square (RMS) and 1·0 m/s (RMS), respectively. By contrast, those of UKF are 0·65 km and 0·6 m/s, respectively. For HMO, the position and velocity errors of EKF are about 3·5 km (RMS) and 0·25 m/s (RMS), respectively. By contrast, those of UKF are 2·5 km and 0·18 m/s, respectively.

Similarly to Scenario-1, when the update rate of measurement data is set to be either ten seconds or one minute, EKF and UKF have similar excellent performance for LMO. For HMO, the navigation accuracy with data update rate of ten seconds is little better than that of one minute. For LMO, the position errors of both EKF and UKF are about 0·8–0·9 km (RMS) and 1·7–2·0 km (Maximum), and the velocity errors are about 0·7–0·8 m/s (RMS) and 1·7–2·1 m/s (Maximum). Meanwhile, for HMO, in the same scenario, the performance in accuracy of UKF is a little better than that of EKF. For data update rate of one minute, the position and velocity errors of EKF are about 2·2 km (RMS) and 0·15 m/s (RMS), respectively. By contrast, those of UKF are about 1·6 km and 0·1 m/s, respectively.

4.5. Simulation conclusion

From the above simulation results, we conclude that

• • No matter whether in Scenario-1 with low inclination orbits or Scenario-2 with high inclination orbits, the autonomous navigation method of Mars probes by combining optical data of viewing Martian moons and Satellite-to-Satellite radio tracking data is practicable.

• • EKF and UKF have similar excellent performance with ten seconds and one minute data update rates. Generally, the performance in accuracy of UKF is a little better than that that of EKF. If the measurement data update rate is ten minutes, UKF navigation accuracy is obviously better than EKF. Due to the dynamics model error, the larger the measurement interval is, the more obviously the strong non-linearity occurs. The EKF algorithm requires linearization to the dynamics model but the UKF algorithm does not.

• • Because there are 12 state variables in the system, 25 sigma points should be created by UT, and the time consumption of UKF is twenty times more than that of EKF. So EKF performs faster.

• • The shorter the measurement data update rate is, the higher the navigation accuracy is. If the measurement data update rate is shorter than one minute, the navigation accuracy is little improved.

• • With data update rate less than one minute, for LMO, the autonomous navigation accuracy can attain about 1·6 km (RMS), 2·0–3·0 km (Maximum) in position, and about 1·0 m/s (RMS), 2·0–3·0 m/s (Maximum) in velocity. And for HMO, the autonomous navigation accuracy can attain about 1·5 km (RMS), 3–4 km (Maximum) in position, and about 0·1 m/s (RMS), 0·2 m/s (Maximum) in velocity.