1. Introduction
Autonomous navigation has become more and more urgent with the need to determine the orbit of interplanetary spacecraft (Sanderson, Reference Sanderson2010; Thisdell, Reference Thisdell2013; Xin et al., Reference Xin, Fang and Ning2013; Karimi and Mortari, Reference Karimi and Mortari2015; Carreau, Reference Carreau2016; Steffes and Barton, Reference Steffes and Barton2017; Jiang et al., Reference Jiang, Yang and Li2018). The flight course of Mars probes usually consists of three stages, namely the cruising stage, the Mars-capturing stage and the Mars-orbiting stage. For autonomous navigation at the Mars-orbiting stage, satellite-to-satellite tracking (SST) in addition to other measurements can usually be observed for some time span to estimate the constellation navigation solution effectively (Liu and Liu Reference Liu and Liu2001; Hill et al., Reference Hill, Lot and Born2006; Hill and Born, Reference Hill and Born2007; Hill and Born, Reference Hill and Born2008; Ma et al., Reference Ma, Jiang and Baoyin2015). A Kalman-like filter is the typical onboard computation scheme, wherein the dynamic models described by the state set consisting of position and velocity components are conventionally utilised in the Cartesian frame (Hill and Born, Reference Hill and Born2007, Reference Hill and Born2008; Ma et al., Reference Ma, Jiang and Baoyin2015; Ning et al., Reference Ning, Wang and Fang2016). At the cruising stage, a probe orbiting the Sun generally possesses a time span of more than several months. For such a long-day cruising flight, the state set in the Cartesian frame above has weak nonlinearity for its slow orbit variation. The weak nonlinearity is equally valid at the Mars-capturing stage, wherein the orbit period is no more than several hours and the flight path can be linearly approximated well in such a short time span. However, the orbit period of Mars probes at the Mars-orbiting stage often varies from several hours to several days, which leads to some intensive nonlinearity (Doody, Reference Doody2009; Martin-Mur et al., Reference Martin-Mur, Kruizinga, Burkhart, Abilleira, Wong and Kangas2014). In the three sets of navigation scenarios above, the extended Kalman filter (EKF) is implemented for effective calculation, which has to rely on the observation of SST combined with other measurements. Unfortunately, the intensive nonlinearity occurring at the Mars-orbiting stage, especially for low-orbit probes, makes the dynamic models expressed by the position–velocity state set above fail in the EKF implementation under the condition where only SST measurements are available. The nonlinearity is caused by the rapid alteration of position and velocity components in one orbit period. The accumulated linearised errors in a relatively large measurement step can make the EKF computation fail to obtain a high-precision navigation solution. The nonlinearity brought by the dynamic models can however be solved by some nonlinear filter techniques, such as the unscented Kalman filter (UKF) method (Ma et al., Reference Ma, Jiang and Baoyin2015, Reference Ma, Wang, Jiang, Mu and Baoyin2017).
An alternative state selection described by the Kepler orbit elements, i.e. semimajor axes, eccentricity, orbit inclination, ascending nodes, argument of perigee and mean anomaly, can effectively mitigate the intensive nonlinearity occurring in the dynamic models at the Mars-orbiting stage as compared with the position–velocity state expression. For Kepler orbit elements, the reduced nonlinearity benefits from the fact that the mean anomaly is the unique rapidly altering variable compared with the six equally altering position–velocity elements. Consequently, the linearised errors are significantly reduced and thus the navigation solution relies less on the computation steps. Unfortunately, as the eccentricity and orbit inclination approach zero, the singularity problem cannot be avoided in dynamic models which are expressed by the standard Kepler orbit elements. Thus revised dynamic models are proposed in this paper by use of singularity-avoiding orbit elements, which have been applied well to Earth-orbiting spacecraft (Chobotov, Reference Chobotov2002; Tapley et al., Reference Tapley, Schutz and Born2004).
Another issue affecting the autonomous navigation of Mars probes is the state observability of all participating probes when only SST measurements are available. Simple two-body dynamic models, namely considering only the central gravitation of Mars, cannot be utilised to estimate the orbit states of a pair of probes in such SST-alone measurement scenarios (Liu and Liu, Reference Liu and Liu2001). To improve the observability of orbit states effectively, more dynamic constraints have to be added to the two-body problem. Hill and Born (Reference Hill and Born2007) found that a spacecraft, orbiting on the L1 and L2 Lagrange points respectively, can experience significant third-body acceleration and thus the orbit states become observable. Xiong et al. (Reference Xiong, Wei and Liu2014) also utilised the crosslink range measurements between two spacecraft, orbiting the Earth and the Moon respectively, to estimate the absolute orbit states. These two trials illustrate that the orbital state observability mainly depends on the assignable third-body attraction in addition to the SST measurements. But the rigorous requirement of remote inter-spacecraft ranges always makes the SST measurements intractable in the onboard calculation implementation, especially for two Mars-orbiting probes. In this sense, the nonspherical perturbation of Mars is a naturally good substitute for third-body attraction as some dynamic constraints for the following two reasons. On the one hand, the components of tesseral harmonics in the nonspherical perturbation have obvious area-related asymmetry which has some perturbation effects on Mars probes. On the other hand, the Mars tesseral harmonic components are much greater than in the case of Earth (Lemoine et al., Reference Lemoine, Smith and Rowlands2001), and thus can be readily sensed by the SST measurements. Besides, this dynamic constraint relaxes the rigorous range requirement between two Mars probes. This paper focuses on utilisation of the tesseral harmonic perturbation of Mars as well as the SST measurements to obtain a better navigation solution for two Mars probes. Another important issue addressed in this paper is the selection of singularity-avoiding orbit elements to address the problems of smaller eccentricities and inclination in the UKF computation.
This paper is organised as follows. Section 2 presents the dynamic models with singularity-avoiding orbit elements associated with the SST measurement models. Section 3 analyses the tesseral harmonic perturbation which is considered as the main dynamic constraint. Simulation tests and the corresponding navigation solution are reported in Section 4. Conclusions and discussion are drawn in Section 5.
2. Autonomous navigation implementation
2.1 Dynamic models using singularity-avoiding orbit elements
For Mars probes, dynamic models are usually established in the J2000.0 Mars-centre Mars-equator orientating frame. These dynamic models are expressed as
\begin{equation}\ddot{\boldsymbol{r}} = {\boldsymbol{a}_0} + {\boldsymbol{a}_{\textrm{ns}}} + {\boldsymbol{a}_{\textrm{nbody}}} + {\boldsymbol{a}_{\textrm{srp}}} + {\boldsymbol{a}_{\textrm{drag}}} + {\boldsymbol{a}_{\varepsilon}}\end{equation}
where 
${\boldsymbol{a}_0}$ denotes the central gravitation of Mars, 
${\boldsymbol{a}_{\textrm{ns}}}$ denotes the nonspherical perturbation of Mars, 
${\boldsymbol{a}_{\textrm{nbody}}}$ denotes the third-body gravitation perturbation from the Sun and other major planets, 
${\boldsymbol{a}_{\textrm{srp}}}$ denotes the solar radiation pressure perturbation, 
${\boldsymbol{a}_{\textrm{drag}}}$ denotes the atmosphere drag perturbation especially for the low-orbit Mars probe, and 
${\boldsymbol{a}_{\varepsilon}}$ denotes the other much smaller perturbation which can be ignored. For more detailed expression of the different gravitation and perturbation, see Ma et al. (Reference Ma, Jiang and Baoyin2015).
In the dynamic models above, the state set is obviously selected as the position and velocity components. However, the intensive nonlinearity manifested by all the six elements above, readily makes the linearised errors of dynamic models calculated by the conventional EKF computation exceed the weak perturbation of tesseral harmonics or the third-body gravitation. By this consideration, the state set can be selected as the Kepler orbit elements and some other combination of Kepler orbit elements. The classical Kepler orbit elements are defined by six elements, i.e. semimajor axes, eccentricity, orbit inclination, right ascending nodes, argument of perigee and mean anomaly, respectively denoted by the following state vector (Chobotov, Reference Chobotov2002; Tapley et al., Reference Tapley, Schutz and Born2004).
\begin{equation}{\boldsymbol{\sigma }_0}:a,e,i,\Omega ,\omega ,M\end{equation}These six elements cannot describe the exact orbit states well, especially for the two following singularity problems. Firstly, the argument of perigee is ambiguous for near-circular orbits with smaller eccentricities. Secondly, the ascending nodes also become inconclusive as the orbit inclination approaches zero. Considering these two kinds of singularities, the state selection of classical Kepler orbit elements cannot be effectively utilised in general to describe the dynamic models.
An alternative expression is the selection of singularity-avoiding orbit elements, which has been well applied to the approximating-circle orbits for Earth-orbiting spacecraft. Singularity-avoiding orbit elements are defined as follows (Liu, Reference Liu1977).
\begin{equation}{\boldsymbol{\sigma }_1}:a,i,\Omega ,\xi = e\cos \omega ,\eta = e\sin \omega ,\lambda = \omega + M\end{equation}The state expression in Equation (3), which avoids the singular eccentricity problem, can also be equivalently used to describe the dynamic models of Mars probes. Besides, the singular orbit inclination problem can be successfully avoided by selecting some reference plane other than the conventional Mars-equator-orientating plane, such as the Earth-equator-orientating plane. Even as the orbit plane of Mars probes approaches the Earth-equator plane, the reference plane can be subsequently altered by selecting the Mars-centred Mars-equator-orientating plane or the Mars-centred ecliptic plane instead of that in the singular cases.
The dynamic models for each Mars probe, described by the singularity-avoiding orbit elements, are established as follows
\begin{equation}{\dot{\boldsymbol{\sigma }}_{1,k}} = f({{\boldsymbol{\sigma }_{1,k}},t} ),({k = 1,2} )\end{equation} The explicit expression for the 
$k - \textrm{th}$ probe is given as the following Gaussian perturbation equation (Liu, Reference Liu1977).
\begin{align}&\left\{ \begin{array}{l} \dfrac{{d{a_k}}}{{dt}} = \dfrac{2}{{n\sqrt {1 - e_k^2} }}[{{a_{k,Z}}({{\xi_k}\sin {u_k} + {\eta_k}\cos {u_k}} )+ {a_{k,X}}({1 + {e_k}\cos {f_k}} )} ]\\ \dfrac{{d{i_k}}}{{dt}} = \dfrac{{{r_k}\cos {u_k}}}{{{n_k}a_k^2\sqrt {1 - e_k^2} }}{a_{k,Y}}\\ \dfrac{{d{\Omega _k}}}{{dt}} = \dfrac{{{r_k}\sin {u_k}}}{{{n_k}a_k^2\sqrt {1 - e_k^2} \sin {i_k}}}{a_{k,Y}} \end{array} \right.\end{align}
\begin{align}&\left\{ \begin{array}{l} \dfrac{{d{\xi_k}}}{{dt}} = \dfrac{{\sqrt {1 - e_k^2} }}{{{n_k}{a_k}}}\\ \quad \left\{ {{a_{k,Z}}\sin {u_k} + \dfrac{{{a_{k,X}}}}{{\sqrt {1 - e_k^2} }}\left[ {\left( {1 + \sqrt {1 - e_k^2} } \right)\cos {u_k} - \dfrac{{{\xi_k}}}{{1 + \sqrt {1 - e_k^2} }}({{\xi_k}\cos {{\tilde{u}}_k} + {\eta_k}\sin {{\tilde{u}}_k}} )} \right]} \right\}\\ \qquad +\, \eta_{k}\cos {i_k}\dfrac{{d{\Omega _k}}}{{dt}}\\ \dfrac{{d{\eta_k}}}{{dt}} = \dfrac{{\sqrt {1 - e_k^2} }}{{{n_k}{a_k}}}\\ \quad \left\{ {{a_{k,Z}}\cos {u_k} + \dfrac{{{a_{k,X}}}}{{\sqrt {1 - e_k^2} }}\left[ { - \left( {1 + \sqrt {1 - e_k^2} } \right)\sin {u_k} + \dfrac{{{\eta_k}}}{{1 + \sqrt {1 - e_k^2} }}({{\xi_k}\cos {{\tilde{u}}_k} + {\eta_k}\sin {{\tilde{u}}_k}} )} \right]} \right\}\\ \qquad +\, \xi_{k}\cos {i_k}\dfrac{{d{\Omega _k}}}{{dt}}\\ \dfrac{{d{\lambda_k}}}{{dt}} = {n_k} - \dfrac{{\sqrt {1 - e_k^2} }}{{{n_k}{a_k}}}\Bigg\{ {2{a_{k,Z}}\sqrt {1 - e_k^2} \left( {\dfrac{{{r_k}}}{{{p_k}}}} \right) + \dfrac{1}{{1 + \sqrt {1 - e_k^2} }}\Bigg[{{a_{k,Z}}({{\xi_k}\cos {u_k} + {\eta_k}\sin {u_k}} )} } \\ \qquad { { -\, {a_{k,X}}\left( {1 + \dfrac{{{r_k}}}{{{p_k}}}} \right)({{\xi_k}\sin {u_k} - {\eta_k}\cos {u_k}} )} \Bigg]} \Bigg\} - \cos {i_k}\dfrac{{d{\Omega _k}}}{{dt}} \end{array} \right.\end{align}
where 
${a_{k,X}}$ denotes the in-plane forward perturbation, 
${a_{k,Y}}$ denotes the normal-plane perturbation and 
${a_{k,Z}}$ denotes the in-plane radial perturbation. Correspondingly, the orbital frame is defined with each axis pointing to the direction of perturbation acceleration above.
For the 
$k - \textrm{th}$ Mars probe, the position vector, velocity vector and acceleration vector are denoted respectively by 
${\boldsymbol{r}_k}$, 
${\dot{\boldsymbol{r}}_k}$ and 
${\ddot{\boldsymbol{r}}_k}$. The absolute acceleration of Mars probes in the inertial frame equals the sum of centre-body gravitation and all of the perturbation acceleration as follows
\begin{equation}{\ddot{\boldsymbol{r}}_k} ={-} \frac{{{\mu _{\textrm{MARS}}}}}{{\boldsymbol{r}_k^3}} + \boldsymbol{C}_k^T{\left[ {\begin{array}{*{20}{c}} {{a_{k,X}}}&{{a_{k,Y}}}&{{a_{k,Z}}} \end{array}} \right]^T}\end{equation}Thus the perturbation acceleration, which needs to be substituted into the dynamic models, can be readily obtained as
\begin{equation}{\left[ {\begin{array}{*{20}{c}} {{a_{k,X}}}&{{a_{k,Y}}}&{{a_{k,Z}}} \end{array}} \right]^T} = {\boldsymbol{C}_k}\left( {{{\ddot{\boldsymbol{r}}}_k} + \frac{{{\mu_{\textrm{MARS}}}}}{{\boldsymbol{r}_k^3}}} \right)\end{equation}
where the direction cosine matrix 
${\boldsymbol{C}_k}$ transforms the perturbation acceleration in the inertial frame to that in the orbital frame. The row elements of 
${\boldsymbol{C}_k}$ can be readily computed as
\begin{equation}\left\{ \begin{array}{l} {\boldsymbol{C}_k}(3,:) ={-} {{{\boldsymbol{r}_k}} / {|{{\boldsymbol{r}_k}} |}}\\ {\boldsymbol{C}_k}(2,:) ={-} {{{\boldsymbol{r}_k} \times {{\dot{\boldsymbol{r}}}_k}} / {|{\boldsymbol{r}_k} \times {{\dot{\boldsymbol{r}}}_k}|}}\\ {\boldsymbol{C}_k}(1,:) = {\boldsymbol{C}_k}(2,:) \times {\boldsymbol{C}_k}(3,:) \end{array} \right.\end{equation}For more details about other notations in Equations (5.1) and (5.2), readers are recommended to refer to the explanation by Liu (Reference Liu1977).
2.2 Observation models based on SST measurements
The SST measurements are usually selected as the relative range and Doppler measurements (
$\rho$ and 
$\dot{\rho }$) between two Mars probes. For the SST measurements, the line of sight between the two probes may be obscured by Mars in some time span. Hence, the SST measurements can only be utilised in the visible observation segments. The observation equation is readily given as
\begin{equation}\boldsymbol{y} = {\left[ {\begin{array}{*{20}{c}} \rho &{\dot{\rho }} \end{array}} \right]^T} = h({{\boldsymbol{\sigma }_{1,k}},t} ),\quad ({k = 1,2} )\end{equation}where the relative range and Doppler measurements are described as follows (Ma et al., Reference Ma, Jiang and Baoyin2015)
\begin{align}\rho &= ||{{\boldsymbol{r}_1} - {\boldsymbol{r}_2}} ||= \sqrt {{{({x_1} - {x_2})}^2} + {{({y_1} - {y_2})}^2} + {{({z_1} - {z_2})}^2}}\end{align}
\begin{align}\dot{\rho } &= \frac{{({\boldsymbol{r}_1} - {\boldsymbol{r}_2}) \cdot ({{\dot{\boldsymbol{r}}}_1} - {{\dot{\boldsymbol{r}}}_2})}}{\rho } = \frac{{({x_1} - {x_2})({{\dot{x}}_1} - {{\dot{x}}_2}) + ({y_1} - {y_2})({{\dot{y}}_1} - {{\dot{y}}_2}) + ({z_1} - {z_2})({{\dot{z}}_1} - {{\dot{z}}_2})}}{{\sqrt {{{({x_1} - {x_2})}^2} + {{({y_1} - {y_2})}^2} + {{({z_1} - {z_2})}^2}} }}\end{align} The orbit states for each Mars probe, i.e. 
${\boldsymbol{r}_k}$ and 
${\dot{\boldsymbol{r}}_k}$, are obtained from the conventional transformation of classical Kepler orbit elements of 
${\boldsymbol{\sigma }_{0,k}}$, which can be readily solved from the singularity-avoiding orbit elements of 
${\boldsymbol{\sigma }_{1,k}}$ (Chobotov, Reference Chobotov2002; Tapley et al., Reference Tapley, Schutz and Born2004).
2.3 The UKF estimation method
According to the aforementioned dynamic models and SST measurement models, the UKF estimation method is selected to accomplish the autonomous navigation calculation. The state transition matrix of dynamic models and Jacobi matrix of SST measurement models, which are inevitable requisites in EKF, are not needed in UKF. Instead, the mean and covariance of probability distributions are propagated by some appropriate weighted sample points, named as sigma points, which can be obtained by the unscented transformation (UT) (Julier and Uhlmann, Reference Julier and Uhlmann1997, Reference Julier and Uhlmann2004). The intensive nonlinearities occurring in the dynamic models and SST measurement models can be faithfully described by these sigma points. Thus the UKF estimator can yield an estimation performance equivalent to the KF for linear systems but generalise elegancy to nonlinear systems without the linearisation steps required by the EKF.
The primary calculation scheme in UKF is described as follows.
Step 1 Construct the sigma points by UT as follows.
(12)
\begin{equation}\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{\chi }_0} = {\boldsymbol{x}_j}}&{{w_0} = \dfrac{\kappa }{{n + \kappa }}}\\ {{\boldsymbol{\chi }_i} = {\boldsymbol{x}_j} + {{\left( {\sqrt {({n + \kappa } ){\boldsymbol{P}_{xx}}} } \right)}_i}}&{{w_i} = \dfrac{1}{{2({n + \kappa } )}}}\\ {{\boldsymbol{\chi }_{i + n}} = {\boldsymbol{x}_j} - {{\left( {\sqrt {({n + \kappa } ){\boldsymbol{P}_{xx}}} } \right)}_i}}&{\begin{array}{*{20}{c}} {{w_{i + n}} = \dfrac{1}{{2({n + \kappa } )}}} \end{array}} \end{array}} \right.\end{equation}
Assuming that the state of 
$\boldsymbol{x}$ follows Gaussian distribution, a heuristic selection of 
$n + \kappa = 3$ is always recommended. In the filtering implementation, 
$n = 12$ and 
$\kappa ={-} 9$; 
${w_i}$ denotes the weight factor associated with the 
$i - \textrm{th}$ sigma point.
Step 2 Transform the sigma points through the dynamic models by some numerical integration, such as the RK7(8) integration, as follows.
\begin{equation}{\boldsymbol{\chi }_{i,j + 1|j}} = f({{\boldsymbol{\chi }_i}} )\end{equation}Step 3 Propagate the state estimation and the covariance matrix as follows.
\begin{equation}\begin{cases} {\boldsymbol{x}_{j + 1|j}} = \displaystyle\sum_{i = 0}^{2n} {{w_i}} {\boldsymbol{\chi }_{i,j + 1|j}}\\ {\boldsymbol{P}_{xx,j + 1|j}} = \displaystyle\sum\limits_{i = 0}^{2n} {{w_i}({{\boldsymbol{\chi }_{i,j + 1|j}} - {\boldsymbol{x}_{j + 1|j}}} ){{({{\boldsymbol{\chi }_{i,j + 1|j}} - {\boldsymbol{x}_{j + 1|j}}} )}^T} + {Q_j}} \end{cases}\end{equation}
where 
${Q_j}$ denotes the unmoulded process noise.
Step 4 Predict the observation and the covariance matrix as follows.
\begin{equation}\begin{cases}{\boldsymbol{y}_{j + 1|j}} = \displaystyle\sum_{i = 0}^{2n} {{w_i}{\boldsymbol{\varUpsilon }_{i,j + 1|j}} = \displaystyle\sum\limits_{i = 0}^{2n} {{w_i}h({{\boldsymbol{\chi }_{i,j + 1|j}}} )} } \\ {\boldsymbol{P}_{yy,j + 1|j}} = \displaystyle\sum\limits_{i = 0}^{2n} {{w_i}({{\boldsymbol{\varUpsilon }_{i,j + 1|j}} - {\boldsymbol{y}_{j + 1|j}}} ){{({{\boldsymbol{\varUpsilon }_{i,j + 1|j}} - {\boldsymbol{y}_{j + 1|j}}} )}^T} + {R_{j + 1}}} \end{cases}\end{equation}
where 
${R_{j + 1}}$ denotes the observation noise.
Step 5 Calculate the Kalman gain matrix as follows.
\begin{equation}{\boldsymbol{K}_{j + 1}} = {\boldsymbol{P}_{xy,j + 1|j}}\boldsymbol{P}_{yy,j + 1|j}^{ - 1} = \boldsymbol{P}_{yy,j + 1|j}^{ - 1}\sum\limits_{i = 0}^{2n} {{w_i}({{\boldsymbol{\chi }_{i,j + 1|j}} - {\boldsymbol{x}_{j + 1|j}}} ){{({{\boldsymbol{\varUpsilon }_{i,j + 1|j}} - {\boldsymbol{y}_{j + 1|j}}} )}^T}}\end{equation}Step 6 Update the estimated state and the covariance matrix as follows.
\begin{equation}\left\{ \begin{array}{l} {\boldsymbol{x}_{j + 1}} = {\boldsymbol{x}_{j + 1|j}} + {\boldsymbol{K}_{j + 1}}({{\boldsymbol{y}_{j + 1}} - {\boldsymbol{y}_{j + 1|j}}} )\\ {\boldsymbol{P}_{xx,j + 1}} = {\boldsymbol{P}_{xx,j + 1|j}} - {\boldsymbol{K}_{j + 1}}{\boldsymbol{P}_{yy,j + 1|j}}\boldsymbol{K}_{j + 1}^T \end{array} \right.\end{equation}3. Analysis of Mars nonspherical perturbation
The potential function of the Mars nonspherical perturbation is usually expressed by the following spherical harmonic equation (Chobotov, Reference Chobotov2002; Tapley et al., Reference Tapley, Schutz and Born2004), which is similar to the Earth nonspherical perturbation.
\begin{align} V & = \dfrac{{{\mu _{\textrm{Mars}}}}}{r}\left\{ {\sum\limits_{n = 2}^\infty {{{\bar{C}}_{n,0}}{{\left( {\dfrac{{{R_{\textrm{Mars}}}}}{r}} \right)}^n}{{\bar{P}}_n}({\sin \varphi } )} } \right.\notag\\ & \quad + \left. {\sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^n {{{\left( {\dfrac{{{R_{\textrm{Mars}}}}}{r}} \right)}^n}} {{\bar{P}}_{nm}}({\sin \varphi } )({{{\bar{C}}_{nm}}\cos m\lambda + {{\bar{S}}_{nm}}\sin m\lambda } )} } \right\} \end{align}
where 
${\bar{P}_n}(\sin \phi )$ and 
${\bar{P}_{nm}}(\sin \phi )$ denote the normalised Legendre and association of Legendre function respectively, 
${\bar{C}_{n,0}}$ denotes the normalised coefficients of zonal harmonics, 
${\bar{C}_{nm}}$ and 
${\bar{S}_{nm}}$ denote the normalised coefficients of tesseral harmonics, 
${R_{\textrm{Mars}}}$ denotes the Mars equatorial radius, 
$\lambda$ and 
$\varphi$ denote the longitude and latitude respectively in the Mars-centred Mars-fixed frame, and r denotes the distance from the centre of Mars to the probe.
According to the potential function in Equation (18), the Mars nonspherical perturbation can be calculated by the following expression (Chobotov, Reference Chobotov2002; Tapley et al., Reference Tapley, Schutz and Born2004).
\begin{equation}{\boldsymbol{a}_{ns}} = \frac{{\partial V}}{{\partial \boldsymbol{r}}} = \left[ \begin{array}{l} \dfrac{{\partial V}}{{\partial x}}\\ \dfrac{{\partial V}}{{\partial y}}\\ \dfrac{{\partial V}}{{\partial z}} \end{array} \right] = \left[ \begin{array}{l} \dfrac{{\partial V}}{{\partial r}} \cdot \dfrac{{\partial r}}{{\partial x}} + \dfrac{{\partial V}}{{\partial \varphi }} \cdot \dfrac{{\partial \varphi }}{{\partial x}} + \dfrac{{\partial V}}{{\partial \lambda }} \cdot \dfrac{{\partial \lambda }}{{\partial x}}\\ \dfrac{{\partial V}}{{\partial r}} \cdot \dfrac{{\partial r}}{{\partial y}} + \dfrac{{\partial V}}{{\partial \varphi }} \cdot \dfrac{{\partial \varphi }}{{\partial y}} + \dfrac{{\partial V}}{{\partial \lambda }} \cdot \dfrac{{\partial \lambda }}{{\partial y}}\\ \dfrac{{\partial V}}{{\partial r}} \cdot \dfrac{{\partial r}}{{\partial z}} + \dfrac{{\partial V}}{{\partial \varphi }} \cdot \dfrac{{\partial \varphi }}{{\partial z}} + \dfrac{{\partial V}}{{\partial \lambda }} \cdot \dfrac{{\partial \lambda }}{{\partial z}} \end{array} \right]\end{equation}The explicit expression for each partial derivative above is respectively given as follows.
\begin{align}\frac{{\partial r}}{{\partial x}} &= \frac{x}{r},\;\frac{{\partial r}}{{\partial y}} = \frac{y}{r},\;\frac{{\partial r}}{{\partial z}} = \frac{z}{r}\end{align}
\begin{align}\frac{{\partial \varphi }}{{\partial x}} &={-} \frac{{xz}}{{{r^2}\sqrt {{x^2} + {y^2}} }},\;\frac{{\partial \varphi }}{{\partial y}} ={-} \frac{{yz}}{{{r^2}\sqrt {{x^2} + {y^2}} }},\;\frac{{\partial \varphi }}{{\partial z}} = \frac{{\sqrt {{x^2} + {y^2}} }}{{{r^2}}}\end{align}
\begin{align}\frac{{\partial \lambda }}{{\partial x}} &= \frac{{ - y}}{{\sqrt {{x^2} + {y^2}} }},\;\frac{{\partial \lambda }}{{\partial y}} = \frac{x}{{\sqrt {{x^2} + {y^2}} }},\;\frac{{\partial \lambda }}{{\partial z}} = 0\end{align}
\begin{align}\frac{{\partial V}}{{\partial r}} &={-} \frac{{{\mu _{\textrm{Mars}}}}}{{{r^2}}}\left[ {\sum\limits_{n = 2}^\infty {({n + 1} ){{\bar{C}}_{n,0}}{{\left( {\frac{{{R_{\textrm{Mars}}}}}{r}} \right)}^n}{{\bar{P}}_n}(\sin \varphi ) + \sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^n {({n + 1} ){{\bar{P}}_{nm}}({\sin \varphi } ){T_{nm}}} } } } \right]\end{align}
\begin{align}\frac{{\partial V}}{{\partial \varphi }} &= \frac{{{\mu _{\textrm{Mars}}}}}{r}\left[ {\sum\limits_{n = 2}^\infty {{{\bar{C}}_{n,0}}{{\left( {\frac{{{R_{\textrm{Mars}}}}}{r}} \right)}^n}\frac{{\partial {{\bar{P}}_n}({\sin \varphi } )}}{{\partial \varphi }} + \sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^n {\frac{{\partial {{\bar{P}}_{nm}}({\sin \varphi } )}}{{\partial \varphi }}{T_{nm}}} } } } \right]\end{align}
\begin{align}\frac{{\partial V}}{{\partial \lambda }} &= \frac{{{\mu _{\textrm{Mars}}}}}{r}\sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^n {{{\bar{P}}_{nm}}({\sin \varphi } )\frac{{\partial {T_{nm}}}}{{\partial \lambda }}} }\end{align}
\begin{align}{T_{nm}} &= {\left( {\frac{{{R_{\textrm{Mars}}}}}{r}} \right)^n}({\bar{C}_{nm}}\cos m\lambda + {\bar{S}_{nm}}\sin m\lambda )\end{align}
\begin{align}\frac{{\partial {T_{nm}}}}{{\partial \lambda }} &= m{\left( {\frac{{{R_{\textrm{Mars}}}}}{r}} \right)^n}[{{{\bar{S}}_{nm}}\cos m\lambda - {{\bar{C}}_{nm}}\sin m\lambda } ]\end{align} It is indicated from Equation (26) that the item of 
${T_{nm}}$ is related to the tesseral harmonics of 
${\bar{C}_{nm}}$ and 
${\bar{S}_{nm}}$, but has no association with the zonal harmonics of 
${\bar{C}_{n,0}}$. As is subsequently shown from Equations (23)–(25), the perturbation acceleration in the distance and latitude dimensions depends on both the zonal and tesseral harmonics, but in the longitude dimension it depends only on the tesseral harmonics. According to this conclusion, no matter how the time span of SST measurements grows, the ascending nodes of the two probes, which are closely related to the longitude dimension, are still indistinguishable by the zonal harmonics perturbation. On the contrary, the effects of tesseral harmonic perturbation on the probes are related to the SST measurement segments. In other words, the absolute orbit states depend closely on the tesseral harmonic perturbation in different areas on Mars. Owing to this dependence, the ascending nodes may be recognised in theory for two probes.
As is well known, the nonspherical perturbation of Earth has the same expression as that of Mars in Equation (19). For Earth-orbiting spacecraft, the perturbation coefficient of 
${\bar{C}_{2,0}}$, which is the principal zonal harmonics perturbation, is no more than 10−3, but the tesseral harmonic perturbation is less than 10−6 (Vallado, Reference Vallado2013). Obviously, the tesseral harmonic perturbation is much smaller than that in the zonal harmonics case, for which the absolute orbit states of two spacecraft can hardly be obtained by using SST measurements alone in a relatively short time span. Fortunately, the Mars tesseral harmonic perturbation is relatively greater than that of the Earth. A widespread gravitational-field model of Mars called the GMM-2B model (Lemoine et al., Reference Lemoine, Smith and Rowlands2001) is listed in Table 1. As can be seen from Table 1, most of the coefficients below the 4-by-4 degree-and-order are greater than 10−5, which makes the tesseral harmonic perturbation have more obvious effects on the low-orbit Mars probes. Besides, for the high-orbit Mars probes, the tesseral harmonic perturbation will become much smaller for greater distances, but the solar gravitational perturbation is still not to be ignored and can be also selected as an additional dynamic constraint relating to its absolute orbit states. For example, the solar gravitational perturbation on the probe, with an orbit altitude of more than 20,000 km, is greater than 10−5, which is the typical magnitude for the central gravitation of Mars. Thus in the third-body problem including both the gravitation of Mars and the Sun, the absolute orbit states are also observable in the prerequisite of SST-alone measurements. In these two senses above, the observability of absolute orbit states is relatively improved for autonomous navigation by Mars probes, no matter which constellation is designed, i.e., one low-orbit probe and one high-orbit probe, two low-orbit probes or two high-orbit probes. Due to this improved observability, an inevitable benefit of obtaining the precise navigation solution by observing the SST measurements in a short time span becomes a prospective manner of autonomous navigation for Mars probes.
Coefficients of the Mars GMM-2B models
Another troublesome issue, inhibiting the implementation of computation of autonomous navigation effectively, is overcoming the nonlinearity manifested by the dynamic models. In fact, in both cases of nonspherical gravitation perturbation and third-body gravitation perturbation, the absolute orbit states are weakly estimable using the SST measurements alone. As for the EKF estimation, the linearised errors of dynamic models may exceed the nonspherical gravitation perturbation or the third-body gravitation perturbation. The resulting failed EKF computation of orbit states is usually solved by turning to its UKF counterpart. As is well acknowledged (VanDyke et al., Reference VanDyke, Schwartz and Hall2004; Daum, Reference Daum2005), only the first-order approximation of nonlinear models can be achieved by the EKF method, whereas the UKF mode approaches the nonlinear dynamic models with the second-order errors still preserved. By this consideration, the precise navigation solution may be obtained by the UKF computation method while only SST measurements between two Mars probes are available.
Additionally, the infinite degrees and orders of nonspherical perturbation models in Equation (18) need to be truncated for the numerical calculation. For the GMM-2B model, the calculation of the most precise position and velocity states can usually be implemented by the 80-by-80 degree-and-order model at most, which can certainly be selected as the referenced prototype. Too small degree-and-order truncated nonspherical perturbation models may degrade the SST measurements. For instance, the predictive position errors of Mars-orbiting probes with an altitude of 300 km in the time span of two days can reach the magnitude of more than 1 km by the 32-by-32 degree-and-order model, which can be clearly seen from Figure 1.
Predictive position and velocity errors by the 32-by-32 GMM-2B model
Other comparison of predictive velocity and position errors with the referenced prototype, respectively by the 50-by-50, 60-by-60 and 70-by-70 degree-and-order truncated models, is also depicted in Figures 2 and 3. The position errors in the time span of two days are less than 50 m, 40 m and 10 m respectively for the corresponding truncated model. The velocity errors in the same time span are less than 0⋅05 m/s, 0⋅04 m/s and 0⋅01 m/s, respectively. In such scenarios, the least 70-by-70-order model has to be selected for the coincidence with the predictive requirement of position errors less than 10 m in two days.
Predictive position errors by higher-order truncated GMM-2B models
Predictive velocity errors by higher-order truncated GMM-2B models
4. Simulation and results
4.1 Simulation scenarios
Simulation scenarios for validating the proposed method of autonomous navigation are described by two Mars-orbiting probes: one low-orbit probe and one high-orbit probe. The orbit states of these two Mars probes need to be fictionalised to calculate the SST measurements. The initial epoch is set at UTC time 2020-01-01 0:00:0.00. For the low-orbit probe, the orbit parameters of 3,696 km, 0⋅005 and 5° are respectively selected for the semimajor axis, the eccentricity and the inclination. As a comparison, the semimajor axis is set at 20,000 km, the eccentricity is set at 0⋅001 and the inclination is set at 75° for the high-orbit Mars probe. The initial orbit states of these two Mars probes are further described in Table 2.
Initial orbit states of two Mars probes
According to the initial orbit states in Table 2, the referenced orbit evolution of the two Mars probes can be calculated by using high-precision dynamic models dominated by the central gravitation of Mars, the 80-by-80 degree-and-order nonspherical perturbation, the third-body gravitation perturbation, the solar radiation pressure perturbation and the drag pressure perturbation. In calculating the orbit evolution, the RKF7(8) numerical integration algorithm is applied to the orbit propagation for 30 days for the considered dynamic models. Based on the orbit evolution of the two Mars probes, the relative range and Doppler measurements, with a sampling rate of 1 min, can be directly calculated as a temporal series. The root mean square errors of 10 m and 0⋅1 m/s are respectively added to the true SST measurements. The selected random errors are equally suitable for the current SST measurements in the real navigation scenarios.
4.2 Visibility analysis
For the SST measurements, the line of sight between the low-orbit and the high-orbit probe may be obscured by Mars in some time spans. SST measurements are certainly generated in the visible observation segments. The invisible segments may lead to some time delay in the convergence of the navigation solution, but the final navigation solution can still be obtained only by the visible observation segments. The visible segments of the SST measurements in a one-day time span are described in Figure 4.
Visible segments in one-day time span
4.3 Navigation results
With SST measurements alone, autonomous navigation based on the UKF computation is respectively accomplished by two kinds of dynamic models. No matter what dynamic models are utilised, the initial position and velocity errors are respectively selected as 5,000 m and 3 m/s for the low-orbit Mars probe, and 20,000 m and 1 m/s for the high-orbit Mars probe. For the first dynamic model described by the position–velocity state and the 80-by-80 degree-and-order nonspherical gravitation perturbation, the navigation errors are shown in Figures 5 and 6. As is clearly illustrated, no divergence occurred in the observed position and velocity, and the navigation errors are almost identical with the initial errors. Thus the UKF method failed for the position–velocity state selection even though the 80-by-80 degree-and-order model was utilised.
Position errors by the 80-by-80 degree-and-order model
Velocity errors by the 80-by-80 degree-and-order model
For the second dynamic model described by the singularity-avoiding orbit elements, the 32-by-32, 50-by-50, 70-by-70 and 80-by-80 degree-and-order nonspherical perturbation are respectively analysed. Due to the randomness of measurement noises, 10 sets of Monte Carlo simulation runs are implemented respectively for the different degree-and-order models above. For a brief description, only the navigation position errors, which are calculated by averaging the daily position errors in the whole time span of 30 days, are graphically shown as follows.
As can be drawn from Figures 7–14, more accurate navigation solutions can be obtained by the higher degree-and-order models. Select one arbitrary set of Monte Carlo simulation runs respectively by different degree-and-order models above as the error comparison. This conclusion can be graphically shown as in Figures 15–18. Additionally, the numerical statistics of navigation errors in the converged time span which is selected from the third to the 30th days are described in Table 3.
Position errors by the 32-by-32 model for the low-orbit Mars probe
Position errors by the 32-by-32 model for the high-orbit Mars probe
Position errors by the 50-by-50 model for the low-orbit Mars probe
Position errors by the 50-by-50 model for the high-orbit Mars probe
Position errors by the 70-by-70 model for the low-orbit Mars probe
Position errors by the 70-by-70 model for the high-orbit Mars probe
Position errors by the 80-by-80 model for the low-orbit Mars probe
Position errors by the 80-by-80 model for the high-orbit Mars probe
Position errors by different models for the low-orbit Mars probe
Velocity errors by different models for the low-orbit Mars probe
Position errors by different models for the high-orbit Mars probe
Velocity errors by different models for the high-orbit Mars probe
Statistics of navigation errors by different degree-and-order models
As can be drawn from Table 3, the navigation accuracy can be significantly improved as the higher degree-and-order nonspherical perturbation model is utilised for both the low-orbit and high-orbit Mars probes. For the 32-by-32 and 50-by-50 degree-and-order models, insufficient model accuracy leads to relatively larger errors of the navigation solution than those in the other two higher degree-and-order cases. Besides, the differences of position and velocity errors can be nearly negligible for the 70-by-70 and 80-by-80 degree-and-order models. In both cases, the position and velocity errors are respectively less than 100 m and 0⋅1 m/s for the low-orbit Mars probe. The position and velocity errors are respectively smaller than 500 m and 0⋅1 m/s for the high-orbit Mars probe.
Thus, for the purpose of real-time calculation, some tradeoff has to be made between the calculation costs and the model degree-and-order numbers. The average calculation time of 10 sets of Monte Carlo runs is 3,285⋅349 s, 6,789⋅364 s, 12,594⋅907 s and 16,366⋅194 s respectively for the 32-by-32, 50-by-50, 70-by-70 and 80-by-80 degree-and-order models. From this viewpoint, the 70-by-70 degree-and-order model is the optimal candidate. In such an optimal model case, the errors of semimajor axes are less than 1 m and 10 m respectively for the low-orbit and high-orbit Mars probes, which can be typically applied to some special Mars-probing missions. Meanwhile, the orbit inclination and ascending nodes are both less than 0⋅0015°, which can also be regarded as excellent confirmation of the precise navigation solution.
4.4 Some special simulation scenario and navigation results
In the preceding simulation scenario, two Mars probes lie in different orbital planes, namely with different orbit inclinations. The navigation results are thus obtained in a general sense. The constellation's geometrical configuration may have some effects on the navigation solution. As a special geometrical configuration, these two Mars probes are set to lie in the same orbital plane and the navigation solution is further explored. In this simulation trial, the initial states of the low-orbit Mars probe are changed but the high-orbit probe remains unchanged. For the low-orbit probe, the eccentricity is set at 0⋅001 and the inclination is set at 75°. Its semimajor axis is still set at 3,696 km. Thus the low-orbit Mars probe has the same orbit parameters except the flight height as the high-orbit Mars probe. The Monte Carlo simulation runs above are repeatedly implemented in this special navigation scenario. For a brief description, the same conclusion can still be drawn as the aforementioned navigation scenario, which is vividly illustrated from the navigation errors listed in Table 4.
Statistics of navigation errors by different degree-and-order models in the special case
It can be drawn from Tables 3 and 4 that the position and velocity errors in the two simulation scenarios are almost equivalent for the 70-by-70 degree-and-order model. But the errors of semimajor axes in the latter scenario are less than 1 m for both the low-orbit and high-orbit Mars probes. According to the statistics, the constellation configuration, including one low-orbit and one high-orbit Mars probe, has better estimation performance in the special navigation scenario. Besides, as the random errors in simulation scenarios are selected as 10 m/s and 0⋅1 m/s respectively for relative ranges and Doppler measurements, which can be easily fulfilled by the current SST measurements, the navigation results above are expected to be fully obtained even in real scenarios.
5. Conclusion and discussion
Autonomous navigation is vital for the current Mars-probing missions. This paper proposes a promising autonomous navigation manner which mainly focuses on tesseral harmonics in the Mars nonspherical perturbation. The main contribution of this paper is summarised as follows.
1. It reveals that third-body gravitation perturbation and Mars tesseral harmonics perturbation, which are added to the conventional two-body dynamic model, make autonomous navigation with SST measurements alone become weakly observable. The weak observability of absolute orbit states makes the EKF computation fail in obtaining the precise navigation solution. The UKF computation then has to be employed to accomplish the autonomous navigation because the linearised errors calculated by its EKF counterpart may exceed the tesseral harmonics perturbation.
2. Dynamic models expressed by the singularity-avoiding orbit elements are firstly proposed for the Mars probes. The selection of singularity-avoiding orbit elements can attenuate the intensive nonlinearity of dynamic models which cannot be effectively addressed by the position–velocity state expression. A reliable navigation solution can only be obtained by utilising such dynamic models. Through the tradeoff between the calculation costs and the model degree-and-order numbers, the 70-by-70 degree-and-order dynamic model is optimal for the real-time calculation onboard.
3. The effects of the geometrical configuration of the constellation on the autonomous navigation performance are investigated. When two Mars probes lie on the same orbital plane, which is a special flight scenario, a more accurate navigation solution can be obtained.
The simulation results reveal that the navigation solution in the time span of 30 days, with semimajor axis errors of less than 10 m in a general sense and less than 1 m in some special scenarios, can be effectively obtained in the case where the SST measurements alone are available. The acceptable navigation solution validates the effectiveness and availability of the newly proposed autonomous navigation manner.
The major limitation that deeply affects the navigation accuracy is the dynamic model of the nonspherical perturbation of Mars. For more recent dynamic models relying on the rapidly increasing observation data at a larger scale (Hirt et al., Reference Hirt, Claessens and Kuhn2012; Genova et al., Reference Genova, Goossens and Lemoine2016), the nonspherical perturbation of Mars may provide more accurate tesseral harmonics and thus better navigation solutions than in the GMM-2B case. Additionally, although the simulation results indicate that the nonlinearity can be overcome well by the depiction of converged navigation errors, theoretical analyses as some further consideration may still be needed to approve the weak observability of orbit states in the UKF computation framework.
