Navigation performance analysis of Earth– Moon spacecraft using GNSS, INS, and star tracker

Global Navigation Satellite System (GNSS) can provide an approach for spacecraft autonomous navigation in earth– moon space to make up for the insufficiency of earth-based tracking, telemetry, and control systems. However, its weak power and poor observation geometry near the moon causes new problems. After the GNSS signal characteristics and satellite visibility were evaluated in Phasing Orbit and Lunar Transfer Orbit, we proposed an adaptive Kalman filter based on the Carrier-to-Noise ratio ( C / N 0 ) and innovation vector to weaken the influence of GNSS accuracy attenuation as much as possible. The experimental results show that the spacecraft position and velocity accuracy are better than 10 m and 0.1 m/s near the Earth, and better than 50 m and approximately 0.2 m/s near the moon use GNSS with the proposed adaptive algorithms. Additionally, because of the deterioration of navigation performance based on the orbit filter during orbital maneuvering, we used accelerometer data to compensate for the dynamic model to maintain navigation performance. The results of the experiment provide a reference for subsequent studies.


Introduction
Space research institutions worldwide have been investing significant energy into the moon, with numerous exploration missions planned to the moon as the destination.The 2018 and 2022 Global Exploration Roadmap supplement described the latest architecture for lunar surface missions, and contained a map of human mission planning for lunar exploration over the next decade (International Space Exploration Coordination Group, 2018, 2022).In recent years, notable lunar exploration missions have included the Chinese Chang'e-V, American Artemis, and Japanese HAKUTO-R.Lunar exploration has emerged as the initial stride toward future human deep space exploration missions and will captivate the spotlight in the upcoming decade.
It is well known that the determination of the spacecraft position, velocity, and attitude is an important part of the exploration mission, and is currently the responsibility of the ground tracking, telemetry, and control group.The group forms a deep space network by building several monitoring stations on the ground, and uses the unified S-band system and very long baseline interferometry to determine the position and velocity of spacecraft (Erhu et al., 2023;Liu et al., 2022;Ulvestad, 1992).Although high-precision position and velocity can be obtained, the limitations of ground and on-board computing resources make it difficult to achieve the realtime tracking, telemetry, and control of large numbers of spacecraft.Therefore, spacecraft real-time autonomous navigation has been a research direction with much attention (Ely et al., 2022).
Since the construction of the Global Navigation Satellite System (GNSS), many scholars have come to understand that it not only serves terrestrial users and Low Earth Orbit (LEO) satellites by using main lobe signals but also extends its service to Medium Earth Orbit (MEO), High Earth Orbit (HEO), and even ultra-HEO through its side lobe signals.The full Space Service Volume (SSV) of GNSS is shown in Fig. 1.The region from the Earth's surface to 3000 km altitude is usually called the Terrestrial Service Volume (TSV), the area from 3000 to 36,000 km is named SSV, and the area above 36,000 km is called Extended SSV (ESSV) (Bauer et al., 2017).In 2012, Work Group B (WG-B) of the International Committee on GNSS (ICG) identified the advantages of an interoperable GNSS SSV for the space user community.Subsequently, they have been working with GNSS service providers and space agencies to advance the SSV and consider its applications within and beyond cis-lunar space.Then WG-B officially released two booklets titled "The Interoperable Global Navigation Satellite Systems Space Service Volume" in 2018 and 2021, and the SSV characteristics of each constellation and service performance in space missions were analyzed (United Nations, 2021).
Other agencies and scholars have also focused on the performance of GNSS services in earth-moon space, and the characteristics and navigation performance of Global Position System (GPS) and Galileo satellite navigation system (Galileo) signals have been analyzed more deeply.The European Space Agency is interested in lunar exploration and motivated an experimental platform to evaluate the navigation performance of GPS L1C/A, L5Q, Galileo E1C, and E5aQ signals in Lunar Transfer Orbit (LTO), Low Lunar Orbit (LLO), and descent and landing (D&L) in detail (Lopes et al., 2014).One of the most important points is satellite visibility, which needs to be modeled and simulated as accurately as possible.Shehaj et al. (2017) found that the Three-Dimensional (3D) GNSS antenna pattern model obviously affects satellite visibility in ultra-high orbit compared with the two-dimensional model.Accurate GPS and Galileo 3D antenna patterns were used to analyze the Carrier-to-Noise ratio (C/N 0 ), satellite availability, and Position Dilution of Precision (PDOP) of receivers in various orbits around the Moon.Delépaut et al. took into account the navigation message's C/N 0 demodulation threshold, and an "ephemeris-based visibility" approach was proposed to determine signal visibility (Delépaut et al., 2019(Delépaut et al., , 2020)).Also, some studies used GPS signals and inter-satellite links to determine the orbit and the time synchronization of the Lagrange orbit near the Moon, which verified the availability of GNSS in the global space of the Earth and moon (Iiyama et al., 2021;Qi & Oguri, 2023;Sirbu & Leonardi, 2023).As for BeiDou Navigation Satellite System (BDS), Lin et al. (2019) analyzed the propagation characteristics of the SSV signal in the ionosphere.Then Lin et al. (2020) used the 3D antenna pattern of BeiDou-3 Navigation Satellite System (BDS-3) satellites to simulate and analyze the Doppler shift, C/N 0 , and measurement noise of the navigation signal for various space missions.The SSV service performance of BDS-3 and the advantages of Inclined Geo-Synchronous Orbit (IGSO) and Geostationary Orbit (GEO) satellites in SSV service were analyzed (Guan et al., 2022;Ma et al., 2023).
However, the performance of GNSS decreases with the distance from the Earth, so that it cannot provide high-precision navigation services for users in lunar space.Therefore, scholars have explored the integration of multi-sensors for more accurate spacecraft estimation state.Capuano et al. (2014) has been attempting to implement integrated navigation in ultra-high orbit since 2014, and proposed the integration of the Inertial Navigation System (INS), GNSS, star tracker, and orbital filter to obtain higher precision and robust navigation performance.An autonomous navigation system was designed and used to provide real-time autonomous navigation and attitude determination for LEO, GEO, and even higher HEO missions (Capuano et al., 2017(Capuano et al., , 2018)).However, accelerometers can only measure non-gravitational forces in space and cannot provide a complete dynamic model for spacecraft (Wei et al., 2023), some conclusions may need to be refined more clearly.
As a common attitude determination sensor, star tracker can provide a high-precision attitude for the vehicle.Therefore, some multi-sensors integration methods Fig. 1 Signal reception geometry in the space service volume and adaptive algorithms have been proposed, including decentralized data fusion methods such as federated filtering, in addition to centralized extended Kalman filtering (Gao et al., 1993).At present, it widely integrated with INS and GNSS to provide autonomous navigation services for ground and near-earth vehicles, even include lunar and Mars surface probes (Gao et al., 2018;Hu et al., 2023;Ning et al., 2011).Different from the general scenario, star tracker is more commonly integrated with an Earth sensor and orbital filter to provide the position and velocity for spacecraft (called starlight angle navigation), in addition to providing the spacecraft with a high-precision attitude, and is widely used in deep space probes because of its autonomy (Bhaskaran et al., 2000;Yu et al., 2021).However, it is limited by the performance of the Earth sensor, which results in low navigation performance.Many scholars have attempted to improve navigation accuracy by combining laser ranging, accelerometer, and other methods (Gui et al., 2021;Li et al., 2022a;Wang & Zhang, 2013).With the improvement of the accuracy of the earth sensor, the integration of starlight angle navigation and GNSS will be the main autonomous navigation method for the future Earth-moon space spacecraft and requires further study.
After studying various difficulties in the autonomous navigation of earth-moon spacecraft, we believe that the primary problem for GNSS is to overcome precision attenuation as the distance increases.Additionally, the spacecraft will suffer from inevitable dynamic disturbance, such as orbital maneuver, descent, and landing, which will also cause a significant decline in navigation performance.
To solve these problems, we first simulate the observations of the GNSS receiver, IMU, and star tracker based on the Phasing Orbit (PHO) and LTO of the lunar probe, and analyze the signal characteristics and visibility of GPS and BDS in earth-moon space.Then we propose an adaptive algorithm based on the innovation vector to reflect the change of observation noise more accurately to suppress the precision attenuation of GNSS near the moon and improve spacecraft navigation performance.Finally, we use accelerometers during orbital maneuvers to compensate for the dynamic model, and integrate INS, GNSS, and star tracker to maintain navigation performance.
This paper is structured as follows: the models and algorithms were described for the adaptive Kalman filter in detail in "Models and algorithms" section.Then the modeling and simulation approach for the sensors were introduced in "Observation simulation and characteristic analysis" section.Finally, we analyze the solution results of various schemes used to evaluate the proposed approaches in "Experiments and analysis" section.And the conclusions and discussions are drawn in "Conclusions and discussion" section.

Models and algorithms
First, to respond to the various states of earth-moon spacecraft, an autonomous navigation method based on orbital dynamics, assisted by GNSS, INS, and star tracker, is shown Fig. 2. Given the clear understanding of orbital dynamics, we use GNSS based on orbital filtering as the primary navigation method.Considering the onboard accelerometer can only measure non-gravitational acceleration so that the non-gravitational model can be replaced by accelerometer data.We only use the accelerometer observations when the spacecraft is subjected to additional and unmodeled non-gravitational accelerations acceleration (e.g., the thrust of the thruster).Regarding the star tracker, we only use it as an attitude determination system in our experiments.

Setting of the Kalman filter
For nonlinear stochastic model, the dynamic model and measurement model are shown in Eq. ( 1) and (2) as follow: where For conventional orbital filtering, the unknown state parameters X k of the k-th epoch can be expressed as which comprises the 3D position and velocity vectors r and v in Earth-centered inertial reference frame (i-frame); one solar radiation pressure coefficient C R ; the receiver clock offset parament dt in seconds (Li et al., 2022b); and the speed of light c.
The state transfer matrix can be calculated using the orbit dynamic model: (1) is the state transition matrixes of orbital dynamics; t is the time interval and I is a unit matrix.
As mentioned above, on-board accelerometers have sufficient accuracy but only can measure non-gravitational acceleration.The forces on the spacecraft during orbital maneuvers can be expressed as The gravitational force F G include the Earth's gravity F earth , Sun's gravity F sun , moon's gravity F moon , and grav- ity of other celestial bodies F n .The non-gravitational force F NG mainly includes solar radiation pressure F srp , thruster thrust F man , and other unmodeled forces F ε .The maneuvering time is usually short and the gravitational force is relatively stable, whereas the non-gravitational force can be measured by the accelerometer.Therefore, the pseudo-observations of IMU can be constructed and the spacecraft state during orbital maneuvering can be estimated using the integrated navigation system.
In general, the on-board accelerometer is installed parallel to the spacecraft body frame (b-frame), and the non-gravitational force F NG under the i-frame can be expressed as where C i b is the rotation matrix from b-frame to i-frame, which can be obtained from the spacecraft attitude.The dynamic equation of INS in the i-frame can be described as where r and v denote position and velocity vectors in the i-frame, respectively; f b ib is the specific force vector in the b-frame; F G is the gravitational forces in the i-frame which is related to the position r and can be obtained by precise dynamic model.Ω b ib is the skew-symmetric matrix of ω b ib , while ω b ib is the rotation angular rate of gyroscopes output.
Usually, in order to suppress the error dispersion of INS, in addition to integrating with GNSS, it is also integrated with star tracker to form an INS, GNSS and Star tracker integrate navigation system.And the accurate position, velocity and attitude can also be obtained even if the dynamic model provided by orbital mechanics is ( 8) The parameters to be estimated for the integration system can be expressed as where, in addition to 3D position vector r , velocity vector v , receiver clock offset dt ; additional parameters include spacecraft 3D attitude ψ , the accelerometer bias b a and gyros bias b g .The state transition matrix Φ INS k,k−1 of the INS in the i-frame can be expressed as where N g = ∂ � F g /∂� r , ∧ is the symbol for calculate the skew-symmetric matrix, I 3 is a 3-order unit matrix and 0 3 is a 3-order zero matrix.
After we get the state transfer matrix for orbital filtering Φ k,k−1 and that for INS Φ INS k,k−1 , respectively, the predicted state Xk , covariance matrix Pk , and epoch time k can be expressed as where X k−1 and P k−1 are the estimated state and corre- sponding covariance matrix, respectively.Q k is the covar- iance matrix of system noise and represents the accuracy of the dynamic model.
As for the orbit filter, the pseudo-range ρ of GNSS as the observation, if n satellites are observed at epoch time k, the measurement matrix H k , and measurement inno- vation δz k can be expressed as and ( 11) where the predicted value and measured value of * are represented by * and * , respectively.e i is the unit vector of the line-of-sight direction from the GNSS receiver to the i-th satellite.z k is the observation at epoch time k.
As for INS, GNSS and Star tracker integrate navigation system, if m observations of the Star tracker also be considered, the measurement matrix and measurement innovation can be rewritten as and where l j is the unit vector of the j-th star in the b-frame; and r S j is the unit vector of the j-th star in the i-frame.
In addition, according to Eq. ( 10), the attitude dynamics does not depend on other parameters except rotation angular rate, so the attitude of the spacecraft can be independently calculated according to the outputs of the gyroscope and the star tracker.
Then the measurement is updated, and gain matrix K k , estimated state X k , and covariance matrix P k can be calcu- lated using where gain matrix K k determines the weight of obser- vation information when the state is updated.R k is the covariance matrix of observation noise and is a diagonal matrix composed of independent observation variance σ 2 at epoch time k.If observation noise increases, it is necessary to reduce the weight of observation information by adjusting the R matrix and reducing matrix K k .( 17) e 1 0 3 0 3 0 1 . . . . . . . . . . . . . . .

GNSS fault detection and exclusion base on integrated system
As we all know, GNSS observations will inevitably have failures in practical situation, especially in the environment of weak signals in deep space, will make the wrong state parameters output by the estimator and reduce the accuracy of autonomous navigation.Therefore, Fault Detection and Exclusion (FDE) must be applied before the integrated solution.
In the Kalman filtering of integrated system, the predicted state X k and measurement vector z k are inte- grated.The unknown parameters are estimated through the least squares estimation, and the measurement model is given by (Aghapour & Farrell, 2019;Hewitson & Wang, 2007;Jiang et al., 2020): where A k denotes the design matrix, V k denotes the residual vector, and where V z k denotes the residual vector of the measure- ment, V ỸK denotes the residual vector of the predicted state, C L k denotes the covariance matrix corresponding to the measurement model (Hewitson & Wang, 2007).Based on the least squares estimation, the estimation of state parameters Ỹ k and its covariance matrix Q XK are obtained by Residual vector of the measurements V k and its covari- ance matrix Q V k are derived by: ( 22) In the integrated navigation system, the test statistics k is defined based on the residual and its covariance matrix, namely, k obeys the Chi-square distribution with the freedom degree of m-n, namely, where m is number of observations and n is number of unknown state parameters.It indicates no failures if the statistics satisfy the Eq. ( 32), once k is bigger than χ 2 m−n (α) , the fault is detected.
After the failure identification is applied to locate failure and a fault has been detected with the global detection algorithm, the ω-test can then be used to identify the corresponding measurement (Teunissen, 1998), where the test statistic is defined as where, T is a zero vector whose ith element is set to 1. ω i obeys the normal distribution when no failures exist, namely, ω i ∼ N (0, 1) .In the pres- ence of outliers, the test statistics would be bigger than the threshold T D = N α/2 (0, 1) ( α here is set to 0.1%, and the corresponding threshold T D = 3.2905 ), and the i-th measurement is thought failure.Accordingly, the failure identification is completed.

Adaptive estimation algorithm
As mentioned above, the real noise v of the observa- tion is difficult to obtain; hence, matrix R is generally set to an empirical constant.However, different from the ground environment, when a spacecraft flies away from the Earth, the measurement noise of the GNSS signal increases significantly because of the free space propagation loss and the positioning geometry becomes worse inevitably.Therefore, noise covariance matrix R needs to be adjusted adaptively to the approximate actual measurement noise.In this section, we introduce two adaptive estimation algorithms: "Measurement Noise-Based Adaptive Estimation" (MBAE) and "Innovation-Based Adaptive Estimation" (IBAE).

MBAE
The measurement noise is mainly caused by sign free space loss in deep space.So that in addition to setting the initial measurement noise, we also use a formula related to C/N 0 in every epoch to make the noise covariance matrix approach the actual measurement noise more accurately.The covariance matrix of the measurements can be given by where σ ρ i k is the pseudo-range measurement noise of the i-th satellite, and is set to 10 m.The coefficient m can be calculated using a modified sigma-Δ model and calculated from the receiver C/N 0 (Brunner et al., 1999), which can be described as where the coefficients a and b are obtained using a priori fitting based on different navigation systems and different frequencies.Using the above equations, the adaptive estimation of the measurement noise matrix at low C/N 0 can be achieved, which results in an approximation of the actual measurement noise.

IBAE
Based on the MBAE algorithm above, we can obtain the covariance matrix of the theoretical measurement noise.However, obviously, the actual measurement noise is caused by various factors, such as multipath effects, in addition to signal free space loss, Eq. ( 35) cannot fully represent the real measurement noise.Therefore, additional algorithms are needed to further improve the navigation accuracy with weak signals.
When the dynamic model is accurate and observation noise increases, the measurement innovation δz k obviously increases, which is the difference between the observation value z k and predicted value H k Xk .Accord- ing to this characteristic, we construct an adaptive algorithm based on measuring innovation.
According to Eq. ( 3) and ( 5), Eq. ( 3) can be rewritten as In addition, covariance matrix Pk in the Kalman filter is the mathematical covariance of the difference between predicted state Xk and truth value X k , which can be expressed as Therefore, the squared uncertainty of the vector H k Xk predicted by the predicted state Xk is Consequently, according to Eq. ( 16), The mathematical expectation and covariance of the measurement innovation δz k can be expressed as So far, we can find that the measurement innovation δz k obeys a Gaussian normal distribution with a mean of zero and a variance of Generally, in the process of state estimation using Kalman filter, if the variance R k+1 of noise w k+1 increases for the measurement z k+1 at epoch time k + 1 and the dynamic model is accurate and stable, the variance H k+1 Pk+1 H T k+1 + R k+1 of the measurement innovation δz k will also increase.In this case, it is necessary to use the variance expansion method to reduce the weight of observation information in the state update.
According to Eq. ( 36) and (37), a standardized innovation vector S δz,k at epoch k can be defined as follows where the standardized innovation vector Sδz,k obeys the normal distribution when the covariance R k can be given accurately.It is obviously impossible and is set to an empirical constant in the case of high quality and stable observations.But it cannot be applied to deep space environments with weak signals.
Therefore, the sigma-Δ model is used to obtained the prior covariance matrix R MBAE k , which is more accuracy than the empirical constant.Then the covariance matrix R k in Eq. ( 42) can be rewritten as R MBAE k .We can see that the standardized innovation vector can reflect the change of observation noise.Therefore, adaptive factor ω k can be constructed and the observation noise covariance matrix R MBAE k of the k-th epoch can be reconstructed as Rk .The calculation method for adaptive factor ω k and the reconstruction method for the meas- urement noise covariance matrix, respectively, are given as follows: (39) where k 0 and k 1 are two empirical constants, which are usually selected in the ranges of 2.0-3.0 and 4.5-8.5, respectively (Yang et al., 2010).In this paper, we set the values to 2 and 7.5, respectively.Equation ( 42) shows that when innovation vector δz k is small, the quality of obser- vations is relatively stable; and an increase in vector δz k means that the noise of the observation increases, which results in the increase of the vector Sδz,k .At this time the adaptive factor ω k need to be used to expand covariance matrix R MBAE k , and reduce the weight of observation in the state update.However, when observation noise is too large, the weight of observation needs to be minimized.The covariance matrix Rk obtained by the above method is substituted into Eq.( 19) for subsequent measurement updates to obtain estimated parameters.( 43) The IBAE algorithm can adaptively estimate the measurement noise covariance matrix according to the measurement innovation, which is suitable for earth-moon space, where GNSS measurement noise changes greatly; however, the initial value must be as accurate as possible.Hence, it is necessary to use the MBAE algorithm to determine the noise covariance matrix before the IBAE algorithm is applied.The adaptive Kalman filter solving process with the FDE, IBAE and MBAE algorithms is shown in Fig. 3.
We used simulation data to verify the effectiveness of above algorithms.In "Observation simulation and characteristic analysis" section, we focus on the simulation of observed data from multiple sensors and the GNSS signal characteristics was also analyzed.The navigation performance of earth-moon spacecraft was analyzed with the simulated observations and the adaptive algorithms in "Experiments and analysis" section.

Parameters of the simulation orbit
We consider a typical earth-moon space exploration project: "Chang'e-I" (named as CE-1) as an example.As shown in Fig. 4, the orbit of the CE-1 spacecraft consists Fig. 3 Adaptive Kalman filtering algorithm procedure of PHO and LTO.PHO includes multiple HEOs with 16-h, 24-h, and 48-h periods.The spacecraft makes several adjustments in PHO before entering LTO and finally reaching the Moon.The 48-h PHO and LTO for the data simulation and analysis were selected and the detailed initial state parameters and dynamic model settings used in the simulation are shown in Tables 1 and 2, respectively.
We simulated the spacecraft attitude using "Satellite Tool Kit" software and set the attitude mode to "Nadir alignment with ECI velocity constraint." The body Z axis points towards the nadir, and the X axis is constrained by the direction of inertial velocity.Then the observations from multi-sensors can be simulated.

GNSS
Celestial eclipse and signal propagation loss are the main factors that affect satellite visibility and the signal measurement accuracy of deep space GNSS receivers.
Obviously, the main eclipse body for PHO and LTO in earth-moon space is the Earth, and for LLO and D&L, the Moon is also included.
We use the geometric approach to determine whether the satellite signal is blocked by the Earth.As shown in Fig. 5, θ main is the threshold of the available range of the  GNSS main lobe signal, beyond which the side lobe signal is considered to be received.is the off-boresight angle of the receiver antenna, which is related to the receiver antenna gain.Additionally, eclipse angle θ Earth is also defined, and the spacecraft-satellite-Earth angle is marked as θ .According to geometry, satellite visibil- ity and the type of signal received can be determined.
The earth eclipse angle and main lobe angle of BDS-3 and GPS are shown in Table 3 (Rathinam & Dempster, 2016).Signal propagation loss and receiver receiving power are key parameters in the simulation.We choose BDS-3 B1C and GPS L1 as the frequency of simulated observations.The simulation parameter settings of the observations are shown in Table 4.
The simulation methods to calculate the C/N 0 of the GNSS receiver is shown below: where P T is the output power of the transmitter ampli- fier; G T is the transmitter antenna gain, which relates to the elevation and azimuth; G R is the receiver antenna gain, which relates to the off-boresight angle; and L S rep- resents the free space propagation loss.L A is the influ- ence of the atmosphere and we set it to 0 dB Hz in our simulation because of the small influence.N 0 is the noise power spectrum density and can be calculated as where k = 1.3806505 × 10 −23 J/K is the Boltzmann con- stant and T s is the system noise temperature.
After we obtain the receiver C/N 0 using Eq. ( 45), the theoretical code noise can be calculated according to the signal modulation method.For GPS L1 signals using Binary Phase Shift Keying (BPSK) modulation, thermal noise code tracking jitter σ DLL can be calculated as where B n = 0.05 is the code loop noise bandwidth in Hz, D = 0.25 is the early-to-late correlator spacing in chips.T = 0.02 is the coherent integration time in seconds.
As we described in the introduction, an accurate 3D antenna model is important for simulations.Therefore, we used the BDS-3 and GPS 3D antenna model to achieve a more accurate simulation.Figures 6 and  7 show the antenna equivalent isotropically radiated power, which includes the power amplifier and transmitting antenna gain.GPS L1 models of all types of satellites were released by Lockheed Martin (Marquis & Reigh, 2015) and the BDS-3 B1C models can be obtained from Lin et al. (2020).
Then we assumed that the spacecraft is equipped with nadir and zenith antennas that can capture GPS and BDS-3 signals.A high-gain antenna is mounted on the nadir receiver for the weak signals.The receiver antenna gains of the two antennas are shown in Fig. 8.The acquisition and tracking threshold of the receiver adopts the recommended value, which is set to 15 dB Hz (United Nations, 2021).Now, we can analyze the signal characteristics of GPS and BDS navigation signals received by spacecraft in PHO and LTO.Additionally, we can refer to the calculation method for the Doppler shift (Lin et al., 2020).The two subplots in Fig. 9 show the Doppler shifts (top), C/N 0 (middle), and number of visible satellites (bottom) for PHO and LTO, respectively.
We can determine that, at the beginning of PHO, the main lobe signal Doppler shift variation range can reach ± 50 kHz and is then below ± 15 kHz after 3 h, and the side lobe signal Doppler shift is approximately ± 25 kHz in the entire orbit.Regarding LTO, the main lobe signal has a lower Doppler shift variation range, which is approximately ± 13 kHz, and the side lobe signal is approximately ± 25 kHz.Additionally, we can see that the main lobe received signal C/N 0 is According to Eq. ( 45), the minimum, maximum, and average values of the thermal noise code tracking jitter can be obtained for PHO and LTO, as shown in Fig. 10.We can see that the disparity of the ranging noise of different satellites is large, the maximum value is more than 10 m, and the minimum value is better than 5 m, which is mainly related to the satellite antenna transmitting power and free space loss.

IMU and star tracker
An aerospace IMU was be used and the parameters of the accelerometer and gyroscope are shown in Tables 5 and  6, respectively (Capuano et al., 2014;Groves, 2013).
Regarding the star tracker, the diagonal field size of the simulated star tracker is set to 10 degrees, the focal length is 85 mm, and the threshold of the stellar magnitude is 7 Mv.Additionally, the Hipparcos catalog is chosen as the navigation catalog (ESA, 1997).Considering the accuracy of the star tracker is mainly affected by the recognition error of starlight in the image coordinate system, we set the recognition error to 1 × 10 −5 degrees to be as close as possible to the actual error.Figure 11 shows the attitude solution error of the simulated star tracker with a time interval of 1 s.The accuracy of the pitch, roll, and yaw are 3.87, 31.43, and 5.31″, respectively, which is consistent with "Blue Canyon Technologies Nano-Star Tracker" mentioned in Capuano et al. (2014).Additionally, it is worth noting that the accuracy of the roll angle is much lower than that of the pitch and yaw, which is also reasonable (Liebe, 2002).

Experiments and analysis
After we obtained the simulate observations of multisensors, the accuracy of orbit determination using multisensors can be analyzed this section.

FDE for GNSS observations
Firstly, in order to maintain the reliability and stability of state solution in low C/N 0 environment, the gross errors

Parameter Value
Accelerometer biases in x, y, z direction ( G 1 = 9.80665 Accelerometer scale factor and cross coupling errors   of GNSS observations must be detected and excluded before state estimation.
The LTO orbit at the 48th hour was selected as the target orbit and the gross errors was added to the observations from three GNSS satellites, with the detailed configuration shown in Table 7.
Then FDE algorithm mentioned above is used to detect and eliminate gross errors, Fig. 12 provides the value of k and Chi-square distribution threshold χ 2 m−n (α = 0.01) in the calculation process.
It can be seen that gross errors of C24 and C25 will cause value of k to increase significantly, which exceeds the corresponding detection threshold and can be easily detected and excluded, but the detection rate of small gross errors of G17 is significantly reduced by FDE algorithm.The recognition and exclusion rates for C24, C25, and G17 are 100%, 100%, and 30%, respectively.The solution results before and after the FDE algorithm used are shown in Figs. 13 and 14, and the statistical results are presented in Table 8.
Obviously, the gross error in the observations makes the solution result extremely worse, and the accuracy of position and velocity determination is obviously decreased, which reduces the reliability of autonomous navigation using GNSS.However, FDE algorithm can excluded gross errors and the solution results can be maintained in 3-sigma confidence intervals, even if the recognition ability of small gross errors is limited.

Navigation use GNSS with the adaptive algorithm
According to the above analysis, GNSS can provide continuous available signals for spacecraft in earthmoon space, which is suitable for spacecraft navigation.
In addition, considering the noise in GNSS observation, the proposed adaptive algorithms are also used in the autonomous navigation solution of the Earth-moon spacecraft.
The dynamic model and observation model parameters for the Kalman filter configuration settings are listed in Table 9.Because of the limited power of the onboard processor unit, and the higher-order gravity term has no significant effect on the accuracy of orbit determination for high and ultra-high orbits.We use gravity models of various orders at different locations in earth-moon space, which significantly reduces the amount of computation while maintaining similar navigation performance.10.We can  conclude that compared with S1, the position accuracy of S2 increased by 16.0% and the velocity accuracy increased by 38.7%.Additionally, more solutions exist in the 3-sigma range.S3 has the highest solution accuracy and the average position error of the solution is 99.97% within the 3-sigma.Compared with S1, the position and velocity accuracy increased by 28.3% and 44.3%, respectively, which proves the effectiveness of the proposed IBAE algorithm.Similarly, we used the same three strategies to determine the position and velocity of LTO.Additionally, we statistically calculated the accuracy of the solutions in three segments according to distance r between the spacecraft and the center of the Earth, where r < 20R E , 20R E < r < 40R E , and 40R E < r < 60R E , and R E is the radius of the Earth.The results are shown in Figs. 18,19 and 20 and Table 11.Similar to PHO, the navigation performance of S3 is also better than that of S2 and S1.Compared with S1, S3 has a position accuracy increase of 36.27%, the average probability of being within the range of 3-sigma is 99.97%, and the average probability of the velocity result within the range of 3-sigma is 100.00%.(The position result probabilities of being within the 3-sigma are 99.56%,99.57%, and 99.40%.Additionally, the velocity results are 99.99%,99.99%, and 100.00%, respectively) space deteriorates significantly, and continuous navigation services cannot be achieved if the threshold is set to 25 dB Hz.Therefore, the spacecraft needs to carry a receiver with a C/N 0 threshold of 15 dB Hz to navigate with GNSS in cis-lunar space, which is also a common view.For PHO, the requirement of the C/N 0 threshold is relatively low, and the navigation performance better than 15 m can be obtained even with a threshold of 0.25 dB Hz.
In conclusion, we found that GNSS can provide navigation services in most regions of earth-moon space through simulation experiments.Additionally, we verified the effectiveness of the proposed adaptive algorithm, which can further improve the performance of the navigation service of GNSS with an accurate dynamic model, and provide a reference for the use of GNSS in earth-moon space.By the way, it can be seen that there are obvious fluctuations in the curves of PHO and LTO solution results, which is caused by GNSS positioning geometry.

Integrated navigation with INS, GNSS and star tracker
During the long-term flight in deep space, the spacecraft will inevitably be disturbed by unknown acceleration that cannot be modeled (e.g., the change of the space environment or orbital maneuver).The accuracy of the dynamics model based on orbital mechanics will be significantly reduced, and the accuracy of the solution will be significantly deteriorated.Spaceborne accelerometers can make  up for the defects of orbital mechanics and provide an accurate dynamic model.We attempt to analyze the navigation accuracy of the dynamic model compensated for by the accelerometer when the spacecraft is subjected to unknown acceleration.We assume that the spacecraft made an orbital maneuver in LTO at 15:41:51.907 on April 24, 2007.An acceleration of (0.281096, − 0.194355, 0.0755011) m/s 2 was applied in the X, Y, and Z directions in b-frame for 10 s.At this time, the known dynamic model failed and the additional acceleration caused by this orbital maneuver can be obtained by on-board accelerometer.In addition, the moment of using accelerometer data can be determined according to the difference between accelerometer data and model acceleration, since normally the difference should not be significant.
Orbits 20 min before and 40 min after the moment of the maneuver were selected, and the accelerometer data were simulated according to the characteristics of the accelerometer in Table 4. Figure 23 shows the assumed real acceleration of the spacecraft, the model acceleration in the orbit filter, and the simulated accelerometer data.
It should be noted that the premise of using the proposed IBAE algorithm is an accurate dynamic model,   12.In Fig. 24.We found that relying only on the orbital dynamic model cannot maintain navigation accuracy during the orbital maneuver, and the position and velocity accuracies decrease rapidly, with position accuracy even below 500 m.However, when the accelerometer data are used instead of the dynamic model, navigation accuracy returns to normal.
Additionally, it is considered that the attitude of spacecraft can be determined using a gyroscope and star tracker.The spacecraft attitude determined by the star tracker alone and the gyroscope, star tracker integrated solution is provided in Fig. 26.
We found that the attitude determination accuracy of the single star tracker is better than 50 arcseconds, and the combined solution of the gyroscope and the star tracker can further improve accuracy and is better than 10 arcseconds.

Conclusions and discussion
Considering the need for autonomous navigation for future spacecraft, we explored the navigation performance of GNSS in earth-moon space.And an autonomous navigation system based on the GNSS receiver, IMU, and star tracker integrated system with orbit filtering is adopted.Additionally, we proposed the MBAE and IBAE algorithms, which we used in the case of GNSS ranging accuracy attenuation in earth-moon space.
Signal characteristics received by PHO and LTO orbiting spacecraft were analyzed, which shows that satellite navigation signals can cover most areas in Earth-moon space.The experiment shows that users can receive more BDS satellites than GPS, which may be due to the stronger power of the BDS sidelobe signal, and the higher orbital altitude of the BDS GEO and IGSO.
The accuracies of PHO in the X, Y, and Z directions using GNSS and orbital filter solving are 1.89, 6.27, and 5.09 m, respectively.Accuracy improved by 16.0% with the MBAE algorithm and further improved by 28.3% with the IBAE algorithm.A similar conclusion also be obtained for LTO, where the accuracies are 35.15, 20.92, and 30.21 m in the three directions when the spacecraft distance from the center of the Earth is greater than 40 Re, which proves that GNSS can provide navigation services in near-lunar space.We verified the effectiveness of the proposed adaptive algorithm.Additionally, we considered spacecraft maneuvers, for which position accuracy was substantially reduce than 500 m when determined using GNSS alone with the orbit filter; hence, we used accelerometers to compensate for the dynamic model, and improve the position and velocity accuracy during orbital maneuvers use INS, GNSS and star tracker integration navigation system.With the integration of a gyroscope and star tracker, an attitude accuracy of better than 10 arcseconds can be obtained.
Based on the above experiments and analysis, we proved the feasibility of autonomous navigation using GNSS for Earth-moon spacecraft, and then proposed an adaptive algorithm to effectively suppress GNSS observations noise, so that high-precision position of spacecraft can still be obtained in lunar space.Finally, the performance of INS, GNSS and star tracker integrated navigation system is verified, and the solution accuracy can be maintained during orbital maneuvers.
In the future, we will consider the integration of more sensors and navigation methods (e.g., the celestial navigation system) to further improve the autonomous navigation performance of spacecraft in and outside earth-moon space.

Fig. 4
Fig. 4 CE-1 lunar mission trajectory: the green and orange trajectories are PHO and LTO, respectively, and the yellow trajectory represents GNSS MEO and GEO.The orbit maneuver point is also marked

Fig. 5
Fig. 5 Geometry between GNSS satellites, spacecraft, and the Earth: spacecraft can receive GNSS signals from both the zenith and nadir directions

  100 ×Table 6 Fig. 11
Fig. 11 Attitude determination error of the simulated star tracker

Fig. 12
Fig. 12 Value of k and Chi-square distribution threshold in the process of solution and the success rate of elimination is counted

Fig. 13
Fig. 13 Results of position (left) and velocity (right) determination without FDE algorithm.(The three direction position result probabilities of being within the 3-sigma are 91.58%,93.19%, and 94.08%.Additionally, the velocity results are 94.61%, 96.67%, and 96.25%, respectively) We design three schemes to evaluate the performance of integrated navigation and the adaptive algorithm proposed in this study: Scheme 1 (short for S1): Navigation solution based on the Kalman filter with the constant covariance matrix (the diagonal elements of the covariance matrix are set to the constant 10 m for pseudo-range).Scheme 2 (short for S2): Navigation solution based on the Kalman filter with the MBAE algorithm.(For GPS, we set the coefficients of a and b in Eq. (35) to 0.1596 and 32.0745, respectively.Additionally, for BDS-3, we set the values to 0.0798 and 16.0373, respectively.)Scheme 3 (short for S3): Navigation solution based on the Kalman filter with the IBAE algorithm.The solution results for the position and velocity for the three schemes of PHO and are shown in Figs. 15, 16 and 17.Additionally, the statistical Root Mean Square (RMS) of the solutions are listed in Table

Fig. 16
Fig. 16 Results of the PHO position (left) and velocity (right) determination of S2.(The three direction position result probabilities of being within the 3-sigma are 99.94%,99.94%, and 99.93%.Additionally, the velocity results in all directions are 100.00%) Additionally, the left and right figures in Fig.21show the position results in PHO and LTO with the C/N 0 thresholds of 15, 20, and 25 dB Hz.And the statistical results of the position accuracy and the number of visible satellites with different C/N 0 threshold are shown in Fig.22.To analyze the performance of GNSS in lunar space, only the results at the end of the LTO ( 40R E < r < 60R E ) were counted.We found that, compared with the navigation performance of spacecraft in PHO, that of spacecraft in LTO is significantly affected by the C/N 0 threshold.If the threshold of C/N 0 increase to 20 dB Hz, the navigation accuracy near lunar

Fig. 21
Fig. 21 Navigation performance of spacecraft in PHO (left) and LTO (right) when the C/N 0 thresholds are set to 15, 20, and 25 dB Hz

Fig. 23
Fig. 23 Assumed real acceleration (green), model acceleration used in the orbital filter (blue), and acceleration measured by the accelerometer (red).Considering the scale limit, the change of the maneuvering acceleration is not displayed

Table 1
Initial state parameter settings for the simulation orbits

Table 2
Dynamic parameter settings for the simulation orbits

Table 3
Earth eclipse angle and main lobe angle of BDS-3 and GPS

Table 4
Parameter settings for the simulation observations

Table 7
Detailed configuration for GNSS gross error

Table 8
Statistical values of solution results with and without FDE algorithm SchemePosition in different direction (m) Velocity in different direction (m/s)

Table 9
Parameter settings ( r is the distance from the Earth's center)

Table 10
Statistical values of the PHO solution results for the three

Table 11
Statistical values of the LTO solution results for the three schemes

Table 12
Statistical results for the two schemes

Scheme Position in different direction (m) Velocity in different direction (m/s)
Fig. 26 Attitude solution results determined by a single star tracker and gyro/star tracker integrated system