2. MARS ENTRY DYNAMIC MODEL AND NAVIGATION SCHEME

2.1. Mars entry dynamic model

The Mars entry phase is defined as the period from when the vehicle enters the Mars atmosphere interface (defined as a radius of 3,552 km from the centre of Mars, and an altitude of 125 km over the surface) to parachute deployment. During the Mars entry process, the entry vehicle will experience a complex and diverse environment. The dynamic model describes the state variables of the entry vehicle, which is required for state prediction. The traditional simplified 3-DOF dynamic model of Mars atmospheric entry is built based on the assumption that the vehicle is regarded as a particle subjected to gravity and aerodynamic forces during the Mars entry phase with the Mars rotation and winds being neglected. Thus, the dynamic model defined with respect to the Mars-centred, Mars-fixed coordinate system (Levesque, Reference Levesque2006) is given by:

(1)$$\eqalign{&{\dot r} = v \sin \gamma \cr & {\dot \theta } = \displaystyle{{v \cos \gamma \sin \psi} \over {r \cos \lambda}} \cr & {\dot \lambda} = \displaystyle{{v \cos \gamma \cos \psi}\over{r}} \cr & \dot {v} = -D-g(r)\sin \gamma \cr & \dot {\gamma} = \left({\displaystyle{{v}\over{r}}-\displaystyle{{g(r)}\over{v}}}\right) \cos \gamma +\displaystyle{{L}\over{v}}\cos \phi \cr & \dot {\psi } = \displaystyle{{v}\over{r}}\sin \psi \cos \gamma \tan \lambda +\displaystyle{{L}\over{v\cos \gamma}}\sin \phi}$$

where r denotes the distance from the centre of Mars to the centre of mass of the entry vehicle, θ and λ denote the longitude and the latitude, respectively, v is the velocity of the entry vehicle, γ is the flight-path-angle (FPA), ψ is the azimuth angle defined as a clockwise rotation angle starting at due north and ϕ is the bank angle defined as the angle about the velocity vector from the local vertical plane to the lift vector, which is considered as the only control variable. To be convenient for simulation, an inverse square gravitational acceleration is used in the model of Mars entry dynamics

(2)$$g(r) = \displaystyle{{\mu}\over{r^2}},\quad \mu =42828.29 \times 10^9\,\hbox{m}^3/{\hbox{s}^2}$$

where μ=GM _Mars is the Mars gravitational constant. D and L are the aerodynamic drag and lift accelerations, which are defined by:

(3)$$D=\displaystyle{P{1}\over{2}} \rho v^2\displaystyle{{C_D S}\over{m}} \quad \quad L = \displaystyle{{1}\over{2}} \rho v^2\displaystyle{{C_L S}\over{m}}$$

where C_D and C_L are the aerodynamic drag and lift coefficients, respectively, which are the functions referring to the angle of attack, the sideslip angle and the Math number, S is the reference surface area, m is the mass of the entry vehicle, L=L/D·D, L/D is the lift-to-drag ratio, B=C _D S/m is called the ballistic coefficient, $\displaystyle{{1}\over{2}}\rho v^{2}$ denotes the dynamic pressure and ρ is the Mars atmospheric density. A simplified exponential Mars atmospheric density model is assumed as:

(4)$$\rho =\rho_0 \exp \left\{{\displaystyle{{(r_0 -r)}\over{h_s }}}\right\}$$

where ρ₀ is the reference density, r ₀=3,437·2 km is the reference radius of Mars and h _s=7·5 km is the constant atmospheric scale height.

Since the true Mars atmosphere density varies with the seasons and temperature, an accurate Mars atmospheric density is still difficult to determine. Equation (4) is only an approximation of the true atmosphere density. The deviation in this model is believed to follow a normal distribution (Li et al., 2014). Moreover, the uncertainties in the aerodynamic coefficients are also a significant source of Mars entry dispersion. Braun et al. (Reference Braun, Powell, Engelund, Gnoffo, Weilmuenster and Mitcheltree2012) analysed the source of uncertainties that affect the Mars Pathfinder flight model. They illustrated that the lack of knowledge of atmosphere density and the computational uncertainties of the aerodynamics analysis were the source of error; also the uncertainties of the aerodynamic drag coefficient and the atmosphere density were time-varying. Striepe et al. (Reference Striepe, Way, Dwyer and Balaram2006) showed that the parameter variations may follow Gaussian and uniform distribution. In the work of Wang and Xia (Reference Wang and Xia2015), it was assumed that the actual values of the three uncertain parameters were subjected to uniform distribution. When the uncertainties in Mars atmosphere density, ballistic coefficient and lift-to-drag ratio are considered, a navigation filter for reducing the adverse effect of the parameter uncertainties is essential.

2.2. Navigation scheme

As mentioned above, due to the existence of the heat shield, several navigation sensors cannot be used during the Mars entry, only the IMU and radiometric measurement are available, but the accuracy of the IMU and radiometric measurement navigation need to be improved. This paper adopts IMU/radio beacon integrated navigation as the navigation scheme. The IMU has the capacity of measuring the accelerator and rotational angular velocity in real time and gaining the state variables of the target without any additional information, which is the key device to implement autonomous navigation. Communication between the vehicle and the radio beacons can provide some important navigation information for the vehicle. Furthermore, research shows that integrated navigation is completely observable when there are not less than three radio beacons (Wang and Xia, Reference Wang and Xia2015). The measurement equations of the vehicle are given by (Li and Peng, Reference Li and Peng2011)

(5)$${\bi z}_k = \left[\matrix{{{\bi a}^m} \cr {\bi R}}\right]+{\bi v}_k$$

where v_k is the zero-mean, white noise with the covariance R_k, a^m is the output vector of accelerometer in the Mars-centred, Mars-fixed coordinate system and R is the radial range between the entry vehicle and the orbiters or the surface beacons.

2.2.1. IMU outputs

An IMU consists of an accelerometer that is used for measuring the non-gravitational linear acceleration along three orthogonal axes and a gyroscope for the rotation angular velocity measurement. The IMU coordinate system is assumed to coincide with the body-fixed coordinate system, thus the acceleration measurement in the body-fixed coordinate system can be given by an accelerometer as follows

(6)$${\tilde{\bi a}}^{b}=\lpar {{\bi I}+{\bi \Gamma}_a}\rpar \lpar {{\bi I}+{\bi S}_a}\rpar \left( {{\bi a}^b +{\bi b}_a +{\bi \xi}_a}\right) $$

(7)$${\bi \Gamma}_a@\left[{\matrix{ 0 & {\gamma_{xz}} & { - \gamma_{xy}} \cr { - \gamma_{yz}} & 0 & {\gamma_{yx}} \cr {\gamma_{zy}} & { - \gamma_{zx}} & 0 \cr } } \right],\quad {\bi S}_a@\left[ {\matrix{ {S_x} & 0 & 0 \cr 0 & {S_y} & 0 \cr 0 & 0 & {S_z} \cr } }\right]$$

where Γ_a is the misalignment/non-orthogonal error matrix, γ_xz, γ_xy, γ_yz, γ_yx, γ_zy, γ_zx are the non-orthogonal and axes misalignment errors, S_a is the scale factor error matrix of the accelerometer, S _x, S _y, S _z are the scale factor errors, a^b denotes the true acceleration along the body-fixed coordinate system, b_a is the unknown acceleration bias and ξ_a is the accelerator stochastic noise, which is modelled as zero-mean, Gaussian-distributed white noise. When the error matrices are small, the acceleration can be formulated as:

(8)$${\tilde{\bi a}}^{b} ={\bi a}^b +\lpar {{\bi \Gamma}_a +{\bi S}_a + {\bi \Gamma}_a {\bi S}_a}\rpar {\bi a}^b +\lpar {{\bi I}+{\bi \Gamma}_a}\rpar \lpar {{\bi I}+{\bi S}_a}\rpar {\bi b}_a +\lpar {{\bi I}+{\bi \Gamma}_a}\rpar \lpar {{\bi I}+{\bi S}_a}\rpar {\bi \xi}_a$$

where ${\bi a}^{b}={\bi T}_{v}^{b} {\bi a}^{v}$, ${\bi a}^{v}=\left[\matrix{ {-D,} & {-D\cdot L / D\cdot \sin \phi}, & {D \cdot L/D \cdot \cos \phi}}\right]^{\rm T}$ is the measurement acceleration vector in the velocity coordinate system. ${\bi T}_{v}^{b}$ is the coordinate transition matrix for the velocity coordinate system to the body-fixed coordinate system.

To be unified with the dynamic model, it is essential for the measurement vector to transfer from the velocity coordinate system to the Mars-centred, Mars-fixed coordinate system:

(9)$${\bi a}^m = {\bi T}_p^m {\bi T}_v^p {\bi a}^v$$

where a^m is the output vector of the accelerometer in the Mars-centred, Mars-fixed coordinate system, ${\bi T}_{p}^{m}$ is the coordinate transition matrix for the navigation coordinate system to the Mars-centred, Mars-fixed coordinate system, and ${\bi T}_{v}^{p}$ is the transition matrix for the velocity coordinate system to the navigation coordinate system, and are given by

(10)$$\eqalign{ & \, {\bi T}_p^m =\left[{\matrix{{\cos \lambda \cos \theta} & {-\sin \theta} & {-\sin \lambda \cos \theta} \cr {\cos \lambda \sin \theta } & {\cos \theta } & {-\sin \lambda \sin \theta} \cr {\sin \lambda } & 0 & {\cos \lambda}}}\right] \cr &{\bi T}_v^p =\left[{\matrix{{\cos \gamma } & {\sin \gamma} & 0 \cr {-}\sin \gamma \cos \psi & {\cos \gamma \cos \psi } &{-\sin \psi} \cr {-\sin \gamma \sin \psi} & {\cos \gamma \sin \psi } & {\cos \psi } \cr}}\right]}$$

2.2.2. Radio ranging measurement

The vehicle communicates with orbiters or surface beacons via Ultra-High Frequency (UHF) radio. By analysing the range signals, the radial range between the entry vehicle and orbiters or surface beacons can be modelled as follows

(11)$${\tilde{\bi R}}={\bi R} + {\bi \xi}_R$$

where ξ_R is zero-mean, Gaussian-distributed white noise.

(12)$${\bi R} = \sqrt{({\bi r}_l - {\bi r}_b)^{\rm T}({\bi r}_l - {\bi r}_b)}$$

(13)$$r_l =\left[ \matrix{ r\cos \lambda \cos \theta \cr r\cos \lambda \sin \theta \cr r\sin \lambda } \right]$$

r_l denotes the position vector of the entry vehicle, r_b denotes the position vector of the orbiter or the surface beacon and both the vectors mentioned above are defined in the Mars centred, Mars-fixed coordinate system.

3. DERIVATIVE-FREE NONLINEAR VERSION OF EXTENDED RECURSIVE THREE-STEP FILTER (DNERTSF)

According to the nonlinear Mars entry dynamical model and measurement equations with parameter uncertainties, the nonlinear systems can be written in the following general form as

(14)$${\bi x}_{k+1} = {\bi f}_{k} ({\bi x}_k ,{\bi u}_k ,{\bi b}_k)+{\bi w}_k$$

(15)$${\bi y}_{k} = {\bi h}_{k} ({\bi x}_k ,{\bi u}_k ,{\bi b}_k)+{\bi v}_k$$

where f(·) and h(·) are nonlinear state transformation and nonlinear measurement functions, respectively. ${\bi x}_{k + 1} \in {\Re}^n$ denotes the state vector, ${\bi y}_k \in {\Re}^m$ denotes the observation vector and ${\bi b}_k \in {\Re}^p$ is the uncertain parameter vector. The process noise w_k and the measurement noise v_k are zero-mean white random noises with known covariance matrices, assumed to be mutually uncorrelated. The noise covariance matrices satisfy the following formula

(16)$$E\left[{\bi w}_k{\bi w}_j^T \right] = {\bi Q}_k\delta_{k,j}, \quad E\left[{\bi v}_k{\bi v}_j^T \right] = {\bi R}_k\delta _{k,j}$$

where δ_k,j = 0, for k ≠ j; δ_k,j=1 for k = j. The initial state estimate ${\hat{\bi x}}_{0}$ is the unbiased mean of the state ${\bi x}_{0} $ with covariance of P₀, which is independent to the noises w_k and v_k.

In this method, the parameters should be estimated firstly by identifying the parameter uncertainties of which the matrices’ numerical approximation (Hsieh, Reference Hsieh2012) is given by

(17)$${\bi F}_k =\left[ {{\bi F}_k^1 \cdots {\bi F}_k^p}\right] ,{\bi F}_k^i = \displaystyle{{{\bi f}_{k} \left( {{\hat{\bi x}}_k ,{\bi u}_k ,{\hat{\bi b}}_k +\Delta_i \times {\bi e}_i^p }\right) -{\bi f}_{k} \left( {{\hat{\bi x}}_k ,{\bi u}_k ,{\hat{\bi b}}_k} \right) }\over{\Delta_i}}$$

(18)$${\bi G}_k =\lsqb {{\bi G}_k^1 \cdots {\bi G}_k^p } \rsqb ,{\bi G}_k^i =\displaystyle{{{\bi h}_k \left( {\hat{{\bi x}}_k^{-} ,{\bi u}_k ,\hat{{\bi b}}_{k-1} +\Delta_i \times {\bi e}_i^p } \right) -{\bi h}_k \left( {\hat{{\bi x}}_k^{-}, {\bi u}_k ,\hat{{\bi b}}_{k-1} } \right) }\over{\Delta_i}}$$

where Equations (17) and (18) are suitable partial derivative matrices, Δ_i,1 ≤ i ≤ p is a carefully chosen small value and $e_{i}^{p}$ is the i-th column of the identity matrix $e^p \in {\Re}^p$. Then the nonlinear discrete-time system given by Equations (14) and (15) can be written as

(19)$${\bi x}_{k + 1} = {\bi f}_k^* \left( {\bi x}_k,{\bi u}_k,{\hat {\bi b}}_k \right) + {\bi F}_k{\bi b}_k + {\bi w}_k$$

(20)$${\bi y}_{k} = {\bi h}_k^* \left( {\bi x}_k,{\bi u}_k,{\hat {\bi b}}_{k-1} \right) + {\bi G}_k {\bi b}_k + {\bi v}_k$$

where ${\bi f}_{k}^{\ast} \left( {{\bi x}_{k} ,{\bi u}_{k} ,\hat{{\bi b}}_{k} } \right) \approx {\bi f}_{k} \left( {{\bi x}_{k},{\bi u}_{k} ,\hat{{\bi b}}_{k} } \right) -\,{\bi F}_{k} \hat{{\bi b}}_{k} $, and ${\bi h}_{k}^{\ast} \left( {{\bi x}_{k} ,{\bi u}_{k} ,\hat{{\bi b}}_{k-1} } \right) \approx {\bi h}_{k} \left( {{\bi x}_{k} ,{\bi u}_{k} ,\hat{{\bi b}}_{k-1} } \right) - {\bi G}_{k} \hat{{\bi b}}_{k-1}$.

Under the general UKF framework, and by applying the ERTSF (Hsieh, Reference Hsieh2009; Hsieh, Reference Hsieh2012) to Equations (19) and (20), we can obtain the DNERTSF for Mars entry navigation as follows

Initialisation:

(21)$$\hat{{\bi x}}_0 =E({\bi x}_0), \quad {\bi P}_0 = E\lsqb {({\bi x}_0 -\hat{{\bi x}}_0 )({\bi x}_0 -\hat{{\bi x}}_0)^{\rm T}} \rsqb $$

Time update:

(22)$${\bi \chi}_{k-1,i} =\left[{{\hat{\bi x}}_{k-1} {\hat{\bi x}}_{k-1} - {\left({\sqrt{(n + \lambda){\bi P}_{k-1}}}\right)_i} {\hat{\bi x}}_{k-1} + {\left({\sqrt{(n +\lambda){\bi P}_{k-1}}}\right)_{i-n}}}\right]$$

(23)$${\bi {\chi}}_{k/k-1,{\rm i}}^\ast = \left\{{\matrix{{\bi f}_{k-1} \left( {{\chi}_{k-1,i}, {\bi u}_{k-1}, {\hat{\bi b}}_{k-1}}\right) \hfill& i=0,1,2,\cdots 2n,\quad k=1 \cr \hfill {\bi f}_{k-1}^\ast \left( {{\bi{\chi}}_{k-1,i},{\bi u}_{k-1},{\hat{\bi b}}_{k-1}}\right) +{\bi F}_{k-1} {\hat{\bi b}}_{k-1,i} & i=0,1,2,\cdots 2n,\quad k\ge 2}}\right.$$

(24)$${\hat{\bi x}}_k^{-} = \sum \limits_{i=0}^{2n} \omega_i^{(m)} {\bi \chi}_{k/{k-1},i}^{\ast}$$

(25)$${\bi P}_k^{-} =\mathop \sum \limits_{i=0}^{2n} \left\{ \omega_i^{\lpar c \rpar } \left[ {\chi}_{k / {k-1,i}}^\ast-\hat{\bi x}_k^{-} \right] \left[ {\chi}_{k/ {k-1,i}}^\ast -\hat{\bi x}_k^{-} \right]^{\rm T} \right\} +{\bi Q}_{k-1}$$

where χ_{k -1,i} is the i-th sampling point, n is the dimension of the state vector, $\left({\sqrt {{\bi P}_{k-1}}} \right)_{i}$ denotes the i-th column of the matrix P_k−1 square root, λ = α²(n + κ)−n is a scaling parameter, α is a constant determining the spread of the sigma points around the mean value, which is usually set to a small positive value, κ is a constant secondary scaling parameter which is usually set to 0 or 3−n and the weights ω_i (Wan and Merwe, Reference Wan and Merwe2000) are given by

(26)$$\eqalign{&{\omega}_0^{(m)} = \lambda/({n+\lambda}),{\omega}_0^{(c)} = \lambda/(n+\lambda)+(1-\alpha^2+\beta), \cr & \omega_i^{(m)} = \omega_i^{(c)} = \lambda/\{2(n+\lambda)\} \quad i=1,2,\cdots 2n}$$

where β is used to incorporate prior knowledge of the distribution of the mean value.

Input estimation and measurement update:

(27)$${\bi \chi}_{k/k-1,i} = \left[{{\hat{\bi x}}^{-}_k \quad {\hat{\bi x}}^{-}_k - {\left({\sqrt{(n + \lambda){\bi P}_k^{-}}}\right)_i} \quad {\hat{\bi x}}^{-}_k +{\left({\sqrt{(n + \lambda){\bi P}_k^{-}}}\right)}_{i-n}}\right]$$

(28)$${\bi y}_{k/k-1,i} = {\bi h}_k^{\ast} \lpar {{\bi \chi}_{k/k-1,i}}\rpar \quad i=0,1,\cdots ,2n$$

(29)$${\hat{\bi y}}_k^{-} = \sum \limits_{i=0}^{2n} \omega_i^{(m)} {\bi y}_{k/k-1,i}$$

(30)$${\bi P}_{\rm yy} = \sum \limits_{i=0}^{2n} \left\{ {\omega_i^{(c)} \left[ {{\bi y}_{k/k-1,i} - {\hat{\bi y}}^{-}_k }\right] \left[ {{\bi y}_{k/k-1,i} -{\hat{\bi y}}_k^{-}}\right]^{\rm T}} \right\} +{\bi R}_k$$

(31)$${\bi P}_{{\rm xy}} = \sum\limits_{i = 0}^{2n} {\left\{ {\omega_i^{(c)} \left[ {{\bi \chi}_{k/k - 1,i} - \hat {\bi x}_k^ - } \right] {\left[ {{\bi y}_{k/k - 1,i} - \hat {\bi y}_k^ - } \right] }^{\rm T}} \right\} }$$

(32)$${\bi {K}}_k ={\bi P}_{{\rm xy}} {\bi P}_{{\rm yy}}^{-1}$$

(33)$${\bi S}_k =\lsqb {\bi G}_k \quad {\bi {H}}_k {\bi F}_{k-1} \lpar {\bi I}-{\bi G}_{k-1}^+ {\bi G}_{k-1} \rpar \rsqb $$

(34)$${\bi \Gamma}_k =\lsqb {\bi {0}} \quad {\bi F}_{k-1} \lpar {\bi I}-{\bi G}_{k-1}^+ {\bi G}_{k-1} \rpar \rsqb $$

(35)$${\bi S}_k^\ast =\left( {{\bi S}_k^{\rm T} {\bi P}_{{\rm yy}}^{-1} {\bi S}_k } \right) ^+{\bi S}_k^{\rm T} {\bi P}_{{\rm yy}}^{-1}$$

(36)$${\bi {K}}_k^b =\lsqb {{\bi G}_{k-1}^+ {\bi G}_{k-1}} \quad {\bi {0} } \rsqb {\bi S}_k^\ast$$

(37)$$\hat{{\bi b}}_{k,i} = {\bi {K}}_k^b \lpar {{\bi y}_k -{\bi y}_{k/k-1,i} } \rpar $$

(38)$$\hat{{\bi b}}_k =\mathop \sum \limits_{i=0}^{2n} \omega _i^{\lpar m\rpar } \hat{{\bi b}}_{k-1,i}$$

(39)$${\bi {L}}_k ={\bi {K}}_k +\lpar {{\bi \Gamma}_k -{\bi {K}}_k {\bi S}_k } \rpar {\bi S}_k^\ast$$

(40)$$\hat{{\bi x}}_k =\hat{{\bi x}}_k^{-} +{\bi {L}}_k \lpar {{\bi y}_k -\hat{{\bi y}}^{-}_k } \rpar $$

(41)$${\bi P}_k ={\bi P}_k^{-} +{\bi {L}}_k {\bi P}_{{\rm yy}} {\bi {L}}_k^{\rm T} -{\bi {L}}_k {\bi P}_{{\rm xy}}^{\rm T} -{\bi P}_{{\rm xy}} {\bi {L}}_k^{\rm T}$$

where H_k is the Jacobian of ${\bi h}_{k}^{\ast} \left( {{\bi x}_{k} ,{\bi u}_{k} ,\hat{{\bi b}}_{k-1} } \right) ,{\bi {H}}_{k} =\left. {{\partial {\bi h}_{k}^{\ast} \left( {{\bi x}_{k}, {\bi u}_{k}, \hat{{\bi b}}_{k-1} } \right) \over \partial {\bi x}_{k} }} \right|_{{\bi x}_{k} =\hat{{\bi x}}_{k}^{-}} = {\bi P}_{{\rm xy}}^{\rm T} \left( {{\bi P}_{k}^{-}} \right) ^{-1}.$

First, the state estimate $\hat{{\bi x}}_{k-1} $ and parameters estimate $\hat{{\bi b}}_{k-1} $ are assumed to be unbiased, thus the propagated estimate $\hat{{\bi x}}_{k}^{-}$ is unbiased as shown in Equation (42). The errors in the propagated state estimate, state estimate and parameter estimate can be written as

(42)$$\tilde{{\bi x}}_k^{-} ={\bi x}_k -\hat{{\bi x}}_k^{-} ={\bi f}_{\,k-1} \left( {{\bi x}_k ,{\bi u}_k ,\hat{{\bi b}}_k } \right) -\mathop \sum \limits_{i=0}^{2n} \omega _i^{\left( m\right) } {\bi f}_{\,k-1} \left( {\bi {\chi}_{k-1,i} ,{\bi u}_{k-1} ,\hat{{\bi b}}_{k-1} } \right) +\,{\bi F}_{k-1} \tilde{{\bi b}}_{k-1} +{\bi w}_{k-1}$$

(43)$$\eqalign{ \tilde{{\bi b}}_k &={\bi b}_k -\hat{\bi b}_k ={\bi b}_k -{\bi K}_k^b \left( {\bi y}_k -{\bi y}_k^{-} \right) \cr &=\left( {\bi I}-{\bi K}_k^b {\bi G}_k \right) {\bi b}_k -{\bi K}_k^b \left[ {\bi h}_k^\ast \left( {\bi x}_k ,{\bi u}_k ,\hat{\bi b}_{k-1} \right) +{\bi v}_k -\mathop \sum \limits_{i=0}^{2n} \omega_i^{\left( m \right) } {\bi h}_k^\ast \left( {\chi}_{k/k-1,i}, {\bi u}_k, \hat{\bi b}_{k-1} \right) \right]}$$

(44)$$\eqalign{ \tilde{{\bi x}}_k &={\bi x}_k -\hat{{\bi x}}_k ={\bi x}_k -\hat{{\bi x}}_k^{-} -{\bi {L}}_k \left( {{\bi y}_k -{\bi y}_{k/k-1} } \right) \cr &=\tilde{{\bi x}}_k^{-} - {\bi {L}}_k \left[ {{\bi h}_k^\ast \left( {{\bi x}_k, {\bi u}_k, \hat{{\bi b}}_{k-1}} \right) -\mathop \sum \limits_{i=0}^{2n} \omega_i^{\left( m \right) } {\bi h}_k^\ast \left( {{\chi}_{k/k-1,i}, {\bi u}_k ,\hat{{\bi b}}_{k-1} } \right) +{\bi {v}}_k +{\bi G}_k {\bi b}_k } \right]}$$

By defining that ${\bi \Pi}_{k} ={\bi {I}}-{\bi {K}}_{k}^{b} {\bi G}_{k}, {\bi \Xi}_{k} ={\bi {K}}_{k}^{b} {\bi {H}}_{k} {\bi F}_{k-1}$ and ${\bi Xi}_{i} {\bi \Pi}_{i-1} ={\bi {0}}$ hold for every k, and substituting the first-order Taylor approximation of Equations (42)–(44) into Equations (43) and (44), we can obtain

(45)$$\eqalign{ E \lsqb {\tilde{{\bi b}}_k } \rsqb &\approx {\bi \Pi}_k {\bi b}_k - \lpar {\bi \Xi}_k {\bi \Pi}_{k-1} \rpar {\bi b}_{k-1} + {\bi \Xi}_k \lpar {\bi \Xi}_{k-1} {\bi \Pi}_{k-2} \rpar {\bi b}_{k-2} \cr & + \cdots + \lpar {-1} \rpar ^k{\bi \Xi}_k \cdots {\bi \Xi}_2 \lpar {\bi \Xi}_1 {\bi \Pi}_0 \rpar {\bi b}_0 = {\bi \Pi}_k {\bi b}_k}$$

(46)$$E\lsqb {\tilde{{\bi x}}_k } \rsqb \approx \lpar {\bi I}-{\bi L}_k {\bi H}_k \rpar {\bi F}_{k-1} {\bi \Pi}_{k-1} {\bi b}_{k-1} -{\bi {L}}_k {\bi G}_k {\bi b}_k$$

In order to obtain the unbiased state estimate $\hat{{\bi x}}_{k-1}$ and the parameter estimate $\hat{{\bi b}}_{k-1} $, the following constraint should be satisfied.

(47)$${\bi K}_k^b \lsqb {\bi G}_k \quad {\bi H}_k {\bi F}_{k-1} {\bi \Pi}_{k-1} \rsqb = \lsqb {\bi I}-{\bi \Pi}_k \quad {\bi 0} \rsqb $$

(48)$${\bi L}_k \lsqb {\bi G}_k \quad {\bi H}_k {\bi F}_{k-1} {\bi \Pi}_{k-1} \rsqb = \lsqb {\bi 0} \quad {\bi F}_{k-1} {\bi \Pi}_{k-1} \rsqb $$

Then, the estimation error covariance matrices are

(49)$${\bi P}_k^b = E\left[ \tilde{{\bi b}}_k \tilde{\bi b}_k^{\rm T} \right] = {\bi K}_k^b {\bi P}_{\rm yy} \lpar {\bi K}_k^b \rpar ^{\rm T}$$

(50)$${\bi P}_k = E \lsqb \tilde{{\bi x}}_k \tilde{\bi x}_k^{\rm T} \rsqb ={\bi P}_k^{-} + {\bi L}_k {\bi P}_{\rm yy} {\bi L}_k^{\rm T} -{\bi L}_k {\bi P}_{\rm xy}^{\rm T} -{\bi P}_{\rm xy} {\bi L}_k^{\rm T}$$

Finally, the gain matrix ${\bi {K}}_{k}^{b} $ is determined by minimising the trace of estimation error covariance matrix Equation (49) under the constraint Equation (47), and ${\bi {L}}_{k} $ is obtained by minimising the trace of estimation error covariance matrix Equation (50) under the constraint Equation (48). Thus, they have the forms of

(51)$${\bi K}_k^b =\lsqb {\bi G}_{k-1}^+ {\bi G}_{k-1} \quad {\bi 0} \rsqb {\bi S}_k^\ast$$

(52)$${\bi L}_k ={\bi K}_k +\lpar {\bi \Gamma}_k -{\bi K}_k {\bi S}_k \rpar {\bi S}_k^\ast$$

where ${\bi S}_{k} =\lsqb {\bi G}_{k} \quad {\bi H}_{k} {\bi F}_{k-1} \lpar {\bi I}-{\bi G}_{k-1}^+ {\bi G}_{k-1} \rpar \rsqb $, ${\bi S}_{k}^{\ast} =\left( {\bi S}_{k}^{\rm T} {\bi P}_{\rm yy}^{-1} {\bi S}_{k} \right) ^+{\bi S}_{k}^{\rm T} {\bi P}_{\rm yy}^{-1}$ and ${\bi K}_{k} ={\bi P}_{\rm xy} {\bi P}_{\rm yy}^{-1}$.

As shown in Equation (22), a set of carefully chosen sample sigma points are used to represent the state distribution such as in a UKF. These sigma points can completely capture the actual mean and covariance of the state, and when propagating through the nonlinear system, they can also capture the posterior mean and covariance accurately. This is one of the advantages for this method in achieving a more accurate navigation for Mars entry. Another advantage is that the distribution of parameter uncertainties is also represented by these sigma points, which avoids obtaining the uncertainties’ error covariance. Moreover, the Jacobian and Hessian computations are not needed any more, which solves the problem that a NERTSF may encounter in the state estimation of a complex nonlinear system.

4. NUMERICAL SIMULATION AND RESULTS

In this section, numerical simulation and analysis of Mars entry navigation using IMU/radio beacon integrated navigation are carried out to illustrate the advantages of the proposed algorithm based on DNERTSF. In the simulation, atmosphere density, ballistic coefficient and lift-to-drag ratio uncertainties are considered. The numerical simulation results are compared with that of UKF. The actual values of the atmosphere density, ballistic coefficient and lift-to-drag ratio come from the investigation of Li et al. (Reference Li, Jiang and Liu2014b) and Wang and Xia (Reference Wang and Xia2015) by assuming that they follow the normal distribution with a standard deviation of 10%. The state variables and nominal values of the parameters are given in Table 1. As mentioned above, three radio beacons are enough to ensure the observability of the system, thus we assume that three surface beacons are available during Mars entry in our simulation. The initial longitudes and latitudes of the surface beacons are given in Table 2. The planned entry span is assumed to be 400 s, which is equal to the period of the vehicle from the atmospheric surface at 125 km to 10km altitude, and the simulation sample step is set to 1 s. The scaling parameters for UKF and DNERTSF are chosen as α = 0.6, β = 2, κ = 3 – n.

Table 1. Initial state variables and physical parameters of entry vehicle (Fu et al., Reference Fu, Yang, Xiao, Wu and Zhang2015; Levesque, Reference Levesque2006)

Table 2. Longitudes and latitudes of surface beacons (Fu et al., Reference Fu, Yang, Xiao, Wu and Zhang2015; Levesque, Reference Levesque2006)

The actual and estimated state variables of the entry vehicle are plotted in Figure 1. The filter is proposed for estimating the state and the parameters of nonlinear stochastic discrete-time varying systems with parameter uncertainties simultaneously. The uncertainties in the atmosphere density, ballistic coefficient and lift-to-drag ratio are considered, while the actual and estimated parameters are shown in Figures 2 and 3. The state variables estimates shown in Figure 1 exhibit good performances as there are almost no differences from the actual values, and Figures 2 and 3 exhibit only a little deviation in the first few seconds, which may be the reason that the algorithm needs time to track the changes of the parameters. Also, the estimate of the product of atmosphere density and ballistic coefficient ρ×B in Figure 2 and the estimate of lift-to-drag ratio L/D in Figure 3 are very close to the actual values after the first 20 second adjustment. Thus, we can conclude that DNERTSF can estimate the state and the parameters accurately.

Figure 1. Actual and estimated state variables of vehicle during Mars entry.

Figure 2. Actual and estimated product of atmosphere density and ballistic coefficient ρ×B.

Figure 3. Actual and estimated lift-to-drag ratio L/D.

Figure 4 depicts the state variables estimated errors from UKF-based and DNERTSF-based IMU/radio beacon integrated navigation. As shown in Figure 4, the state variables estimated errors of DNERTSF-based IMU/radio beacon integrated navigation are much smaller than those of UKF-based integrated navigation. At the beginning of the simulation, the state estimates of UKF are close to the actual values. Since it lacks the capability of estimating parameters, the state estimates of UKF subsequently deviate sharply from the actual values, resulting in the divergence of the estimates. However, the performance of DNERTSF is better. In the first few seconds as shown in Figure 4, the estimate errors of some state variables are larger than those of UKF as it is adjusting. This phenomenon is consistent with the state estimates of DNERTSF deviating from the actual values over the first 20 seconds. When the high precision parameter and state estimates of DNERTSF are provided, the errors are limited to a very low magnitude until the end of the simulation.

Figure 4. State variables estimated errors from UKF and DNERTSF.

To further investigate the performance of DNERTSF, Figure 5 shows the state variables estimated errors of DNERTSF-based integrated navigation with their 3σ covariance bounds.

Figure 5. State variables estimated errors from DNERTSF.

In all the sub-figures of Figure 5, the estimated errors are well captured by the bounds of the covariance of each state to a small magnitude, which shows that the possibility of divergence has been greatly reduced. Finally, 500 Monte Carlo simulations were carried out and the average is represented by Root Mean Square Error (RMSE). The position RMSE and velocity RMSE are depicted in Figures 6 and 7. The subfigures in Figures 6 and 7 show the partially amplified figures of the position RMSE and velocity RMSE of DNERTSF from 200 to 400 seconds. As shown in Figures 6 and 7, the position and velocity RMSE of UKF are much larger than that of DNERTSF. It can be seen from the subfigures that the position and velocity RMSE of DNERTSF are less than 50 m and 5 m/s, respectively, which indicates that this result can fulfil the navigation requirement for future Mars pinpoint landing missions. Even if UKF is adopted for the state estimation with the same simulation condition, it can only approach errors of 1200 m and 50 m/s for position and velocity, respectively. However, the calculation time of DNERTSF is rather greater. In the simulation, UKF only takes 0.7419s but DNERTSF takes 0.9681s. This is because the UKF has not considered the parameter uncertainties and thus saved time in this step of estimation. In all, DNERTSF shows its advantages by improving the accuracy significantly with a small sacrifice in calculation time in dealing with the situation when the system contains parameter uncertainties.

Figure 6. Position RMSE of UKF and DNERTSF.

Figure 7. Velocity RMSE from UKF and DNERTSF.

For the entry vehicle, radio blackout typically lasts a couple of dozen seconds during which IMU/radio beacon integrated navigation is taken as the main navigation scheme. Thus, the radio blackout is considered in our simulation. As Mars Science Laboratory (MSL) experienced an approximate 70 s blackout during the Mars entry phase (Morabito et al., Reference Morabito, Schratz, Bruvold, Ilott, Edquist and Cianciolo2014), we assume that radio blackout phenomenon appears from 150 s to 250 s, for a total of 100 s. The state variables estimated errors of UKF-based and DNERTSF-based integrated navigation with 100 s blackout are shown in Figure 8. This shows the same trend as illustrated in Figure 4 (without the radio blackout). UKF still cannot estimate the state variables accurately, while DNERTSF can track these changes well after a short time adjustment at the beginning. After the first 150 seconds, the IMU and radio ranging measurements can provide effective navigation information for the vehicle. The curves of UKF and DNERTSF are compared in Figure 8. UKF maintains an increasing deviation, while DNERTSF tracks these changes well after a short adjustment period at the beginning. When the blackout happens, the state estimates of DNERTSF deviate a little from the actual values, but they are also captured by their 3σ bounds as shown in Figure 9. We can see from Figure 9 that the state variables estimated errors and their 3σ bounds increase gradually when the radio blackout begins, especially the estimated errors of the longitude and latitude. However, the errors and the 3σ bounds will reduce rapidly to normal levels when the radio blackout ends and radiometric measurements can provide effective ranging information for the vehicle again. Therefore, the method proposed in this paper is suitable for Mars entry navigation even with a radio blackout. It can be concluded that the final navigation accuracy is not decreased when a 100 s radio blackout occurs during Mars entry.

Figure 8. State variables estimated errors from UKF and DNERTSF with radio blackout 150–250 s.

Figure 9. State variables estimated errors from DNERTSF with radio blackout 150–250 s.

In all, the DNERTSF algorithm in this paper can reduce the adverse effects of parameter uncertainties for a nonlinear system effectively, and obtain accurate estimates for Mars entry autonomous navigation, showing good prospects for pinpoint landing.