An Improved Augmented Algorithm for Direction Error in XPNAV

: Recently, X-ray pulsar-based navigation (XPNAV) as a signiﬁcant navigation method has been widely used in deep space exploration. However, the accuracy of XPNAV is limited to the existence of the pulsar direction error. To improve the performance of XPNAV, we have proposed a novel algorithm named “the modiﬁed augmented state extended Kalman ﬁlter” (MASEKF). The algorithm considers the high-order terms of direction error and then adds a more precise direction error into state equation and measurement equation. In the simulation, by comparing the performance of MASEKF, EKF, and ASEKF at the same time, it is found that MASEKF has better performance in the accuracy and stability, and the results also demonstrate that MASEKF algorithm has faster convergence speed. This paper provides a strong reference for other improvements of algorithms towards direction error. The purpose of this study is to establish MASEKF and add the direction error into the measurement equation and the state equation, so as to realize the coordination and symmetry of the algorithm.


Introduction
Over the past decades, the Global Positioning System (GPS), a well-known navigation method, has been performing well in low Earth orbit (LEO) satellite navigation [1]. It provides users with real-time position and speed information. Unaffected by the weather, the positioning accuracy of the single machine is better than 10 m. The dominant application and high precision of GPS have strongly impressed people. However, the GPS signal is designed to transmit information toward the earth; the dominating advantage of GPS cannot be used well in in-depth space exploration. Thus, GPS cannot provide autonomous navigation information services for spacecraft flying in interstellar space [2]. To overcome the above defects, X-ray pulsar-based navigation (XPNAV) is emerged that can work well in solving deep space problems; it uses the X-ray emitted by pulsars in the universe, then obtains the precise information of spacecraft by measuring the Time-of-Arrival (TOA) of X-ray pulsar radiation signals [3]. In XPNAV, the X-ray detector is installed on the spacecraft, the radiated X-ray photons are detected by the pulsar, the photon arrival time is measured, and the pulsar image information is extracted. After the corresponding signal and data processing, the spacecraft determines some parameters, such as the orbit of navigation, time, and attitude. As a matter of fact, XPNAV has attracted plenty of attention around the world.
In the previous researches of XPNAV, a significant research achievement is realized by proposing a recursive extended Kalman filter algorithm (EKF). Through the application of EKF, the accuracy of spacecraft's position error can be controlled within 100 m, and the accuracy of velocity error is around 10 mm/s [4]. On the basis of EKF algorithm, a novel algorithm named the suboptimal multiple fading extended Kalman filter (SMFEKF) is established subsequently, which can be up to 65.43% and 61.70%

The Principle of XPNAV
XPNAV is well suited in dealing with deep space exploration problems, and pulsar used in navigation can emit X-ray photon signals. During the process of the propagation of the signal, X-ray detector on the spacecraft can capture the signal to measure the amount and the arriving time of photon [16,17]. Subsequently, epoch folding method is used (it is the process of counting the recorded photons by means of direct folding or rapid folding and finally, normalizing the photon number and obtaining the pulse contour) [18] to obtain the time t sc of the pulsar arrival at the spacecraft, using a pulsar timing method to get the time t SSB of the pulsar arrival at the center of the solar system (SSB). Then, XPNAV can calculate corresponding time delay between t sc and t SSB [19]. Additionally, the corresponding distance can be acquired, which is the difference between the distance from SSB to pulsar and the distance from spacecraft to the pulsar. In addition, in spite of the tiny changes of pulsars, during the simulation period, when the position of the pulsar and the spacecraft's initial parameters have been determined, a recursive estimation model of spacecraft position and velocity can be established accurately [20,21]. The basic principle of XPNAV is shown in Figure 1.

The Principle of XPNAV
XPNAV is well suited in dealing with deep space exploration problems, and pulsar used in navigation can emit X-ray photon signals. During the process of the propagation of the signal, X-ray detector on the spacecraft can capture the signal to measure the amount and the arriving time of photon [16,17]. Subsequently, epoch folding method is used (it is the process of counting the recorded photons by means of direct folding or rapid folding and finally, normalizing the photon number and obtaining the pulse contour) [18] to obtain the time t of the pulsar arrival at the spacecraft, using a pulsar timing method to get the time of the pulsar arrival at the center of the solar system (SSB). Then, XPNAV can calculate corresponding time delay between and [19]. Additionally, the corresponding distance can be acquired, which is the difference between the distance from SSB to pulsar and the distance from spacecraft to the pulsar. In addition, in spite of the tiny changes of pulsars, during the simulation period, when the position of the pulsar and the spacecraft's initial parameters have been determined, a recursive estimation model of spacecraft position and velocity can be established accurately [20,21]. The basic principle of XPNAV is shown in Figure 1. Following the above analysis, it can be found that and are two fundamental parameters that need to be determined at the beginning, which is usually calculated by the epoch folding method. Another primary parameter ∆ can be easily determined by ∆ = − . The result of the system error is not negligible for the positional accuracy of the spacecraft. Thus, the pulse time model can be modified as follows [4,9].
where is the time from the ℎ pulsar to the SSB, is the time from the ℎ pulsar to the spacecraft, is the speed of light, stands for the real direction vector of the ℎ pulsar, is SSB Pulsar 3 Pulsar 2 Pulsar 1 Earth Spacecraft Figure 1. The principle of X-ray pulsar-based navigation (XPNAV) (the parameters and geometric diagrams used by XPNAV).
Following the above analysis, it can be found that t SC and t SSB are two fundamental parameters that need to be determined at the beginning, which is usually calculated by the epoch folding method. Another primary parameter ∆t can be easily determined by ∆t = t SSB − t SC . The result of the system error is not negligible for the positional accuracy of the spacecraft. Thus, the pulse time model can be modified as follows [4,9]. where t i SSB is the time from the ith pulsar to the SSB, t i SC is the time from the ith pulsar to the spacecraft, c is the speed of light, n i stands for the real direction vector of the ith pulsar, r SC is the vector from SSB to spacecraft, D 0 represents the relative distance between the pulsar and the SSB, b is the position vector of SSB relative to the center of the sun, G is the gravitational constant, M sun is the mass of the sun, and B i is the system bias caused by the direction error of the ith pulsar. On the right-hand side of the expression, the first item is the first-order Doppler delay, which is the propagation time from the detector to the SSB and the pulsar to the SSB. The second item is the annual parallax caused by the annual change in the earth orbit, which together with the first term are collectively referred to as Roemer delay, and the third item is Shapiro delay [22].
As is indicated in Figure 1, since GPS satellite is used in the simulation, the orbital dynamics model is required to establish under the Earth-centered inertial (ECI) reference frame. Thereby the relationship of the vector between spacecraft, earth centroid, and SSB can be expressed below.
where r sc represents the vector from SSB to spacecraft, r i is the vector between the ith spacecraft and earth center of mass, and r e stands for the vector between SSB and earth center of mass, which is usually obtained by JPL ephemeris DE421.

Analysis of System Error
In the process of applying XPNAV, system error was produced by the following aspects: 1. direction error ∆n, in this paper, Roemer delay and Shapiro delay were the main source of direction error; 2. clock error caused by clock drift; 3. the error caused by D 0 ; and 4. position error caused by ephemeris. In general, the ephemeris error is time varying and changes slowly. However, the ephemeris error can increase by 0.25 m within 1 h, which should be taken into account seriously in XPNAV [23]. Among the above four common errors, the error caused by D 0 is tiny, that is because the distance between the pulsar and the SSB is quite large [14,24]. Clock error is another vital item that do affect the accuracy of navigation; in the subsequent discussion, its value and possible changed situation are carefully considered.
In the prior studies, the Time-Differenced method is often used to handle system error. However, due to the inaccuracy of the initial model and without the consideration of direction error, after applying the Time-Differenced method in the Mars probe and Jupiter detector, it is found that the position error can even reach to the level of kilometers. Besides, as a commonly used algorithm UKF, constructing sampling points is easy to diverge and always accompany with poor stability. Some existing methods could not deal with the existence of large system error well [25]. To research the system error more scientifically, the other error terms causing system error should be discussed in detail, and the corresponding algorithm should be optimized to meet the higher accuracy requirements. B i is obtained by the difference on the right-side term of Equation (1). In the difference, the corresponding terms of the real direction vector are different from the terms of the estimated direction vector, and after the elimination of the terms independent of the direction vector, Equation (3) is obtained as below.
wheren i is the estimated unit direction vector. In Equation (3), many terms are defining B i , but it is found that most existing algorithms only consider the effect of the first term n i −n i · r sc . In fact, all terms in B i should be considered seriously; thus, a simple calculation is conducted to get the magnitude of every item in Equation (3); the specific Symmetry 2020, 12, 1059 5 of 23 position in the orbit of the ARGOS vehicles was chosen at time = 2451538.96769266 JD, and the results were only valid for this specific location. The single source used in this analysis was the Crab Pulsar, along with its proper motion of µ α = −17 mas/year and µ δ = 7 mas/year and a distance of 2 kpc [26]. A value of ∆t n = 100 days was selected as a representation of the elapsed time between the measured D 0 and the current time [18].
In Table 1, it can be seen that the value of the first term is quite larger than that of other terms, but that does not mean that the term with smaller value can be discarded casually; it is found that even the third term with the smallest value can produce a position error around 0.45 m (0.000000001493c ≈ 0.45 m). Therefore, all terms in B i need to be handled carefully. In order to capture the subtle changes produced by B i , the derivative of B i versus time t is calculated as is shown below.  In accordance with Equation (4), it can be found that n i −n i plays a vital role in dB i dt , ∆n = n i −n i , and ∆n is the pulsar direction error as well. dr e dt is the earth rotation rate; due to the existence of direction error and the earth rotation, the accompanied deviation can be quite large. For example, if the pulsar direction error is 0.001 , the maximum value of the earth rotation rate is closed to 30.3 km/h and the maximum of deviation is approximate to 1.468 × 10 −4 m/s [27]. Therefore, the algorithm towards the pulsar direction error should be optimized further.

Analysis of Direction Error
According to literature [27], the Taylor series expansion of direction error should be discussed. Taylor's series expansion method is widely used in the above arguments, but its quadratic terms and other terms are usually ignored [28], which is a primary reason because of which the direction error cannot be mainly eliminated.
In this part, the right ascension angle α and the real declination angle δ form the pulsar direction vector; direction vector is usually shown as (α, δ) . Let the estimated position be α,δ , then the pulsar direction error ∆n can be calculated using Usually, the real pulsar direction vector n i is expressed as follows.
In the prior researches, the high order terms in n i are ignored consistently; only the first-order term is saved in the final algorithm. Its expression is usually shown as follows.
In this paper, the high-order expansion terms of direction error are processed further; the second-order term of the Taylor series is shown below.
Actually, the apparent differences between the previous methods and the proposed method need to be shown clearly. We introduced the measurement equation and the state equation into discussions because these two equations can reflect the differences between various ways straightly.

The Improved Algorithm of MASEKF
Extended Kalman filter (EKF) is an efficient recursive filter that estimates the state of a dynamic system based on a series of measurements. It performs the first-order linear truncation of Taylor expansion of the nonlinear function but ignores the remaining higher order terms. Therefore, EKF is applicable for the nonlinear system [29]. However, ignoring the high-order terms does affect the navigation accuracy. In order to overcome the above defect, ASEKF is proposed by considering the state noise and measurement noise in the augmented state.
In the XPNAV, based on the dynamic equation of spacecraft in the Earth-centered inertial, to construct the state equation and measurement equation, several specific states as the position and velocity of the probe are selected to occupy in the equation, thereby the ordinary state equation of ASEKF is calculated as follows [30]. .
where r = (r x , r y , r z ) is the position vector of the spacecraft in ECI, r = r 2 is the velocity vector of the spacecraft in ECI, ∆n i is the norm of pulsar direction error, µ = GM 0 is the gravity of the earth, which can be viewed as a fixed constant. J 2 is the second-order harmonic Symmetry 2020, 12, 1059 7 of 23 coefficient and R e is the radius of the earth. ∆F x , ∆F y , ∆F z are the higher-order perturbation terms, such as solar radiation pressure perturbation, etc., which can be regarded as zero-mean Gaussian white noise. λ 1 , λ 2 , λ 3 are the process noises in the rate of ∆n, which are normally handled as zero-mean Gaussian white noise as well.
It can be found that EKF improves the accuracy of Kalman filtering, then the accuracy of EKF is improved further by ASEKF, whereas ASEKF algorithm is optimized by adding the direction error terms, which has already been widely used in navigation around the world [31,32]. To improve the accuracy of ASEKF further, the other high-order terms existed in direction error should be considered seriously.

The State Equation
In this paper, MASEKF was proposed to meet higher accuracy requirements, which took the high-order terms of direction error into account seriously. As is indicated in Equation (8), the direction error has been expanded in second-order terms by Taylor series. To get the expression of the state equation of MASEKF, the formation of direction error is changed as follows.
where α is the right ascension angle, ∆α i is the first-order differential of right ascension, ∆α i 2 is the second-order differential of right ascension, δ is the real declination angle, ∆δ i is the first-order differential of the real declination angle, and ∆δ i 2 is the second-order of differential of the real declination angle. Then, combining Equation (10), the state equation of MASEKF is subsequently obtained. It can be shown as follows. .
Symmetry 2020, 12, 1059 8 of 23 where F(t) is the Jacobian matrix of the system, J 2 is the second-order harmonic coefficient, R e is the radius of the earth, µ e is the gravity of the earth, λ 1 and λ 2 stand for the change rate of the right ascension and declination of the first pulsar, λ 3 and λ 4 indicate the change rate of the right ascension and declination of the second pulsar, λ 5 and λ 6 show the change rate of the right ascension and declination of the third pulsar, ω is the process noise, and both λ and ω are the zero-mean Gaussian white noise.

The Measurement Equation
The changes in velocity vector and position vector were both considered in the MASEKF algorithm. The measurement model is determined as follows.
, m is the number of the selected pulsars, and V(t) is the measurement noise.
By combining Equation (1) and Equation (12) in the meantime, the measurement equation can be subsequently determined, which is calculated as follows.
where h(X, t) is the right-hand side of Equation (13). When n i =n i + ∆n i , let the second item and the third item on the right side of Equation (18) be high-order term ε. Then, Equation (18) can be transformed as follows. where ∆t is the difference between t i SSB and t i SC . Though the high-order term ε is of a rather small magnitude, in order to avoid the truncation error as much as possible, it still should be linearized subsequently.
In the process of kth filter calculation, the magnitude ofn i is roughly equal to that of n i + ∆n i , hence, the term ε n i is used to replace the term ε n i + ∆n i . Equation (19) can be rewritten as follows.
where ε 1 k n i stands for the second term in Equation (18) and ε 2 k n i stands for the third term in Equation (18); they can be expressed as Since the position of spacecraft in BCRS (Barycentric Celestial Reference System) can be viewed as same as that of the earth in BCRS, compared with the spacecraft's position, the changes of spacecraft's velocity and position were quite small. Thus, it is reasonable to let part of r sc(k) be replaced by r sc(k−1) , then Equations (22) and (23) can be expressed as follows.
Subsequently, the second-order Taylor expansion to r SC(k) at time k − 1 was carried out. In order to make the process of expansion convenient, some simplification is done as follows.

The Procedures of Simulation
In order to further discuss the effectiveness and performance of the proposed algorithm, EKF and ASEKF were also simulated at the same time, and then all simulation results were analyzed in detail. In this paper, the procedures of the simulation were divided into four parts as Figure 2. To make this research more comprehensive, the performance of three different algorithms was verified further. Root mean square error (RMSE) was applied to confirm the simulation precision of three algorithms, which can quantify the deviation between the observed values and real values straightly. It is worth noting that RMSE is extremely sensitive to outliers, such as maximum error and minimum error. The expression of RMSE is as follows.
where is the of the estimated position error, is the RMSE of the estimated velocity error, is the sampling point of the stable filter interval, ∆ is the estimated position error, is the initial moment of filtering stability, and ∆ is the estimated velocity error. In the end, ℎ of the three algorithms was calculated and analyzed cautiously.

Preprocessing
The work of preprocessing should be completed initially, which includes obtaining the parameters of pulsar, spacecraft, etc. The required parameters were determined through some official websites, databases, and software. In particular, the data of the pulsar's position and velocity is obtained in the Australia Telescope National Facility (ATNF). In this study, Pulsar Catalogue was used to generate observation data, which is viewed as the real position as well [14]. It is worth mentioning that the pulsar position used in different algorithms is the real pulsar position plus the set error. Meanwhile, the position vector from the earth center of mass to SSB was obtained from DE421; STK software provided the necessary information of probe, such as position and velocity.

Observation Simulation
The main work of this part was to obtain the proper pulsars used in later simulation; all pulsars' parameters were inputted into Inter Visual Fortran 11. Generally speaking, the higher the signal frequency of the pulsar, the greater the radiation energy of the photon, the smaller the background radiation energy, the narrower the pulse width in the periodic contour, the larger the number of pulses, the larger the quality factor of the pulsar, and the higher the positioning accuracy of the pulsar. Under the constraints of positioning accuracy and quality factor, Weighted Dilution Of Precision (WDOP) method and observation analysis to acquire the most suitable pulsars were used [21]. Specifically, after selecting pulsars, the most important thing was to determine the time delay. The number of photons reaching the detector can be obtained by simulating the number of photons generated in Poisson distribution. The pulse profile of simulated observation can be obtained by adding the number of photons repeatedly. Then, using the FFT method to cross-correlate the simulated observation profile with the standard pulse profile, the time delay can be obtained in the end.

Analysis of System Error
In this research, the direction error caused by Roemer delay and Shapiro delay was analyzed in detail.

4.
Orbit Determination The core work involved in this part was conducted based on the MASEKF module, combining orbital dynamics to process the simulated observation data. The state of the spacecraft was obtained, and the precision statistics of the calculated results were carried out. When the precision evaluation exceeds the threshold value of Cramer-Rao Lower Bound [33], the simulation was finished, otherwise, the orbit simulation was simulated again. The experiment cannot stop simulating until the accuracy meets the requirements.
To make this research more comprehensive, the performance of three different algorithms was verified further. Root mean square error (RMSE) was applied to confirm the simulation precision of three algorithms, which can quantify the deviation between the observed values and real values straightly. It is worth noting that RMSE is extremely sensitive to outliers, such as maximum error and minimum error. The expression of RMSE is as follows.  (32) where RMSE r is the RMSE of the estimated position error, RMSE v is the RMSE of the estimated velocity error, M is the sampling point of the stable filter interval, ∆r is the estimated position error, t ko is the initial moment of filtering stability, and ∆v is the estimated velocity error. In the end, the RMSE of the three algorithms was calculated and analyzed cautiously.

Initial Conditions
Actually, there are 140 pulsars that can be used in navigation; in this paper, WDOP method was applied with the constraint of quality factor to select 40 pulsars, which is the process of preliminary selection as well. If the selected 40 pulsars have better performances, then the most suitable pulsars are chosen through observation analysis. Finally, three pulsars (B1821-24, B1937+21, and B0531+21) were adopted as the navigation sources. The corresponding parameters of the three pulsars are shown in Table 2. In the simulation, the research object was the earth-orbiting spacecraft. The parameters are presented in Table 3. Meanwhile, it is well accepted that GPS plays a vital role in the navigation system, and its data is available [13,25]; many prior studies have selected GPS as the research object due to its importance and extensive application. Thus, GPS_BII-10 was adapted to be further researched in this paper; its corresponding data are provided by STK software. The simulation time was set between 2017.10.1.16.00.00 and 2017.10.2.16.00.00. The measurement time step was determined as 246.8 s, whose effectiveness has been confirmed many times in the existing pieces of literature [3,19]. Usually, the filtering step runs immediately after the pulsar observation step [34]. The total time of the simulation was 86,400 s; the simulation parameters of EKF, ASEKF, and MASEKF were set as follows.  In addition, the observation noise covariance can be determined as follows [18].
where, w is the width of the pulsar, B x is the X-ray background radiation flux, which is about 0.005 ph/cm 2 /s in the prior researches, F x is the radiation photon flux from the pulsar, P f is the ratio of the pulse radiation flux to the average radiation flux in one pulsar period, S is the effective area of the X-ray detector, its value is 1 m 2 , t obs is the observation period, t obs = 246.8 s, and d is the ratio of w to P. The pulsar period parameters are shown in Table 2. Hence, the covariance of observation noise can be calculated. The result of R is given by R = diag 232.23 2 , 247.02 2 , 77.78 2 (34) where diag is an abbreviation for diagonal matrix.

Simulation and Analysis
In the simulation, in order to confirm the performance of the proposed algorithm clearly, ASEKF and EKF were simulated to be compared with MASEKF, conducting the simulation under three different conditions. In this paper, the factor of the direction error is an essential changeable condition of the simulation, since the system error includes the direction error, clock error, ephemeris error, and other possible existing errors, the settings of the above errors are presented below.

2.
Ephemeris error According to reference [35], the ephemeris error is usually less than 1 × 10 −13 in the other planets' ephemerides. However, the error in the earth's ephemeris has been verified to be the biggest, which could have specific influences on the performance of XPNAV. Since the simulation time is from 2017.10.1.16.00.00 to 2017.10.2.16.00.00, and ephemeris is time varying. Consequently, the changeable tendency of ephemeris error within the simulation period can be determined in Figure 3. Symmetry 2020, 12, x FOR PEER REVIEW 14 of 21

Direction error
Owing to the limited space and the previous attempts in other works of literature [34], the simulation was conducted with the initial error condition of (0 mas, 0 mas), (1 mas, 1 mas), and (−1 mas, −1 mas); it should be noted that mas as a unit represents milli-arcsecond.

Cosmic background noise
In the observation and simulation stage, the photon arrival model under the Poisson distribution to simulate the distribution of cosmic background noise was construct; it was found that the pulse observation profile of the photon arrival model under the Poisson distribution was closest to the actual pulse profile. In the orbit determination stage, cosmic background noise was added to the observation noise, and then it was subjected to the Kalman filter algorithm for processing.

Other errors
According to the prior researches, other errors may mainly exist in the following aspects: (1) the distance between SSB and spacecraft, the greater the distance, the higher the system error; (2) the value of the pulsar frequency derivative, the larger the value, the larger the system error, and (3) the size of proper motion and the age of the position epoch, the greater the appropriate motion, the older the epoch, the greater the system error. Considering the above possible existing errors, other errors were set to 10 .
Besides, the simulated observation data was combined with random noise. By adding the random noise during the simulation process of obtaining and , whose value is usually less than 10 , the corresponding position error can be controlled within 30 . Besides, the uncertainty information (the triple standard deviation) was added to prove the reliability of the estimated results. The data are indicated in Figures 4-6.

Direction error
Owing to the limited space and the previous attempts in other works of literature [34], the simulation was conducted with the initial error condition of (0 mas, 0 mas), (1 mas, 1 mas), and (−1 mas, −1 mas); it should be noted that mas as a unit represents milli-arcsecond.

4.
Cosmic background noise In the observation and simulation stage, the photon arrival model under the Poisson distribution to simulate the distribution of cosmic background noise was construct; it was found that the pulse observation profile of the photon arrival model under the Poisson distribution was closest to the actual pulse profile. In the orbit determination stage, cosmic background noise was added to the observation noise, and then it was subjected to the Kalman filter algorithm for processing.

5.
Other errors According to the prior researches, other errors may mainly exist in the following aspects: (1) the distance between SSB and spacecraft, the greater the distance, the higher the system error; (2) the value of the pulsar frequency derivative, the larger the value, the larger the system error, and (3) the size of proper motion and the age of the position epoch, the greater the appropriate motion, the older the epoch, the greater the system error. Considering the above possible existing errors, other errors were set to 10 −8 m.
Besides, the simulated observation data was combined with random noise. By adding the random noise during the simulation process of obtaining t SC and t SSB , whose value is usually less than 10 −7 s, the corresponding position error can be controlled within 30 m. Besides, the uncertainty information (the triple standard deviation) was added to prove the reliability of the estimated results. The data are indicated in Figures 4-6. Besides, the simulated observation data was combined with random noise. By adding the random noise during the simulation process of obtaining and , whose value is usually less than 10 , the corresponding position error can be controlled within 30 . Besides, the uncertainty information (the triple standard deviation) was added to prove the reliability of the estimated results. The data are indicated in Figures 4-6.  In Figures 4-6, it can be found that the estimated result of the position error fluctuates greatly, especially in the case of (0 mas, 0 mas), the estimated result of EKF's position error cannot satisfy the requirement of the triple standard deviation around 50 epochs, which indicates the poor performance of EKF. In contrast, the estimated results of ASEKF and MASEKF's position error are acceptable in  In Figures 4-6, it can be found that the estimated result of the position error fluctuates greatly, especially in the case of (0 mas, 0 mas), the estimated result of EKF's position error cannot satisfy the requirement of the triple standard deviation around 50 epochs, which indicates the poor performance of EKF. In contrast, the estimated results of ASEKF and MASEKF's position error are acceptable in scenarios of (0 mas, 0 mas) and (1 mas, 1 mas), it is notable that the estimated result of velocity error In Figures 4-6, it can be found that the estimated result of the position error fluctuates greatly, especially in the case of (0 mas, 0 mas), the estimated result of EKF's position error cannot satisfy the requirement of the triple standard deviation around 50 epochs, which indicates the poor performance of EKF. In contrast, the estimated results of ASEKF and MASEKF's position error are acceptable in scenarios of (0 mas, 0 mas) and (1 mas, 1 mas), it is notable that the estimated result of velocity error of three algorithms can satisfy the requirement of the triple standard deviation, which can strongly illustrate the effectiveness of ASEKF and MASEKF under any simulation condition. Table 4 displays the statistical simulation results, when there is no initial error; it is apparent that EKF has the worst performance. It is well accepted that EKF performs Kalman filtering by linearizing the nonlinear system, ignoring the direction error term unconsciously. This means that EKF cannot eliminate the corresponding measurement noise and state noise effectively; the above defects of EKF result in its lower accuracy. Through analyzing Figure 4, it was found that the velocity error of EKF diverges at 300~350 epochs, which was rather unstable, thereby it was tough for the system to converge on the whole. What's more, different from EKF algorithm, ASEKF algorithm linearizes the nonlinear system after adding the direction error into the state vector, considering the first-order term of the direction error in Taylor series expansion, which could certainly improve the accuracy in a way. Therefore, the accuracy of ASEKF has improved a lot than that of EKF, but the performance of ASEKF seems to be a little bit of unstable, as indicated in Figures 5 and 6. It could be seen that the accuracy of ASEKF decreases quickly in scenarios of (−1 mas, −1 mas); the accuracy difference between (−1 mas, −1 mas) and (1 mas, 1 mas) is much bigger than that of between (0 mas, 0 mas) and (1 mas, 1 mas), and it seems that ASEKF is more suited in scenarios of (0 mas, 0 mas) and (1 mas, 1 mas).
Furthermore, based on the theory of the ASEKF algorithm, the proposed MASEKF has considered the second-order terms of direction error, adding the more accurate direction error into the state equation and measurement equation subsequently. Therefore, the MASEKF algorithm tends to have a better convergence effect and stronger stability. As is shown in Table 4, it is evident that MASEKF has the best accuracy performance under three different conditions; the accuracy of ASEKF is close to the accuracy of MASEKF in scenarios of (1 mas, 1 mas). Besides, MASEKF can converge quickly to 121 m in position error, whereas the velocity error of MASEKF converges rapidly to about 0.729 m/s after about 50 epochs. It seems that the improvement and optimization of MASEKF only bring a small increase in accuracy in scenarios of (1 mas, 1 mas). In fact, when the simulation are conducted in situations of (−1 mas, −1 mas), the advantages of MASEKF stand out, the accuracy of MASEKF is quite higher than the other two algorithms, and its accuracy difference is quite small between (−1 mas, −1 mas) and (1 mas, 1 mas), which verifies the better stability of MASEKF.
What's more, for the purpose of comparing the performance of different algorithms further, simulations were conducted to obtain the RMSE under different conditions; RMSE can reflect the difference of the accuracy straightly. Figures 7-9 show the performance comparison of EKF, ASEKF, and MASEKF under different initial conditions. between (−1 mas, −1 mas) and (1 mas, 1 mas), which verifies the better stability of MASEKF.
What's more, for the purpose of comparing the performance of different algorithms further, simulations were conducted to obtain the RMSE under different conditions; RMSE can reflect the difference of the accuracy straightly. Figures 7-9 show the performance comparison of EKF, ASEKF, and MASEKF under different initial conditions.  For clarity, Table 5 lists the necessary information of Figures 7-9. Regardless of whether there was system bias, MASEKF always had the fastest converge speed and the best navigation accuracy. The performance of ASEKF was much better than that of EKF, which illustrated that the optimization of the high-order terms of the direction error significantly improved the accuracy of the algorithm. Besides, it can also be observed that the proposed algorithm considered second-order terms of Taylor series of direction error, which outperforms the algorithm that only considered the first-order terms of Taylor series of direction error. Furthermore, as is shown in Table 6, the accuracy of EKF and ASEKF was analyzed compared to that of MASEKF. It was found that EKF has poor accuracy, especially in scenarios of (0 mas, 0 mas), and the accuracy of the velocity error of MASEKF was improved by 43.06% than that of EKF. Similarly, the accuracy of the velocity error of MASEKF was improved by 36.11% than that of ASEKF. Meanwhile, though the accuracy of ASEKF improved a lot than the traditional EKF algorithm, its performance cannot stay stable in all conditions, such as in scenarios of (1 mas, 1 mas); the accuracy of position error and velocity error of ASEKF were worse than that of EKF. All values were positive in Table 6, which means that the accuracy of EKF and ASEKF were worse than that of MASEKF. In fact, MASEKF always performed well, compared with traditional EKF algorithm, and it increased the average accuracy by 8.84% and 20.53% in position error and velocity error, respectively, under the three conditions. For clarity, Table 5 lists the necessary information of Figures 7-9. Regardless of whether there was system bias, MASEKF always had the fastest converge speed and the best navigation accuracy. The performance of ASEKF was much better than that of EKF, which illustrated that the optimization of the high-order terms of the direction error significantly improved the accuracy of the algorithm. Besides, it can also be observed that the proposed algorithm considered second-order terms of Taylor series of direction error, which outperforms the algorithm that only considered the first-order terms of Taylor series of direction error. Furthermore, as is shown in Table 6, the accuracy of EKF and ASEKF was analyzed compared to that of MASEKF. It was found that EKF has poor accuracy, especially in scenarios of (0 mas, 0 mas), and the accuracy of the velocity error of MASEKF was improved by 43.06% than that of EKF. Similarly, the accuracy of the velocity error of MASEKF was improved by 36.11% than that of ASEKF. Meanwhile, though the accuracy of ASEKF improved a lot than the traditional EKF algorithm, its performance cannot stay stable in all conditions, such as in scenarios of (1 mas, 1 mas); the accuracy of position error and velocity error of ASEKF were worse than that of EKF. All values were positive in Table 6, which means that the accuracy of EKF and ASEKF were worse than that of MASEKF. In fact, MASEKF always performed well, compared with traditional EKF algorithm, and it increased the average accuracy by 8.84% and 20.53% in position error and velocity error, respectively, under the three conditions. Explanation: ∆ represents the difference between the compared algorithm and MASEKF, "Percentage" means convert ∆ to a percentage, the larger its percentage value, the lower its accuracy.
To confirm the stability of the proposed algorithm in the whole simulation period, we compared the simulated orbit with the real one. The simulated trajectories and their partial enlargement are shown in Figure 10 [25]. It is evident that the MASEKF algorithm has the closest approximation between the simulated track and the real track at any simulation time. All results demonstrated that MASEKF not only had obvious superiority in the overall result but also ensured its high accuracy and stability at any time.  Explanation: ∆ represents the difference between the compared algorithm and MASEKF, "Percentage" means convert ∆ to a percentage, the larger its percentage value, the lower its accuracy.
To confirm the stability of the proposed algorithm in the whole simulation period, we compared the simulated orbit with the real one. The simulated trajectories and their partial enlargement are shown in Figure 10 [25]. It is evident that the MASEKF algorithm has the closest approximation between the simulated track and the real track at any simulation time. All results demonstrated that MASEKF not only had obvious superiority in the overall result but also ensured its high accuracy and stability at any time.

Conclusions
In this paper, MASEKF algorithm is proposed based on the ASEKF algorithm, which expands the second-order term of the direction error to reduce the high-order truncation errors. Adding the

Conclusions
In this paper, MASEKF algorithm is proposed based on the ASEKF algorithm, which expands the second-order term of the direction error to reduce the high-order truncation errors. Adding the more precise direction error into the state equation and measurement equation, thus, MASEKF can reduce the non-Gaussian error propagation in the nonlinear system. Finally, simulation is conducted by comparing the performance of MASEKF, ASEKF, and EKF. It is found that MASEKF is more versatile and efficient, compared with traditional EKF algorithm, and MASEKF can even increase the average accuracy by 8.84% and 20.53% in position error and velocity error, respectively, in the three simulation conditions. In the future, scholars can research the reliable algorithm of pulsar position error or other errors according to the algorithm proposed in this paper.