A New Celestial Navigation Method for Spacecraft on a Gravity Assist Trajectory

A practical and reliable capability for autonomous navigation needs to reduce operation cost, to improve operational efficiency, and to increase mission safety. Celestial navigation is a very attractive autonomous navigation solution for deep space spacecraft. There are mainly two kinds of celestial navigation methods: the direct calculation method and the filter method. The accuracy of the direct calculation method is low and very sensitive to the measurement noise.The filter method can provide a better navigation performance if a high accuracy dynamical model is available. However, the main practical problem existing in the autonomous celestial navigation of spacecraft on a gravity assist trajectory is that the accuracy of trajectory model is not enough to be used in the real navigation sometimes, which may introduce large estimation error and even cause filter divergence. To solve this problem, a new celestial navigation method is proposed in this paper, which effectively combines the direct calculation method and the filter method using an interacting multiple model unscented Kalman filter (IMMUKF). The ground experimental results demonstrate that this method can provide better navigation performance and higher reliability than the traditional direct calculation method and filter method.


Introduction
The navigation of spacecrafts is performed primarily by the ground tracking, whose navigation accuracy can reach a few kilometers by the Doppler tracking, very long baseline interferometry (VLBI), and delta very long baseline interferometry (DVLBI) of deep space network (DSN) [1,2].However, this earthbound tracking suffers from many limitations, such as huge cost, long time delay, discontinuity, and easy disturbance.Autonomous navigation has become a key requirement for reliability and safety.A variety of autonomous navigation methods for near earth satellites have been proposed and explored, including the magnetometerbased method [3,4], global position system (GPS) [5,6], intersatellite link-based method, and celestial navigation method, but for deep space explorers celestial navigation is the only feasible way.
Generally, celestial navigation methods for spacecraft can be classified into two categories: the direct calculation method and the filter method.The direct calculation method uses celestial measurements to directly calculate the position of the spacecraft according to the geometric relations between the spacecraft and celestial bodies [7].This method is very simple and is independent of the system model, but its navigation accuracy is relatively low and very sensitive to the measurement error.The filter method uses an optimal filter combined with the celestial measurements and a system model to estimate the position of the spacecraft.Because this method can deal with measurement noise, the navigation accuracy is higher than that of the direct calculation, but a highly accurate dynamical model is needed [8][9][10][11].
In deep space missions, a gravity assist trajectory is often used, which uses the gravity of a planet (or other celestial body) to alter the path and speed of a spacecraft.This technique allows to reach destinations which would not be accessible with current technology or to reach targets with significantly reduced propulsion requirements or in a reduced travel time.Many spacecrafts such as Voyager, Galileo, and Cassini use the gravity assist technique to achieve their targets.The two Voyager spacecrafts provide a classic example.Voyager 2 launched in August 1977 took one boost from Jupiter, one from Saturn, later from Uranus, and then climbed all the way to Neptune and beyond.Voyager 1 launched the following month obtained assists from Jupiter and Saturn, respectively.Galileo took one kick from Venus and two from Earth, while orbiting the Sun en route to its destination, Jupiter.Cassini passed by Venus twice, then Earth, and finally Jupiter on the way to Saturn [12][13][14].
In a gravity assist trajectory, angular momentum is transferred from the orbiting planet to a spacecraft approaching from behind the planet in its progress about the sun.The value of spacecraft's speed relative to planet is not changed during a gravity assist flyby, but its direction is changed.However, both value and direction of spacecraft's speed relative to the sun are changed during a gravity assist flyby.Planet's sun-relative orbital velocity is added to the spacecraft's velocity, and the spacecraft does not lose this component on its way out.Instead, the planet itself loses the energy.The massive planet's loss is too small to be measured, but the tiny spacecraft's gain can be very great.That is to say that the orbit of spacecraft has changed before and after a gravity assist flyby.
When the spacecraft is on such a swing-by trajectory, its orbit changes very rapidly which will make it very difficult to build an accurate dynamical model.If the inaccurate dynamical model is used in the filter method, it will introduce large estimation error and even cause the filter divergence.So, this paper presents a new celestial navigation method combining the direct calculation method and the filter method for spacecraft on a gravity assist trajectory.When an accurate dynamical model is available, the filter method is used to reduce the estimation error caused by measurement error.When such an accurate dynamical model is not available, the result of direct calculation method is used to correct the estimation error caused by inaccurate dynamical model and prevent the filter from divergence.An interacting multiple model unscented Kalman filter (IMMUKF) is used to deal with the data association problem.Compared to the traditional celestial methods, this new method can efficiently improve the accuracy and reliability in this case.The feasibility of this new method is validated by ground experiments using a spacecraft on the Mars gravity assist trajectory.
This paper is systematized in five sections.After this introduction, the principle of the direct calculation method is outlined in Section 2. Then the principle of the filter method is described, followed by the state and measurement models in Section 3. Simulation results in Section 4 demonstrate the performance improvement.Conclusions are drawn in Section 5.

Direct Calculation Method
Celestial navigation method is based on the fact that the position of celestial body (Sun, Earth, Mars, and stars) in inertial frame at certain time is known and their position observed on spacecraft in the spacecraft body frame is a function of spacecraft's position.The angle between the vector pointing toward Mars and the vector pointing toward the star from the spacecraft determines a cone.The apex of the cone is located at the position of Mars, the symmetry axis points to

Phobos Mars
Spacecraft the star, and the opening angle is twice the measured angle.The spacecraft must be on the surface of the cone.Observations of two other stars can determine two other cones.The spacecraft must be on the line (L 1 ) where all three cones intersect, as illustrated in Figure 1.L 1 can be determined by solving the following: From an observation of the Phobos, a second line of position (L 2 ) can be determined.The position of spacecraft can be specified uniquely on the point where the two lines intersect, as shown in Figure 2.
The geometry relationship of the spacecraft, Mars, and Phobos can be expressed as follows: where r is the position vector of the spacecraft, R 1 and R 2 are the position vectors of Mars and Phobos, respectively,  1 and  2 are the distances of spacecraft to Mars and Phobos, respectively.Solving  1 and  2 from (3), the position of spacecraft can be determined by (2).The limitation of the direct calculation method is that the measurement error has a great impact on the calculation error and it can not provide the velocity information.

Filter Method
The filter method exploits the dynamics of the spacecraft to remove the effects of the noise and obtain a good estimation of the position and velocity of the spacecraft at the present time using the differences between the measured and estimated celestial data.In the filter method, the state (dynamic) and measurement models are needed.The dynamical model describes how the state variables change over time, and the measurement model describes the relationship between the measurement variables and the state variables.
The patched-conic technique [15,16] is a useful tool for designing a fuel-saving gravity assist trajectory.In this technique, a two-body orbit about the Sun is first kinematically patched onto a two-body hyperbolic orbit about the planet.Then the spacecraft proceeds along the planetcentered hyperbolic orbit during the flyby, and finally a new two-body orbit about the Sun is patched onto the postflyby hyperbolic orbit.For a spacecraft on Mars gravity assist trajectory, its trajectory can be broadly divided into two parts.
When the spacecraft is outside Mars, sphere of influence, its orbital motion can be described as a multibody problem with the sun as the central body.Its dynamical model is expressed in the sun-centered inertial frame (J2000.0)as follows: where ( 1 ,  1 ,  1 ), ( 2 ,  2 ,  2 ), and (, , ) are the position of Mars, the Earth, and the spacecraft in the sun-centered inertial frame.r se and r sm are the position vectors of the Earth and Mars, relative to the Sun.r ps , r pe , and r pm are the position vectors of the spacecraft relative to the Sun, the Earth and Mars, respectively.  ,   , and   are the gravitational constants of the Sun, the Earth, and Mars respectively. = [  ,   ,   ,  V  ,  V  ,  V  ]  is the system noise.
When the spacecraft is in Mars sphere of influence, its orbital motion can be described as a perturbed twobody problem with Mars as the central body.The dynamical model in Mars-centered inertial frame (J2000.0) is written as follows: All parameters used in (5) have the same meaning as those in (4), except that these parameters are expressed in Mars-centered inertial frame instead of sun-centered inertial frame.
Because these two models are in different frames, the final results have to be transformed to the sun-centered inertial frame.
The usually used celestial measurements include (1) the directions of near celestial bodies (the sun and the planet), (2) the planet-planet angles, and (3) the star-planet angles.Here we take the last one as an example and the corresponding measurement equation is given by the following: where s is the position vector of the navigation star in the sun-centered inertial frame, which is obtained from the star catalog by star identification.Assume measurement , the measurement model can be presented as follows: which measurement noise covariance matrix  1 is presented as When the spacecraft flies far from the planets, the starplanet angles change very slowly.The star-planet angles alone cannot supply enough information for position estimation without the help of an accurate dynamical model in the filter method.However, when the spacecraft is on a gravity assist trajectory, its orbit changes very fast and an accurate dynamical model is not available, which may introduce large estimation error and even cause the filter divergence.

New Method Based on IMMUKF
For solving problems of above two methods, a new celestial navigation method is presented for the flyby spacecraft.In this new method, the direct calculation method is used to get a preliminary position estimation result, and then this preliminary result is used as a supplementary measurement in the filter.The difficulty existing in the new method is that the measurement noise covariance matrix is not constant because the position estimation error of the direct calculation method changes over time due to the movement of spacecraft.An IMMUKF method is used to deal with this problem and to enhance the adaptive capability of the filter.

Model probability update 𝜇(𝑘)
Combination In this new method, the position estimation result obtained by the direct calculation method is used as a supplementary measurement besides the star-planet angle measurement defined by (7).The measurement used in this new method can be described as follows: where [  ,   ,   ]  is the position estimation result obtained by the direct calculation method. is a scale factor, which is used to make the magnitude of position measurement match that of star-planet angle measurement.The measurement model can be presented as follows: Because the measurement noise covariance matrix of this model is unknown and time varying, a set of different matrices is used to approximate it.Then  measurement submodels are presented as  2 () =  2 [(), ] +  2 (), which measurement noise covariance matrixes  2 are as follows: where   is the experimental measurement noise covariance matrix of the direct calculation method.In this paper,   is set as diag {4.6 × 10 7 1.9 × 10 7 8.4 × 10 5 }.

IMMUKF Method.
The interacting multiple model (IMM) algorithm is a very effective estimation algorithm for systems which can be described by a mixing of continuous states and discrete modes [17,18].The switching between submodels is governed by a Markov chain.The IMMUKF method includes interaction, filtering, model probability update, and combination of four steps [19][20][21], as shown in Figure 3. . . .
where   is the a priori probability of transition from model  to model .The initial value of  in this paper is set as   = 0.91 with  =  and   = 0.01 with  ̸ = .
The mixing probabilities at time  − 1 are computed as follows: where   ( − 1) is the predicted model probability,   = ∑     ⋅   ( − 1) is normalization constant.The mixed initial conditions for each submodel  are as follows: 4.2.2.Filtering.Because the state model and measurement model described above are nonlinear, the optimal state and covariance estimates   () and   () for each sub-model  are calculated using UKF [22,23].

Model Probability Update.
Consider the following: where   () and   () are innovation vector and its covariance, respectively.Λ  () is Model likelihood.() is Gaussian distribution.

Combination.
Combined state () and its covariance () are calculated from the weighted state estimates   () and covariances   () as follows:

Simulation Results
This section presents simulation results to show the performance improvement of this new method compared with the traditional direct calculation method and the filter method.The process of simulation is shown in Figure 4.
The ideal trajectory of spacecraft is created by the STK (Satellite Tool Kit), whose initial value of semimajor axis is 193,216,365.381km, eccentricity 0.236386, inclination 23.455 ∘ , right Ascension of the ascending node 0.258 ∘ , and the argument of periapsis 71.347 ∘ .The near bodies used in the direct calculation method are Mars and Phobos.All measurements are created by the celestial navigation simulation platform [24] using ideal trajectory of STK, ephemeris of DE405, and Tycho-2 catalogue.The ideal trajectory of the spacecraft is shown in Figure 5, which has been checked according to Tisserand's criterion described in (17) [18,20]: The Tisserand parameter  is conserved around −0.2686 before and after the close encounter with Mars.The accuracy of the Sun sensor, Phobos' sensor, the Earth sensor, and Mars sensor are selected to be 5  , 5  , 1.5  , and 0.005 ∘ , respectively.The accuracy of the star sensor is 3  .
Figure 6 shows the performance of the direct calculation method, which is obtained with a 1-minute sampling interval during the 20-day period when the spacecraft is near Mars.As it can be seen from Figure 6, the navigation accuracy of this method mainly depends on the accuracy of measurement, so it is very sensitive to the measurement noise.The geometry relationship of the spacecraft, Mars, and Phobos also has an impact on the navigation performance.The shorter the distance between spacecraft and Mars (or Phobos), the higher the navigation accuracy is.
Figure 7 shows the performance comparison between this new method and the traditional filter method, which only uses the star-planet angles as measurement, under the same simulation conditions.As demonstrated in Figure 7, the traditional filter method can provide a good navigation performance when the spacecraft is far from Mars and a highly accurate model is available, but its estimation error increases greatly when the spacecraft is close to Mars.This is because the star-planet angles alone cannot supply enough   information to correct the large estimation error introduced by the inaccurate dynamical model, which even cause the divergence of the filter after then.The innovation sequence of the two methods which has subtracted the measurement noise is shown in Figure 7.The figure of the innovation sequence also demonstrates the divergence of the traditional filter method.However, the new method presented in this paper shows a stable and satisfying navigation performance.The details are given in Table 1.
From these results, we see that this IMMUKF-based new celestial navigation method overcomes the shortcomings of the direct calculation method and the traditional filter method.On one hand, when an accurate dynamical model is available, the impact of measurement error on the navigation performance is reduced by the dynamical model.On the other hand, when an accurate dynamical model is not available, the impact of system error on the navigation performance is reduced by the result of direct calculation method.The navigation accuracy and reliability are greatly improved by the effective combination of the direct calculation method and the filter method.These results demonstrate that this new method is a particularly suitable celestial navigation method of spacecrafts on a gravity assist trajectory.

Conclusions
A new autonomous celestial navigation method for spacecrafts on a gravity assist trajectory is studied, which effectively combines the direct calculation method and the filter method based on IMMUKF.The navigation performance of this new method is tested by the ground experiments, and a position estimation error within 20 km is obtained.Compared with the traditional methods, this method can achieve a higher accuracy and stability.These results verify that this method is a promising and attractive scheme of autonomous navigation for the spacecraft using gravity assist.A limitation of this study is the significant increase in computational cost.This method should have a broad application potential in autonomous navigation of deep space probes, especially in cases when the accurate dynamical model is not available.How to enhance the adaptive capability and to reduce the computational cost is under further investigation.

o n e o f p o s i t i o n 1 The cone of position 2 Th e c o n e o f p o s it io n 3 1 Figure 1 :
Figure 1: The principle of direct calculation method.

Figure 2 :
Figure 2: The geometry relationship of the spacecraft, Mars, and Phobos.

Figure 4 :
Figure 4: The process of simulation.

Figure 5 :
Figure 5: The ideal trajectory of the spacecraft.

Figure 6 :
Figure 6: The position estimation error of the direct calculation method.
2.1.Interaction.Suppose the Markov model transition probability matrix  is as follows:

Table 1 :
Performance comparison of different methods.