1. INTRODUCTION
One of the significant engineering challenges for precisely landing a vehicle on the surface of Mars involves the entry, descent and landing (EDL) phase of the mission. Traditional Mars entry vehicles, such as the Viking and Mars exploration rovers, adopted the inertial measurement unit (IMU)-based dead reckoning navigation mode and the unguided ballistic trajectory entry without aerodynamic lift control, which leads to a larger landing error ellipse in the order of several hundred kilometres (Braun, Reference Braun2007; Lu et al., Reference Lu, Rong and Wu2012; Brand et al., Reference Brand, Fuhrman, Geller, Hattis, Paschall and Tao2004). Future Mars missions, such as a manned Mars landing and Mars base programme, need to achieve a pinpoint landing on Mars to within tens of metres to 100 m of a pre-selected target site. Therefore, high-precision entry navigation and active aerodynamic lift control are required (Martin-Mur et al., Reference Martin-Mur., Kruizinga, Burkhart, Wong and Abilleira2012; Li and Zhang, Reference Li and Zhang2009; Chu, Reference Chu2006).
The main difficulty in achieving high-precision Mars entry navigation comes from the following two aspects. Firstly, the entry dynamics models and relevant parameters usually include larger uncertainties and errors, which greatly degrades the performance of a navigation filter. The most significant sources leading to a larger Mars entry dispersion include the uncertainties in the atmospheric density and aerodynamic coefficients, the accumulated state estimation errors at the atmospheric entry point, and the random winds and gusts (Wolf et al., Reference Wolf, Graves, Powell and Johnson2005; Lévesque, Reference Lévesque2006). Secondly, the available navigation sensors are extremely limited in the Mars atmospheric entry phase due to the existence of a heat shield against an extremely adverse thermal environment. The IMU is almost the only sensor that has been available for Mars entry navigation for many years (Li and Peng, Reference Li and Peng2011; Williams et al., Reference Williams, Menon and Demcak2012).
There have been several new Mars entry navigation concepts and algorithms published in the last decade. To reduce the adverse effects of uncertainties in the Martian atmospheric density, the adaptive sigma point Kalman filter bank, a hierarchical mixture of experts' architecture and multiple model adaptive estimation were adopted to achieve spacecraft precision entry navigation in the presence of a highly dynamic environment with noises and unknown forces (Heyne and Bishop, Reference Heyne and Bishop2006; Ely et al., Reference Ely, Bishop and Dubois-Matra2001; Zanetti and Bishop, Reference Zanetti and Bishop2007). One common character of the research works mentioned above is that only IMU (accelerometer and gyro) data is considered as external measurements and processed by the navigation filter. However, the inertial constant biases and drifts cannot be completely removed and are incorporated into navigation observations instead, which thus degrades the performance of the navigation filter. The performance of dead-reckoning navigation alone usually degrades with time due to the inertial constant biases and drifts. One feasible solution is utilizing external measurements, such as computer vision and Light Direction and Ranging (LIDAR), to correct these inertial biases and drifts, and then improve the landing navigation accuracy (Li et al., Reference Li, Cui and Cui2007; Li et al., Reference Li, Peng, Lu, Zhang and Liu2010). However, most external measurements are not available during the Mars atmospheric entry phase due to the existence of a heat shield. Recent research shows that the ionizing plasma around the entry body has little effect on ultra-high frequency (UHF) band (300∼3000 Mhz) radio communication, which can be utilized in real-time to significantly improve the on board state knowledge during the Mars atmospheric entry phase (Williams et al., Reference Williams, Menon and Demcak2012; Burkhart et al., Reference Burkhart, Ely, Duncan, Lightsey, Campbell and Mogensen2005). Li and Peng (Reference Li and Peng2011) have preliminarily discussed the issue of Mars entry navigation using IMU and orbiting/surface radio beacons (Li and Peng, Reference Li and Peng2011). Lévesque and Lafontaine (Reference Lévesque and Lafontaine2007) studied the navigation performance and observability of four measurement scenarios based on radio ranging during Mars entry (Lévesque and Lafontaine, Reference Lévesque and Lafontaine2007). However, neither Li nor Lévesque discussed the issue of high-precision Mars entry navigation under large uncertainties. The navigation measurements were processed using an unscented Kalman filter and an extended Kalman filter respectively, which lack robust adaptive capability and cannot achieve a higher navigation accuracy in the presence of larger state errors and parameter uncertainties (Li and Peng, Reference Li and Peng2011; Lévesque and Lafontaine, Reference Lévesque and Lafontaine2007).
Desensitised optimal control (DOC) methodology, originally proposed by Seywald and Kumar (Reference Seywald and Kumar1996), is an efficient approach to overcome the unfavourable effect of model and parameter uncertainties. The basic concept is to embed the sensitivity penalty into the performance index by weighting approach and then gain the sensitivity reduction at the cost of sacrificing a small part of the original performance. Karlgaard and Shen extended the main characteristic of the DOC approach to the robust filter design problem such that the performance sensitivity of the filters with respect to the model parameter uncertainties can be reduced (Karlgaard and Shen, Reference Karlgaard and Shen2011; Karlgaard and Shen, Reference Karlgaard and Shen2013). Desensitised state estimates were obtained by minimizing a cost function consisting of the posterior covariance matrix trace penalized by a weighted norm of the state estimate error sensitivities. They continued to apply the DOC methodology to desensitised unscented Kalman filtering (DUKF), in which the cost function of the standard unscented Kalman filter (UKF) is augmented to include a penalty on the sensitivities of the posterior state estimates with respect to the uncertain parameters (Shen and Karlgaard, Reference Shen and Karlgaard2012).
The purpose of this paper is to develop a high-precision integrated navigation algorithm under large uncertainties for Mars atmospheric entry based on a desensitised extended Kalman filter. In order to completely and accurately describe the state variables of an entry vehicle, a new six degree-of-freedom (6-DOF) Mars entry dynamics model is derived based on the angular velocity outputs of a gyro, in which the attitude dynamics are replaced by the outputs of a gyro free of modelling errors. Both the accelerometer outputs and radio measurement information between the entry vehicle and the orbiting beacons are utilized as the observations embedded into a navigation filter to perform state estimation and suppress the measurement noise. At the same time, a desensitised extended Kalman filter, exhibiting the desirable characteristics of efficiently reducing the sensitivity with respect to model and parameter uncertainties, is adopted to overcome the adverse impacts of initial state errors and parameter uncertainties during Mars atmospheric entry and further improve the entry navigation accuracy.
The rest of this paper is organized as follows. Section 2 defines a new 6-DOF Mars entry dynamics model. Section 3 introduces the navigation measurement models used in the subsequent section of the integrated navigation algorithm. The desensitised extended Kalman filter and integrated navigation algorithm are developed at length in Section 4. In Section 5, parameter settings for numerical simulation are defined and simulation results are discussed in detail. Finally, the conclusions and suggestions regarding future research are summarized in Section 6.
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:
$$\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,ry,rz]T is the position vector from the centre of Mars to the vehicle's centre of mass, v=[v x,vy,vz]T is the velocity vector of the entry vehicle, e=[φ,ϑ,φ]T is the tri-axial attitude angle. TVI is the coordinate transformation matrix from the velocity coordinate system to the Mars centred inertial coordinate system, TGI is the coordinate transformation matrix from the geographic coordinate system to the Mars centred inertial coordinate system. The definition of coordinate transformation matrices TVI and TVI can be found in Li and Peng (Reference Li and Peng2011), and is thus not repeated here. aV is the aerodynamic acceleration and can be easily constructed according to Equation (5). gG 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=[ω 1,ω 2,ω 3]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.
$${\bf g}_M = \left[ {0,0, - \displaystyle{{\mu _M} \over {\left| {\bf r} \right|^2}}} \right]^T $$
$${\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]$$
$${\bf \tilde \omega} _B = {\bf \omega} _B + {\bf b}_\omega + {\bf \xi} _\omega $$
$${\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
$$L = \displaystyle{1 \over 2}\rho v^2 \displaystyle{{C_L S} \over m}$$
$$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:
$$\rho = \rho _0 \exp ( - \beta (h))$$where ρ 0=559·351005946503/c 1c 2, c 1=188·95110711075, c 2=1·4×10−13h 3−8·85×10−9h 2−1·245×10−3h+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:
$$\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:
$${\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, aB is the true linear acceleration, ba is the acceleration bias, ξa is the white Gaussian acceleration output noise.
Then, accelerometer measurement model is defined as follows:
$${\bf y}_1 = {\bf \tilde a}_B = T_V^B {\bf a}_V^ * + {\bf b}_a + {\bf \xi} _a $$where
$${\bf a}_V^ * = [\matrix{ { - D^*} & { - L^* \sin \sigma} & {L^* \cos \sigma} \cr} ]^{\rm T} $$
$$L^* = \displaystyle{1 \over 2}\rho ^* v^2 \displaystyle{{C_L S} \over m}$$
$$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
$$\displaystyle{{d{\bf r}_o} \over {dt}} = {\bf v}_o $$
$$\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:
$$\tilde R_i = R_i + \zeta _{Ri} $$
$$R_i = \sqrt {({\bf r} - {\bf r}_i )^T ({\bf r} - {\bf r}_i )} $$where r is the position vector of entry vehicle, ri is the position vector of the ith 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
$$\tilde V_i = V_i + \zeta _{Vi} $$
$$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, vi is the velocity vector of the ith 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:
$${\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 1, …, R m]T, V=[V 1, …, 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:
$${\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:
$${\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:
$${\bf \dot X}_k = {\bf f}{\rm (}{\bf X}_{k - 1} {\rm )} + {\bf \Gamma w}_{k - 1} $$
$${\bf y}{\rm (}t_k {\rm )} = {\bf h}{\rm (}{\bf X}_k {\rm )} + {\bf \upsilon} _k $$where Γ is the model noise input matrix. wk 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
$$p({\bf w})\sim N(0,{\bf Q})$$
$$p({\bf \upsilon} )\sim N(0,{\bf R})$$
$${\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
$${\bf \Phi} _{k/k - 1} = {\bf I} + {\bf F}_{k/k - 1} \cdot \Delta t$$
$${\bf \dot \Phi} (t,\tau ) = {\bf F}(t){\bf \Phi} (t,\tau )$$
$${\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
$${\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 Hk
$${\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
$${\bf \hat X}_0 = E[{\bf X}_0 ]$$
$${\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:
$${\bf \hat X\dot}_{k/k - 1} = {\bf f}(k,{\bf \hat X}_{k - 1} )$$
$${\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:
$${\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]$$
$${\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, yk is the observed quantity, Kk is the filtering gain, Pk / k is the error variance matrix, Φk/k− 1 is the state transition matrix. Qk is the process noise variance matrix, defined in Equation (40), and Rk is measurement noise variance matrix.
$${\bf Q}_k = {\bf Q}(t)T$$ The key of the DEKF is to derive the desensitised optimal gain matrix Kk 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):
$$\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):
$$\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 (∂ yk/∂ρ)=0 in this formulation. At the same time, the value of (∂ Kk/∂ρ) is assumed to be zero here. As any non-zero value of (∂ Kk/∂ρ) 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 Kk, 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
$$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 Kk term by term yields
$$\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/∂ Kk)=0, then the desensitised optimal gain matrix can be solved as follows
$${\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 Sk=HkPk/k− 1HkT+Rk.
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 ba and gyro bias bω applying for the simulation are assumed to be [3×10−3, 3×10−3, 3×10−3] m/s 2 and [5×10−6, 5×10−6, 5×10−6] rad/s respectively, the covariance matrix of white Gaussian noise ξa and ξω are set to [5×10−8, 5×10−8, 5×10−8] m 2/s 4 and [9×10−8, 9×10−8, 9×10−8] (rad/s)2 respectively. Radio range and velocity measurement noise variance in navigation measurement model are assumed to be 100 m2 and 0·1 m2/s 2 respectively. Initial state error covariance matrix P0, system process noise variance matrix Q and penalty matrix P0, 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−7I9×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×104 km 3/km 3s 2, Mars rotating angular rate ω M=7·0882×10−5rad/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.
6. CONCLUSIONS
Pinpoint landing capability is the cornerstone for future Mars sample return and Mars base missions. High-precision entry navigation is considered as one of the key technologies required to achieve pinpoint landing. This paper addresses the development of an innovative high-precision Mars entry integrated navigation algorithm based on a desensitizing extended Kalman filter. The accelerometer outputs and the radio measurements from the entry vehicle to the orbiting beacons are adopted as the observations embedded into the navigation filter to simultaneously perform state estimation and suppress the navigation measurement noise. A desensitizing extended Kalman filter is adopted to effectively reduce the adverse impacts of the uncertainty disturbance during Mars atmospheric entry and then further improve the entry navigation accuracy. Numerical simulations show that the DEKF-based integrated navigation algorithm developed in this paper can achieve 100 m position error and 6·5 m/s velocity error in the presence of larger uncertainties, which can meet the navigation requirement of future pinpoint Mars landing missions.
In our current work, an approximation (∂ Kk/∂ρ) is assumed to derive the desensitised optimal gain matrix, which inevitably leads to an approximate sensitivity rather than the true sensitivity. How to determinate the true sensitivity without this approximation is still an open problem. At the same time, preliminary simulation analysis shows that the relative geometry between the entry vehicle and orbiting radio beacons has a significant impact on the integrated navigation accuracy. These two issues lie beyond the scope of the present paper, though, and will be left for future work.
ACKNOWLEDGEMENTS
The work described in this paper was supported by the National Natural Science Foundation of China (Grant No. 61273051 and 60804057), National High Technology Research and Development Program of China (Grant No. 2012AA121601 and 2011AA7034057E), Qing Lan Project of Jiangsu Province, Innovation Funded Project of Shanghai Aerospace Science and Technology (Grant No. SAST201213), Fundamental Research Funds for the Central Universities, Foundation of Graduate Innovation Center in NUAA (Grant No. kfjj130135).
