2. MARS ENTRY DYNAMIC EQUATIONS

System dynamic equations are necessary to predict the state variables of a Mars entry vehicle, therefore the applicability of an integrated navigation filter heavily depends on the availability of sufficiently accurate dynamic equations of Mars entry. The traditional 3-DOF dynamic model only represents the Mars atmospheric entry translational dynamics, excluding the attitude dynamics and kinematics (Vinh et al., Reference Vinh, Busemann and Culp1980; Lockwood et al., Reference Lockwood, Powell, Graves and Carman2001). As only aerodynamic lift is used to reduce the downrange footprint dispersion during Mars entry, the bank angle of entry vehicles is considered as the only control variable. Therefore, the attitude angle should be included into Mars entry dynamic equations to completely and accurately describe the state variables of entry vehicles. However, the aerodynamic torque is very difficult to precisely model due to the uncertainties in vehicle aerodynamics and inertia matrices. In this paper, only the attitude kinematics equation is embedded into the system dynamics model, the angular velocity information is directly obtained from a gyro free of modelling errors. Therefore, the attitude dynamics are no longer needed. Then, the 6-DOF Mars entry dynamic equations can be represented in the Mars centred inertial coordinate system as follows:

(1)$$\eqalign{{\bf \dot r} & = {\bf v} \cr {\bf \dot v} & = {\bf T}_V^I {\bf a}_V + {\bf T}_G^I {\bf g}_G \cr {\bf \dot e} & = {\bf K\tilde \omega} _B} $$

where r=[r _x,r_y,r_z]^T is the position vector from the centre of Mars to the vehicle's centre of mass, v=[v _x,v_y,v_z]^T is the velocity vector of the entry vehicle, e=[φ,ϑ,φ]^T is the tri-axial attitude angle. T_V^I is the coordinate transformation matrix from the velocity coordinate system to the Mars centred inertial coordinate system, T_G^I is the coordinate transformation matrix from the geographic coordinate system to the Mars centred inertial coordinate system. The definition of coordinate transformation matrices T_V^I and T_V^I can be found in Li and Peng (Reference Li and Peng2011), and is thus not repeated here. a_V is the aerodynamic acceleration and can be easily constructed according to Equation (5). g_G is the Mars gravity acceleration and is defined in Equation (2), where μ_M=GM _Mars is the Mars gravitational constant. Coefficient matrix K is the function of attitude angle φ,ϑ and defined in Equation (3). ${\bf \tilde \omega} _B $ is the angular velocity outputs from a gyro defined in Equation (4), ω_B=[ω ₁,ω ₂,ω ₃]^T is the real triaxial angle velocity described in the body-fixed coordinate system, b_ω is the angular rate bias, and ξ_ω is the white Gaussian angular rate output noise.

(2)$${\bf g}_M = \left[ {0,0, - \displaystyle{{\mu _M} \over {\left| {\bf r} \right|^2}}} \right]^T $$

(3)$${\bf K} = \displaystyle{1 \over {\cos \vartheta}} \left[ {\matrix{ {\cos \vartheta} & {\sin \vartheta \sin \varphi} & {\sin \vartheta \cos \varphi} \cr 0 & {\cos \vartheta \cos \varphi} & { - \cos \vartheta \sin \varphi} \cr 0 & {\sin \varphi} & {\cos \varphi} \cr}} \right]$$

(4)$${\bf \tilde \omega} _B = {\bf \omega} _B + {\bf b}_\omega + {\bf \xi} _\omega $$

(5)$${\bf a}_V = [\matrix{ { - D} & { - L\sin \sigma} & {L\cos \sigma} \cr} ]^{\rm T} $$

where L and D are the aerodynamic lift and drag accelerations, defined by

(6)$$L = \displaystyle{1 \over 2}\rho v^2 \displaystyle{{C_L S} \over m}$$

(7)$$D = \displaystyle{1 \over 2}\rho v^2 \displaystyle{{C_D S} \over m}$$

where ρ is the reference Mars atmospheric density defined in Equation (8), C _L and C _D are the aerodynamic lift and drag coefficients respectively, S represents the vehicle reference surface area, and m is the mass of the entry vehicle.

The high-precision engineering-level Mars atmosphere model Mars-GRAM is adopted in this paper (Martin et al., Reference Martin, Wong and Lee2013; Tolson and Prince, Reference Tolson and Prince2011). According to the Mars-GRAM model, the reference Mars atmospheric density is defined as follows:

(8)$$\rho = \rho _0 \exp ( - \beta (h))$$

where ρ ₀=559·351005946503/c ₁c ₂, c ₁=188·95110711075, c ₂=1·4×10⁻¹³h ³−8·85×10⁻⁹h ²−1·245×10⁻³h+205·3645, β(h)=−0·000105 h.

As the true Mars atmospheric density varies with the different seasons and atmospheric temperature, the value obtained from Equation (8) is only an approximation of the true density. When we perform high-fidelity navigation analysis of Mars entry, the deviation in the atmospheric density must be taken into account. The true Mars atmospheric density can be formulated as follows:

(9)$$\rho ^ * = \rho \cdot (1 + \Delta )$$

where Δ denotes the percentage of deviation in Mars atmospheric density.

It should be noted that the aerodynamics lift acceleration L and drag acceleration D in Equations (5), (6) and (7) are closely related to the reference atmospheric density ρ, therefore, the dynamics equations in Equation (1) will inevitably introduce model and parameter uncertainties in the case that there are large uncertainties on the reference Martian atmospheric density.

3. NAVIGATION MEASUREMENT MODELS

3.1. Accelerometer measurement models

An accelerometer is designed to measure the linear acceleration along three orthogonal axes. The acceleration measured by an accelerometer is represented as follows:

(10)$${\bf \tilde a}_B = {\bf a}_B + {\bf b}_a + {\bf \xi} _a $$

where ${\bf \tilde a}_B $ is the linear accelerometer output along body axes, a_B is the true linear acceleration, b_a is the acceleration bias, ξ_a is the white Gaussian acceleration output noise.

Then, accelerometer measurement model is defined as follows:

(11)$${\bf y}_1 = {\bf \tilde a}_B = T_V^B {\bf a}_V^ * + {\bf b}_a + {\bf \xi} _a $$

where

(12)$${\bf a}_V^ * = [\matrix{ { - D^*} & { - L^* \sin \sigma} & {L^* \cos \sigma} \cr} ]^{\rm T} $$

(13)$$L^* = \displaystyle{1 \over 2}\rho ^* v^2 \displaystyle{{C_L S} \over m}$$

(14)$$D^* = \displaystyle{1 \over 2}\rho ^* v^2 \displaystyle{{C_D S} \over m}$$

Since the accelerometer is utilized to measure the real aerodynamic acceleration exerted on the entry vehicle, it must be stressed that the Mars atmospheric density involved in Equations (13) and (14) should be true atmospheric density ρ*, not reference atmospheric density ρ. At the same time, the basis and noise terms in Equation (11) are considered as measurement noises and included into the subsequent navigation filter.

3.2. Radio measurement models

The navigation system for Mars entry incorporates an assumed conglomeration of orbiters previously launched into Mars orbit. The orbiter usually includes two on board navigation sensors: two-way range radio and two-way Doppler radio. The two-way range radio will take range measurements from the entry vehicle to each of the orbiting radio beacons, provided that they are within line-of-sight. The two-way Doppler will take velocity measurements with the orbiting beacons (Boehmer, Reference Boehmer1998).

As the whole process of Mars entry only lasts a relatively short time, the influence of perturbation on the orbit of the orbiter can be neglected here. Thus, the position and velocity of the orbiter can be obtained in the Mars centred inertial coordinate system according to the simple two-body gravity theory

(15)$$\displaystyle{{d{\bf r}_o} \over {dt}} = {\bf v}_o $$

(16)$$\displaystyle{{d{\bf v}_o} \over {dt}} = - \displaystyle{{\mu _M} \over {\left| {{\bf r}_o} \right|^3}} {\bf r}_o $$

The distance between the entry vehicle and an orbiter can be reconstructed as follows:

(17)$$\tilde R_i = R_i + \zeta _{Ri} $$

(18)$$R_i = \sqrt {({\bf r} - {\bf r}_i )^T ({\bf r} - {\bf r}_i )} $$

where r is the position vector of entry vehicle, r_i is the position vector of the i^th orbiting radio beacon. ζ _R is the range measurement noise. Both vectors mentioned above are defined in the Mars centred inertial coordinate system.

The Doppler measurement, or the rate of change of range measurement, from the entry vehicle to an orbiter is designated as

(19)$$\tilde V_i = V_i + \zeta _{Vi} $$

(20)$$V_i = \displaystyle{{({\bf r} - {\bf r}_i )^T ({\bf v} - {\bf v}_i )} \over {R_i}} $$

where v is the velocity vector of entry vehicle, v_i is the velocity vector of the i^th orbiting radio beacon. ζ _V is the velocity measurement noise. Both vectors mentioned above are also defined in the Mars centred inertial coordinate system.

Then, the radio measurement model can be constructed as follows:

(21)$${\bf y}_2 = \left[ \matrix{{{\bf \tilde R}} \cr {{\bf \tilde V}} } \right] = \left[ \matrix{{\bf R} \cr {\bf V} } \right] + \left[ {\matrix{ {{\bf \zeta} _R} \cr {{\bf \zeta} _V} \cr}} \right]$$

where R=[R ₁, …, R _m]^T, V=[V ₁, …, V _m]^T, ζ_R=[ζ_R1, …, ζ_Rm]^T, ζ_V=[ζ_V1, …, ζ_Vm]^T, and subscript m denotes the number of radio orbiting used in the navigation scheme.

4. DESENSITISED EXTENDED KALMAN FILTER DESIGN

To estimate the state variables of entry vehicles and suppress navigation measurement noises, an integrated navigation filter is designed by use of a desensitised extended Kalman filter (DEKF). Moreover, DEKF is adopted to efficiently reduce the sensitivity of state variables with respect to model and parameter uncertainties.

4.1. System equations

If we select ${\bf X}{\rm (}t{\rm )} = [{\bf r}^T, {\bf v}^T, {\bf e}^T ]_{9 \times 1}^T $ as the state variables, then the Mars entry dynamic equations Equation (1) can be rewritten as follows:

(22)$${\bf \dot X}(t) = \left[ {\matrix{ {\bf v} \cr {{\bf T}_V^I {\bf a}_V^* + {\bf T}_G^I {\bf g}_G} \cr {{\bf K\tilde \omega} _B} \cr}} \right]_{9 \times 1} = {\bf f}({\bf X}(t),t) + {\bf w}$$

Similarly, navigation observation equations defined in Equations (11) and (21) can also be written as follows:

(23)$${\bf Y}(t) = \left[ \matrix{{\bf y}_1 \cr {\bf y}_2 } \right] = \left[ \matrix{{\bf \tilde a}_B \cr {{\bf \tilde R}} \cr {{\bf \tilde V}} } \right]_{2m \times 1} = {\bf h}({\bf X}(t),t) + {\bf \upsilon} $$

where m stands for the number of radio beacons used in the integrated navigation.

The system equations utilized in the subsequent navigation filter can thus be defined as follows:

(24)$${\bf \dot X}_k = {\bf f}{\rm (}{\bf X}_{k - 1} {\rm )} + {\bf \Gamma w}_{k - 1} $$

(25)$${\bf y}{\rm (}t_k {\rm )} = {\bf h}{\rm (}{\bf X}_k {\rm )} + {\bf \upsilon} _k $$

where Γ is the model noise input matrix. w_k and υ_k represent the system process noise and measurement noise respectively, they are assumed to be independent of each other, white and with normal probability distributions

(26)$$p({\bf w})\sim N(0,{\bf Q})$$

(27)$$p({\bf \upsilon} )\sim N(0,{\bf R})$$

(28)$${\rm Cov(}{\bf w}_k, {\bf \upsilon} _j {\rm ) = E(}{\bf w}_k {\bf \upsilon} _j^{\rm T} {\rm )} = 0$$

Defining state transfer matrix Φ_{k+1/ k}

(29)$${\bf \Phi} _{k/k - 1} = {\bf I} + {\bf F}_{k/k - 1} \cdot \Delta t$$

(30)$${\bf \dot \Phi} (t,\tau ) = {\bf F}(t){\bf \Phi} (t,\tau )$$

(31)$${\bf \Phi} (t,t) = {\bf I}$$

where F is the Jacobian matrix of partial derivatives of f with respect to state variable X, that is

(32)$${\bf F}_{k/k - 1} = \left. {\displaystyle{{\partial {\bf f}({\bf X},k)} \over {\partial {\bf X}}}} \right|_{{\bf X} = {\bf X}_{k - 1}}=\left[ { .8pt\matrix{ {{\bf 0}_{3 \times 3}} & {{\bf I}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} \cr {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {\displaystyle{{\partial ({\bf T}_V^I {\bf a}_V^ * + {\bf T}_G^I {\bf g}_G )} \over {\partial {\bf e}}}} \cr {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} & {\displaystyle{{\partial ({\bf K\tilde \omega} _B )} \over {\partial {\bf e}}}} \cr}} \right]_{9 \times 9} $$

Defining sensitivity matrix H_k

(33)$${\bf H}_k = \left. {\displaystyle{{\partial {\bf h}{\rm (}{\bf X}{\rm )}} \over {\partial {\bf X}}}} \right|_{{\bf X} = {\bf \hat X}_{k/k - 1}} = \left[ {.8pt\matrix{ {{\bf 0}_{3 \times 3}} & {\displaystyle{{\partial {\bf \tilde a}_B} \over {\partial {\bf v}^{\rm T}}}} & {{\bf 0}_{3 \times 3}} \cr {\displaystyle{{\partial {\bf \tilde R}} \over {\partial {\bf r}^{\rm T}}}} & {{\bf 0}_{3 \times 3}} & {{\bf 0}_{3 \times 3}} \cr {\displaystyle{{\partial {\bf \tilde V}} \over {\partial {\bf r}^{\rm T}}}} & {\displaystyle{{\partial {\bf \tilde V}} \over {\partial {\bf v}^{\rm T}}}} & {{\bf 0}_{3 \times 3}} \cr}} \right]_{9 \times 9} $$

4.2. Desensitised extended Kalman filter

The standard extended Kalman filter can be considered as a specified form of desensitised extended Kalman filter in a certain sense, so the formulation of DEKF is almost exactly the same as for a standard EKF except for the definition of the gain matrix.

Initialization: for k = 0, set

(34)$${\bf \hat X}_0 = E[{\bf X}_0 ]$$

(35)$${\bf P}_0 = E\left[ {({\bf X}_0 - E[{\bf X}_0 ])({\bf X}_0 - E[{\bf X}_0 ])^{\rm T}} \right]$$

DEKF time update equations:

(36)$${\bf \hat X\dot}_{k/k - 1} = {\bf f}(k,{\bf \hat X}_{k - 1} )$$

(37)$${\bf P}_{k/k - 1} = {\bf \Phi} _{k/k - 1} {\bf P}_{k - 1} {\bf \Phi} _{k/k - 1}^{\rm T} + {\bf Q}_{k - 1} $$

DEKF measurement update equations:

(38)$${\bf \hat X}_{k/k} = {\bf \hat X}_{k/k - 1} + {\bf K}_k \left[ {{\bf y}(t_k ) - {\bf h}({\bf \hat X}_{k/k - 1} )} \right]$$

(39)$${\bf P}_k = ({\bf I} - {\bf K}_k {\bf H}_k ){\bf P}_{k/k - 1} ({\bf I} - {\bf K}_k {\bf H}_k )^{\rm T} + {\bf K}_k {\bf R}_k {\bf K}_k^{\rm T} $$

where ${\bf \hat X}_{k/k - 1} $ is one step predicted value of the state at k moment, ${\bf \hat X}_{k/k} $ is the state estimation at k moment, y_k is the observed quantity, K_k is the filtering gain, P_{k / k} is the error variance matrix, Φ_{k/k− 1} is the state transition matrix. Q_k is the process noise variance matrix, defined in Equation (40), and R_k is measurement noise variance matrix.

(40)$${\bf Q}_k = {\bf Q}(t)T$$

The key of the DEKF is to derive the desensitised optimal gain matrix K_k according to the DOC theory. Firstly, the predicted sensitivity of state variables with respect to uncertainties ${\bf \hat \sigma} _{k/k - 1} = \displaystyle{{\partial {\bf \hat X}_{k/k - 1}} \over {\partial \rho}} $ in differential form can be derived by taking partial derivatives of the state prediction equation in Equation (36) (Karlgaard and Shen, Reference Karlgaard and Shen2011):

(41)$$\eqalign{{\bf \hat \sigma\dot}_{k/k - 1} & = \displaystyle{{\partial {\bf \hat X\dot}_{k/k - 1}} \over {\partial \rho}} = \displaystyle{{\partial {\bf f}(k,{\bf \hat X}_{k - 1} )} \over {\partial \rho}} \cr & = {\bf F}_{k/k - 1} \displaystyle{{\partial {\bf \hat X}_{k/k - 1}} \over {\partial \rho}} + \displaystyle{{\partial {\bf F}_{k/k - 1}} \over {\partial \rho}} {\bf \hat X}_{k/k - 1} \cr & = {\bf F}_{k/k - 1} {\bf \hat \sigma} _{k/k - 1} + \displaystyle{{\partial {\bf f}({\bf X}(t_k ),t_k )} \over {\partial \rho}}} $$

Then, the predicted sensitivity is computed by numerical integration of the equation. In the same way, the corrected sensitivity at the measurement update can be obtained by taking partial derivatives of the corrections equation in Equation (38) (Karlgaard and Shen, Reference Karlgaard and Shen2011):

(42)$$\eqalign{{\bf \hat \sigma} _{k/k} & = \displaystyle{{\partial {\bf \hat X}_{k/k}} \over {\partial \rho}} = \displaystyle{{\partial {\bf \hat X}_{k/k - 1}} \over {\partial \rho}} + \displaystyle{{\partial {\bf K}_k} \over {\partial \rho}} \left[ {{\bf y}(t_k ) - {\bf h}({\bf \hat X}_{k/k - 1} )} \right] + {\bf K}_k \displaystyle{{\partial \left[ {{\bf y}(t_k ) - {\bf h}({\bf \hat X}_{k/k - 1} )} \right]} \over {\partial \rho}} \cr & = {\bf \hat \sigma} _{k/k - 1} + {\bf K}_k \displaystyle{{\partial {\bf y}(t_k )} \over {\partial \rho}} - {\bf K}_k \displaystyle{{\partial {\bf h}({\bf \hat X}_{k/k - 1} )} \over {\partial \rho}} \cr & = {\bf \hat \sigma} _{k/k - 1} - {\bf K}_k \displaystyle{{\partial {\bf h}({\bf \hat X}_{k/k - 1} )} \over {\partial \rho}} \cr & = {\bf \hat \sigma} _{k/k - 1} - {\bf K}_k {\bf B}_k} $$

where ${\bf B}_k = {\bf H}_k {\bf \hat \sigma} _{k/k - 1} + {{\partial {\bf h}({\bf X}(t_k ),t_k )} / {\partial \rho}} $. It should be pointed out that the true measurement sensitivity is (∂ y_k/∂ρ)=0 in this formulation. At the same time, the value of (∂ K_k/∂ρ) is assumed to be zero here. As any non-zero value of (∂ K_k/∂ρ) implies that the solution for the optimal gain is a function of the residual $[{\bf y}(t_k ) - {\bf h}({\bf \hat X}_{k/k - 1} )]$, which obviously violates the previous assumption of the linear update equation given in Equation (39) (Karlgaard and Shen, Reference Karlgaard and Shen2011; Shen and Karlgaard, Reference Shen and Karlgaard2012). This assumption implies that the desensitisation method only penalizes an approximate sensitivity rather than the true sensitivity.

To determine the desensitised optimal gain K_k, the traditional cost function of the Kalman filter is augmented with a penalty function consisting of a weighted norm of the posterior sensitivity, given by

(43)$$J = {\rm trace}({\bf \hat P}_k ) + {\bf \hat \sigma} _{k/k}^{\rm T} {\bf \chi \hat \sigma} _{k/k} $$

where χ is a symmetric positive semi-definite weighting matrix for the sensitivity. Substituting Equations (39) and (42) into Equation (43) and taking the derivative with respect to K_k term by term yields

(44)$$\eqalign{\displaystyle{{\partial J} \over {\partial {\bf K}_k}} & = \displaystyle{{\partial {\rm trace}({\bf \hat P}_k )} \over {\partial {\bf K}_k}} + \displaystyle{{\partial ({\bf \hat \sigma} _{k/k}^{\rm T} {\bf \chi \hat \sigma} _{k/k} )} \over {\partial {\bf K}_k}} \cr & = 2{\bf S}_k {\bf K}_k^{\rm T} - 2{\bf H}_k {\bf P}_{k/k - 1} + 2({\bf B}_k {\bf B}_k^{\rm T} {\bf K}_k^{\rm T} {\bf \chi} - {\bf B}_k {\bf \hat \sigma} _{k/k - 1}^{\rm T} {\bf \chi} )} $$

In order to minimize J, setting (∂J/∂ K_k)=0, then the desensitised optimal gain matrix can be solved as follows

(45)$${\bf K}_k = ({\bf P}_{k/k - 1} {\bf H}_k^{\rm T} + {\bf \chi \hat \sigma} _{k/k - 1} {\bf B}_k^{\rm T} ) \cdot ({\bf S}_k + {\bf \chi B}_k {\bf B}_k^{\rm T} )^{ - 1} $$

where S_k=H_kP_{k/k− 1}H_k^T+R_k.

As expected, the standard EKF gain matrix can be recovered by setting χ =0.

5. NUMERICAL SIMULATION AND ANALYSIS

To confirm the validity of the Mars entry integrated navigation algorithm developed in this paper, numerical simulation in a MATLAB/Simulink environment has been carried out. Conventional 3-DOF Mars entry dynamic equations and simple rigid-body attitude dynamics are combined to produce the nominal state variables (Vinh et al., Reference Vinh, Busemann and Culp1980; Ely et al., Reference Ely, Bishop and Dubois-Matra2001). The entry vehicle model is selected to be similar to the MSL lander and the relative physical parameters are given in Table 1. The initial state parameters of the Mars entry vehicle is provided by the deep space network (DSN), and the initial state variables are described in Table 2 according to the MSL simulation data in Martin et al. (Reference Martin, Wong and Lee2013), Tolson et al. (Reference Tolson and Prince2011), Chen et al. (Reference Chen, Vasavada, Cianciolo, Barnes, Tyler, Rafkin, Hinson and Lewis2010) and Schoenenberger et al. (Reference Schoenenberger, Norman, Dyakonov, Karlgaard, Way and Kutty2013). We assumed that three orbiting beacons are available during the Mars atmospheric entry phase, the initial state variables and the errors of orbiting beacons are given in Table 3. The accelerometer bias b_a and gyro bias b_ω applying for the simulation are assumed to be [3×10⁻³, 3×10⁻³, 3×10⁻³] m/s ² and [5×10⁻⁶, 5×10⁻⁶, 5×10⁻⁶] rad/s respectively, the covariance matrix of white Gaussian noise ξ_a and ξ_ω are set to [5×10⁻⁸, 5×10⁻⁸, 5×10⁻⁸] m ^2/s ⁴ and [9×10⁻⁸, 9×10⁻⁸, 9×10⁻⁸] (rad/s)² respectively. Radio range and velocity measurement noise variance in navigation measurement model are assumed to be 100 m² and 0·1 m²/s ² respectively. Initial state error covariance matrix P₀, system process noise variance matrix Q and penalty matrix P₀, are respectively assumed to be as follows: ${\bf P}_0 = \left[ {\matrix{ {10^{13} {\bf I}_{3 \times 3}} & {} & {} \cr {} & {10^8 {\bf I}_{3 \times 3}} & {} \cr {} & {} & {{\bf I}_{3 \times 3}} \cr}} \right]_{9 \times 9} $, $\quad{\bf Q} = \left[ {\matrix{ {10^6 {\bf I}_{3 \times 3}} & {} & {} \cr {} & {10^4 {\bf I}_{3 \times 3}} & {} \cr {} & {} & {10^{ - 5} {\bf I}_{3 \times 3}} \cr}} \right]_{9 \times 9} $, ${\bf R}_k = \left[ {\matrix{ {10^{ - 4} {\bf I}_{3 \times 3}} & {} & {} \cr {} & {10^2 {\bf I}_{3 \times 3}} & {} \cr {} & {} & {10^{ - 1} {\bf I}_{3 \times 3}} \cr}} \right]_{9 \times 9} $, χ=10⁻⁷I_9×9. As the whole Mars entry process only lasts a short time (about 4 minutes in the MSL mission) (Martin et al., Reference Martin, Wong and Lee2013), the influence of perturbation on the orbit is neglected in our simulation, the position and velocity of each orbiter are obtained by a simple two-body model. The planned entry time span is assumed to be 250 seconds and the simulation sample step is set to 0·5 second. A four-order Runge-Kutta algorithm is selected as the numerical solver of the integral entry dynamic equations. Other relative parameters used in the simulation are set as follows: Mars gravitational constant μ _M=4·282829×10⁴ km ³/km ³s ², Mars rotating angular rate ω _M=7·0882×10⁻⁵rad/s.

Table 1.

Physical parameters of entry vehicle.

Table 2.

Initial state variables of entry vehicle.

Table 3.

Initial state and errors of orbiting beacons.

It is believed that there is a maximum of ±15% deviation in the nominal values from the Mars-GRAM model defined in Equation (8) when compared with the real Mars atmosphere density, and the deviation roughly follows a normal distribution around the nominal values (Chen et al., Reference Chen, Vasavada, Cianciolo, Barnes, Tyler, Rafkin, Hinson and Lewis2010; Schoenenberger et al., Reference Schoenenberger, Norman, Dyakonov, Karlgaard, Way and Kutty2013; Chen et al., Reference Chen, Beck, Brugarolas, Edquist, Mendeck, Schoenenberger and Way2013). In order to simulate and analyse the effect of the uncertain disturbance in Mars atmospheric density, a random deviation, obeying the normal distribution with the standard deviation of 15%, is intentionally included in our simulation.

The performance of both DEKF-based integrated navigation and standard EKF-based integrated navigation were tested for comparison in our simulations. The differences between the referenced state variables of the entry vehicle and estimated state variables, that is navigation errors, are plotted in Figures 1–6. Triaxial position estimation errors and velocity estimation errors from standard EKF-based integrated navigation are shown in Figures 1 and 2 respectively, and triaxial position estimation errors and velocity estimation errors from DEKF-based integrated navigation are also plotted in Figures 3 and 4 respectively. It can be seen from Figures 1 to 4 that both the position estimation errors and velocity estimation errors from DEKF-based integrated navigation can clearly be reduced to a smaller magnitude when compared with those of EKF-based integrated navigation, which indicates the performance of DEKF is superior to that of EKF in the presence of larger model and parameter uncertainties. It should be noted from Figures 1 to 4 that x-axial navigation errors are significantly smaller than those from the other two axes, which could be interpreted as meaning that the state variables along the x-axial direction have a better observability in our navigation geometric configuration.

Figure 1.

Triaxial position estimation errors: standard EKF based integrated navigation.

Figure 2.

Triaxial velocity estimation errors: standard EKF-based integrated navigation.

Figure 3.

Triaxial position estimation errors: DEKF-based integrated navigation.

Figure 4.

Triaxial velocity estimation errors: DEKF-based integrated navigation.

Figure 5.

Total position estimation errors: DEKF and standard EKF-based integrated navigation.

Figure 6.

Total velocity estimation errors: DEKF and standard EKF-based integrated navigation.

In order to further analyse the navigation performance under larger model and parameter uncertainties, total position estimation errors and velocity estimation errors from both DEKF and standard EKF-based integrated navigations are plotted in Figures 5 and 6 for comparison. The sub-figures in Figures 5 and 6 depict the partially amplified figures of the total position estimation error and total velocity estimation error form DEKF-based integrated navigation between 200 to 250 seconds respectively. It can be seen form the Figures 5 and 6 that the estimation errors of DEKF-based integrated navigation are fairly small with position errors less than 100 m and velocity errors less than 6·5 m/s, which can meet the navigation need of future pinpoint Mars landing missions. However, if a standard EKF filtering is adopted instead of DEKF, the performance of the integrated navigation algorithm degrades and can only achieve 800 m position error and 60 m/s velocity error with the same simulation condition. It is interesting to note that the performance of DEKF-based integrated navigation is superior to the results reported in Chen (Reference Chen, Beck, Brugarolas, Edquist, Mendeck, Schoenenberger and Way2013), and slightly inferior to the results depicted in Li and Peng (Reference Li and Peng2011) even if a larger uncertain disturbance in Mars atmospheric density is included in our simulation. Based on the simulation results and analysis mentioned above, it can be safely concluded that DEKF-based integrated navigation can efficiently reduce the sensitivity of state variable with respect to uncertainties in the dynamics system and significantly improve the accuracy of state estimation in the presence of a larger uncertainty disturbance.

As the DEKF adopts an almost identical iterative formula as standard EKF, it is expected that the corresponding computing time and burden should not be significantly increased compared with EKF. In our simulation, the DEKF-based integrated navigation algorithm takes 43 milliseconds, while the standard EKF-based integrated navigation algorithm takes 40 milliseconds. Both simulations were run on a laptop PC of Intel Core (TM) i7-3610QM CPU @ 2.3 GHz. It is thus clear that the computing time of the DEKF algorithm is slightly increased when compared to that of standard EKF algorithm.