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

Abstract

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.

1. 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 magnetometer-based 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–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–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.

2. 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 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 (𝐋₁) where all three cones intersect, as illustrated in Figure 1. 𝐋₁ can be determined by solving the following: 𝐋₁⋅𝐬₁=cos𝐴₁𝐋₁⋅𝐬₂=cos𝐴₂𝐋₁⋅𝐬₃=cos𝐴₃. (1)

Figure 1 

The principle of direct calculation method.

From an observation of the Phobos, a second line of position (𝐋₂) can be determined. The position of spacecraft can be specified uniquely on the point where the two lines intersect, as shown in Figure 2.

Figure 2 

The geometry relationship of the spacecraft, Mars, and Phobos.

The geometry relationship of the spacecraft, Mars, and Phobos can be expressed as follows: 𝐫=𝐑₁+𝜌₁𝐋₁=𝐑₂+𝜌₂𝐋₂,𝜌₁𝐋₁−𝜌₂𝐋₂=𝐑₂−𝐑₁, (2)(3) where 𝐫 is the position vector of the spacecraft, 𝐑₁ and 𝐑₂ are the position vectors of Mars and Phobos, respectively, 𝜌₁ and 𝜌₂ are the distances of spacecraft to Mars and Phobos, respectively. Solving 𝜌₁ and 𝜌₂ 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.

3. 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 planet-centered 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: 𝑥̇=𝑣_𝑥+𝑤_𝑥𝑦̇=𝑣_𝑦+𝑤_𝑦𝑧̇=𝑣_𝑧+𝑤_𝑧𝑣̇_𝑥𝜇=−^𝑠𝑥𝑟³_ps−𝜇_𝑚𝑥−𝑥₁𝑟³_pm−𝑥₁𝑟³_sm−𝜇_𝑒𝑥−𝑥₂𝑟³_pe−𝑥₂𝑟³_se+𝑤_𝑣𝑥𝑣̇_𝑦𝜇=−^𝑠𝑦𝑟³_ps−𝜇_𝑚𝑦−𝑦₁𝑟³_pm−𝑦₁𝑟³_sm−𝜇_𝑒𝑦−𝑦₂𝑟³_pe−𝑦₂𝑟³_se+𝑤_𝑣𝑦𝑣̇_𝑧𝜇=−^𝑠𝑧𝑟³_ps−𝜇_𝑚𝑧−𝑧₁𝑟³_pm−𝑧₁𝑟³_sm−𝜇_𝑒𝑧−𝑧₂𝑟³_pe−𝑧₂𝑟³_se+𝑤_𝑣𝑧, (4) where (𝑥₁,𝑦₁,𝑧₁), (𝑥₂,𝑦₂,𝑧₂), and (𝑥,𝑦,𝑧) are the position of Mars, the Earth, and the spacecraft in the sun-centered inertial frame. 𝐫_se and 𝐫_sm are the position vectors of the Earth and Mars, relative to the Sun. 𝐫_ps,𝐫_pe, and 𝐫_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. 𝑊=[𝑤_𝑥,𝑤_𝑦,𝑤_𝑧,𝑤_𝑣𝑥,𝑤_𝑣𝑦,𝑤_𝑣𝑧]^𝑇 is the system noise.

When the spacecraft is in Mars sphere of influence, its orbital motion can be described as a perturbed two-body problem with Mars as the central body. The dynamical model in Mars-centered inertial frame (J2000.0) is written as follows: 𝑥̇^′=𝑣^′_𝑥+𝑤^′_𝑥𝑦̇^′=𝑣^′_𝑦+𝑤^′_𝑦𝑧̇^′=𝑣^′_𝑧+𝑤^′_𝑧𝑣̇^′_𝑥=−𝜇_𝑚𝑥^′𝑟^′3_pm−𝜇_𝑠𝑥^′−𝑥^′₁𝑟^′3_ps+𝑥^′₁𝑟^′3_ms−𝜇_𝑒𝑥^′−𝑥^′₂𝑟^′3_pe+𝑥^′₂𝑟^′3_me+𝑤^′_𝑣𝑥𝑣̇^′_𝑦=−𝜇_𝑚𝑦^′𝑟^′3_pm−𝜇_𝑠𝑦^′−𝑦^′₁𝑟^′3_ps+𝑦^′₁𝑟^′3_ms−𝜇_𝑒𝑦^′−𝑦^′₂𝑟^′3_pe+𝑦^′₂𝑟^′3_me+𝑤^′_𝑣𝑦𝑣̇^′_𝑧=−𝜇_𝑚𝑧^′𝑟^′3_pm−𝜇_𝑠𝑧^′−𝑧^′₁𝑟^′3_ps+𝑧^′₁𝑟^′3_ms−𝜇_𝑒𝑧^′−𝑧^′₂𝑟^′3_pe+𝑧^′₂𝑟^′3_me+𝑤^′_𝑣𝑧𝑟^′=𝑥^′2+𝑦^′2+𝑧^′2. (5) 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: 𝛼_𝑠−𝐫=arccos^ps⋅𝐬𝑟_ps+𝑣_𝛼𝑠𝛼_𝑒−𝐫=arccos^pe⋅𝐬𝑟_pe+𝑣_𝛼𝑒𝛼_𝑚−𝐫=arccos^pm⋅𝐬𝑟_pm+𝑣_𝛼𝑚, (6) where 𝐬 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 𝑍₁=[𝛼_𝑠,𝛼_𝑒,𝛼_𝑚]^𝑇and measurement noise 𝑉₁=[𝑣_𝛼𝑠,𝑣_𝛼𝑒,𝑣_𝛼𝑚]^𝑇, the measurement model can be presented as follows: 𝑍₁(𝑡)=𝐻₁[𝑋(𝑡),𝑡]+𝑉₁(𝑡), (7) which measurement noise covariance matrix 𝑅₁ is presented as 𝑅₁=𝐸𝑉₁(𝑡)𝑉₁(𝑡)^𝑇=diag𝜎²_𝑎𝑠,𝜎²_𝑎𝑒,𝜎²_𝑎𝑚. (8)

When the spacecraft flies far from the planets, the star-planet 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.

4. 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.

4.1. System Model

As mentioned above, the orbital motion of the spacecraft can be described by two dynamical models, which are defined by (5) and (6). They can be written in general as 𝑋̇(𝑡)=𝑓₁(𝑋(𝑡),𝑡)+𝑊₁(𝑡) and 𝑋̇(𝑡)=𝑓₂(𝑋(𝑡),𝑡)+𝑊₂(𝑡). When the spacecraft is outside Mars’ sphere of influence (𝑟_pm>5.773×10⁶), Model 1 is used. When the spacecraft is in Mars’ influence sphere (𝑟_pm≤5.773×10⁶), Model 2 is used. The system noise covariance matrices are 𝐸[𝑊₁(𝑡)𝑊₁(𝑡)^𝑇]=𝑄₁,𝐸[𝑊₂(𝑡)𝑊₂(𝑡)^𝑇]=𝑄₂, respectively.

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: 𝑍₂=𝑥_𝑛𝑘,𝑦_𝑛𝑘,𝑧_𝑛𝑘,𝛼_𝑠,𝛼_𝑒,𝛼_𝑚^𝑇, (9) 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: 𝑍₂(𝑡)=𝐻₂[𝑋(𝑡),𝑡]+𝑉₂(𝑡). (10) 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𝑖(𝑡), which measurement noise covariance matrixes 𝑅_2𝑖 are as follows:𝑅_2𝑖=𝐸𝑉_2𝑖(𝑡)𝑉_2𝑖(𝑡)^𝑇=𝑅_𝑐𝑖𝑅₁𝑅_𝑐𝑖=𝑅_𝑐𝑖(𝑖=1,2,…,𝑁) (11) where 𝑅_𝑐 is the experimental measurement noise covariance matrix of the direct calculation method. In this paper, 𝑅_𝑐 is set as diag4.6×10⁷1.9×10⁷8.4×10⁵.

4.2. 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–21], as shown in Figure 3.

Figure 3 

The IMMUKF algorithm.

4.2.1. Interaction

Suppose the Markov model transition probability matrix 𝑃 is as follows: 𝑝𝑃=¹¹𝑝₁₂⋯𝑝_1𝑁𝑝₂₁⋮𝑝₂₂⋯⋮𝑝_2𝑁⋮𝑝_𝑁1𝑝_𝑁2⋯𝑝_{𝑁𝑁𝑁𝑖}𝑝_𝑖𝑗=^𝑁_𝑗𝑝_𝑖𝑗=1, (12) 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: 𝜇_𝑖/𝑗()=1𝑘,𝑘−1𝑐_𝑗⋅𝑝_𝑖𝑗⋅𝜇_𝑖(𝑘−1), (13) where 𝜇_𝑖(𝑘−1) is the predicted model probability, 𝑐_𝑗=∑^𝑁_𝑖𝑝_𝑖𝑗⋅𝜇_𝑖(𝑘−1) is normalization constant. The mixed initial conditions for each submodel 𝑗 are as follows: 𝑋_𝑗(𝑘,𝑘−1)=^𝑁_𝑖𝑋_𝑖(𝑘−1)⋅𝜇_𝑖/𝑗()𝑃𝑘,𝑘−1_𝑗(𝑘,𝑘−1)=^𝑁_𝑖𝑃_𝑖()+𝑋𝑘−1_𝑖(𝑘−1)−𝑋_𝑗(⋅𝑋𝑘,𝑘−1)_𝑖(𝑘−1)−𝑋_𝑗(𝑘,𝑘−1)^𝑇×𝜇_𝑖/𝑗(𝑘,𝑘−1). (14)

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].

4.2.3. Model Probability Update

Consider the following: 𝑟_𝑖(𝑘)=𝑍(𝑘)−𝑍_𝑖()𝑆𝑘,𝑘−1_𝑖(𝑘)=𝐸𝑟_𝑖(𝑘)𝑟_𝑖(𝑘)^𝑇Λ_𝑖(𝑘)=𝑁𝑟_𝑖(𝑘),0,𝑆_𝑖(𝑘)=exp−(1/2)𝑟_𝑖(𝑘)^𝑇𝑆_𝑖(𝑘)𝑟_𝑖(𝑘)2𝜋𝑆_𝑖(𝑘)𝑐=^𝑁_𝑖Λ_𝑖(𝑘)𝑐_𝑖𝜇_𝑖(𝑘)=1𝑐Λ_𝑖(𝑘)^𝑁_𝑖𝑝_𝑖𝑗⋅𝜇_𝑖()=1𝑘−1𝑐Λ_𝑖(𝑘)𝑐_𝑖, (15) where 𝑟_𝑖(𝑘) and 𝑆_𝑖(𝑘) are innovation vector and its covariance, respectively. Λ_𝑖(𝑘) is Model likelihood. 𝑁(𝑥) is Gaussian distribution.

4.2.4. Combination

Combined state 𝑋(𝑘) and its covariance 𝑃(𝑘) are calculated from the weighted state estimates 𝑋_𝑖(𝑘) and covariances 𝑃_𝑖(𝑘) as follows: 𝑋(𝑘)=^𝑁_𝑖𝑋_𝑖(𝑘)⋅𝜇_𝑖(𝑘)𝑃(𝑘)=^𝑁_𝑖𝑃_𝑖(𝑘)+𝑋_𝑖(𝑘)(𝑘⋅𝑋−𝑋)_𝑖(𝑘)−𝑋(𝑘)^𝑇⋅𝜇_𝑖(𝑘). (16)(17)

5. 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.

Figure 4 

The process of simulation.

The ideal trajectory of spacecraft is created by the STK (Satellite Tool Kit), whose initial value of semimajor axis is 193,216,365.381 km, 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]: 𝑟𝑇=^mars𝑎+2𝑎1−𝑒²𝑟_marscos𝑖. (18) The Tisserand parameter 𝑇 is conserved around −0.2686 before and after the close encounter with Mars.

Figure 5 

The ideal trajectory of the spacecraft.

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 6 

The position estimation error of the direct calculation method.

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.

Table 1 

Performance comparison of different methods.

Figure 7 

Results comparison between IMMUKF and UKF.

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.

6. 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.

Acknowledgments

The research presented in this paper has been supported by the National Natural Science Foundation of China (61233005), Program for New Century Excellent Talents in University (NCET) (NCET-11-0771), and Aerospace Science and Technology Innovation Fund (10300002012117003). The authors wish to express their gratitude to all members of the Science and Technology on Inertial Laboratory and Fundamental Science on Novel Inertial Instrument & Navigation System Technology Laboratory for their valuable comments.