X-Ray Pulsar-Based Navigation Considering Spacecraft Orbital Motion and Systematic Biases

The accuracy of X-ray pulsar-based navigation is greatly affected by the Doppler effect caused by the spacecraft orbital motion and the systematic biases introduced by the pulsar directional error, spacecraft-borne clock error, etc. In this paper, an innovative navigation method simultaneously employing the pulse phase (PP), the difference of two neighbor PPs (DPP) and the Doppler frequency (DF) of X-ray pulsars as measurements is proposed to solve this problem. With the aid of the spacecraft orbital dynamics, a single pair of PP and DF relative to the spacecraft’s state estimation error can be estimated by using the joint probability density function of the arrival photon timestamps as the likelihood function. The systematic biases involved to the PP is proved to be nearly invariant over two adjacent navigation periods and the major part of it is eliminated in the DPP; therefore, the DPP is also exploited as additional navigation measurement to weaken the impact of systematic biases on navigation accuracy. Results of photon-level simulations show that the navigation accuracy of the proposed method is remarkably better than that of the method only using PP, the method using both PP and DF and the method using both PP and DPP for Earth orbit.


Introduction
Pulsars are strongly magnetized and rapidly rotating neutron stars emitting signals that are unique and highly periodical [1,2]. Along with their spatial diversity about the galactic disk, pulsars have been suggested as a natural lighthouse for spacecraft navigation. Among all pulsars, X-ray pulsar are more suitable for spacecraft navigation due to the much smaller detector size required compared to that for radio or optical sources [3]. X-ray pulsar-based navigation (XPNAV) has been suggested as a potential approach for autonomous spacecraft navigation [4][5][6]. Unlike some current applications for spacecraft navigation, such as Deep Space Network and Global Navigation Satellite System, which suffer from low performance outside their effective coverage and rely extensively on ground-based operations, the XPNAV has the same accuracy throughout the solar system and has much more autonomy [7,8]. It can also be used to augment the current navigation systems to improve their performance by introducing pulsar measurements [9,10].
In XPNAV, a key task is the estimation of pulse phase (PP), which requires conversion of the measured photon time-series to the coordinate time at the Solar System Barycenter (SSB) using the initial position estimate of the spacecraft and various ephemeris parameters [6,11,12]. Many studies of XPNAV, such as [13,14], have assumed the spacecraft position is invariable within a filtering period and have no consideration of the ephemeris errors and the satellite-borne clock error. Reference [13]

Signal Models of X-Ray Pulsars
The fundamental data provided to the XPNAV system are the observed TOAs of X-ray photons from the pulsar sources and cosmic X-ray background. The foregone photon arrival model of X-ray pulsars is established at a hypothetical reference location usually chosen as the SSB. The photon arrival model at the spacecraft can be obtained by transferring the one at the SSB to the spacecraft according to the spacecraft's position and velocity.

Photon Arrival Rate Function at the SSB
Assume there exists an X-ray detector at the SSB. Let (t 0 , t f ) be the observation interval. Let t i be the TOA of the i-th photon, and the set {t 0 , t 1 , · · · , t K }, denoted by {t i } K i=1 , be a random sequence in increasing order. The number of the detected photons, K, is also a random variable. Let t 0 = 0, N 0 = 0, and N t be the number of detected photons in the interval (0, t). The detected X-ray photon event timestamps {N t , t > 0} can be modeled as the arrival times of a Non-Homogeneous Poisson Process (NHPP) with a time varying rate function λ SSB (t) ≥ 0 [2]. In other words, {N t , t > 0} satisfies the following conditions [2]: (1) The probability of detecting one photon in a time interval ∆t is λ SSB (t)∆t when ∆t → 0 .
(2) The probability of detecting more than one photon in ∆t is 0 when ∆t → 0 .
The rate function λ SSB (t) align includes all the arriving photons from the X-ray pulsar and the background. It has the form: λ SSB (t) = β SSB + α SSB h(φ SSB (t)) (ph/s) (2) where α SSB ≥ 0 and β SSB ≥ 0 are respectively the known effective source and background arrival rates in unit of photons per second, and h is the normalized pulsar profile, which is unique to a particular pulsar, defines the characteristic shape of the rate function, and is a continuous and nonnegative real function, periodic with period 1 (h(φ + n) = h(φ) for ∀n ∈ Z) and with unit area 1 0 h(φ)dφ = 1 [2]. φ SSB (t) is the observed phase at the SSB with respect to the coordinate time t seen at the SSB. It is expressed as the Taylor series: where f , .
f and ..
f are respectively the source frequency and its first and second derivatives, and φ 0 the initial phase at a reference time t 0 .

Photon Arrival Rate Function at the Spacecraft
Let r(t) and v(t) respectively denote the spacecraft's coordinate and velocity in the SSB coordinate frame at the spacecraft proper time t, and n is the directional vector from the SSB to the pulsar. The photon arrival rate function at the spacecraft λ(t) can be written as: with τ(t) the offset of proper time a photon arrives at the spacecraft compared to the arrival SSB-based coordinate time of the same photon at the SSB. For Earth-orbiting spacecraft, assuming that the spacecraft proper time is chosen as the Terrestrial Time (TT), τ(t) can be expressed as [29,30]: where c is the speed of light, D is the distance between the Sun and the pulsar, b is the position of SSB with respect to the Sun, µ s is the gravitational constant of Sun, ν E/SSB (t) is the velocity vector of Earth with respect to the SSB and P represents the total periodic correction items. The first term on the right-hand side of Equation (5) is the geometric time delay between the spacecraft and the SSB along the line-of-sight to the pulsar, the second term is due to the effects of parallax, and the third term is the Sun's Shapiro delay. The last two terms are the Einstein delay, which accomplishes the conversion of spacecraft proper time to SSB-based coordinate time. It is noteworthy that the Sun's Shapiro delay and the Einstein delay in the Equation (5) is referenced to the ones employed by the Heasoft v6.11.1 software platform provided by NASA [31]. In practical navigation, letting r(t) and v(t) respectively denote the position and velocity of the spacecraft predicted by the orbit dynamics and letting ∆r(t) and ∆ν(t) respectively denote the corresponding position and velocity errors, the accurate position and velocity of spacecraft can be expressed as: Substituting Equations (6)-(7) into Equation (4) and regarding the pulsar frequency as invariant within a time duration of n·∆r(t) c , we can expand the estimated phase model at the spacecraft φ SSB (t + τ(t)) as [20]: in which: is the pulsar phase model at the spacecraft derived by the spacecraft state with errors, p( r(t) + ∆r(t)), s( r(t) + ∆r(t)) and e( v(t) + ∆ν(t)) respectively represent the parallax delay, Sun's Shapiro delay and Einstein delay. The effects of spacecraft state error on them are extremely small and therefore they have been neglected.

Dynamic and Measurement Models of Proposed Navigation Method
The spacecraft's position x(t) and velocity ε(t) in the J2000.0 Earth Centered Inertial (ECI) coordinate system are selected as the state vector to be estimated. Letting [t s , t e ] denote the time interval of a filtering period, the navigation measurements are the estimated phaseφ(t e ) and Doppler φ(t e ) at the time t e as seen at the spacecraft, and the difference ofφ(t e ) and the estimated phase of the previous filtering period,φ(t s ). In what follows, the Dynamic model and measurement models of the proposed navigation method are formulated in detail.

Dynamics Model
In the Earth-centered inertial coordinate system, the dynamic model of earth-orbit spacecraft can be expressed as: where the process noise which can be characterized as a zero-mean Gaussian white noise process, and the acceleration of spacecraft a can be written as: considering the two-body effect of Earth a earth , the third-body effects from the Sun and the Moon, respectively denoted by a sun and a moon , and the Earth non-spherical perturbation acceleration a J2 [10].

PP and DF Measurements
Within a filtering period seen at the spacecraft, regarding ∆ν(t) and . φ SSB (t + n· r(t) c + p( r(t)) + s( r(t)) + e( v(t))) as invariant and equal to the one at the time t s , Equation (8) can be further simplified as: where ∆p = f s and: φ SSB (t s + n · r(t s ) c + p( r(t s )) + s( r(t s )) + e( v(t s ))). (13) ∆p and ∆ f are respectively the pulse phase delay and Doppler frequency offset which can be estimated from the detected photon event timestamps. The PP measurement equation is given as: n·(x(t e )+r E/SSB (t e )) c + p(x(t e ) + r E/SSB (t e )) +s(x(t e ) + r E/SSB (t e )) + e(ε(t e ) + ν E/SSB (t e ))) + u p (t e ) , with r E/SSB (t) being the Earth position vector relative to the SSB, T obs = t e − t s the integration interval of navigation filter and u p (t e ) the phase noise.φ(t e ) = φ(t e ) + ∆p + ∆f · T obs is the PP measurement, where ∆p and ∆f are respectively the estimated pulse phase delay and Doppler frequency offset. The DF measurement equation is obtained by differentiating Equation (14) and is written as: .
φ(t e ) + ∆f = f s (1 + n · ν(t e )/c) + ∆f is the DF measurement, in which the differentials of p( r(t)), s( r(t)) and e( v(t)) are neglected since they are extremely small compared to 1 and n · ν(t e )/c.

DPP Measurement
Taking the pulsar directional error, pulsar distance error and spacecraft-borne clock error into consideration, the PP measurement equation can be rewritten as: where n is the measured pulsar direction vector, D is the measured distance between the Sun and the pulsar, ∆n the pulsar direction error, ∆D the pulsar distance error, ∆t the spacecraft-borne clock error. E(∆n), E(∆D) and E(∆t) are the corresponding systematic biases caused by ∆n, ∆D and ∆t, respectively. E(∆n) and E(∆D) can respectively be approximated by the first order Taylor series: and: The derivation of E(∆n) satisfies: where v max is the spacecraft's maximum velocity in the SSB coordinate frame. Taking the PSR B1821-24 for example and assuming the maximum velocity of the spacecraft v max is 30 km/s, . E(∆n) is smaller than 2 × 10 −9 , meaning that the increment of E(∆D) is smaller than 7 × 10 −6 within 1 h, which is obviously ignorable. The derivation of E(∆D) satisfies: It can be known based on the JPL DE405 that the maximums of r E/SSB (t e ), v E/SSB (t e ) and b are respectively 1.5 × 10 8 km, 30 km/s, and 5 × 10 5 km. In addition, based on the study of [32], the minimum and uncertainty of D are respectively assumed as 2 kpc and 30%. Therefore, for PSR B1821-24, . E(∆D) ≤ 1 × 10 −10 , suggesting that the increment of E(∆D) within two adjacent navigation periods can be neglected. Furthermore, it has been studied that the spacecraft-borne clock error increases approximately by 36 ns/h [9], thus the variation of E(∆t) = f s ∆t can also be neglected.
All the above analysis suggests that the DPP measurement can reduce the major part of E(∆n) E(∆D) and E(∆t). Following from Equation (14), the DPP measurement equation can be written as: It is noteworthy that since the DPP measurement can only provide information about the relative position, the PP measurement is necessary to acquire information about the absolute position.

Measurements Estimation and Fusion Filtering
In Section 3.2.1, we have modeled the arrival rate function at the spacecraft with respect to the spacecraft's foregone state information and state error, seen as Equation (12). It turns out that the derivation of the rate function can be parameterized by a pulse phase delay ∆p and a Doppler frequency offset ∆ f , which are respectively caused by the initial position error and velocity error of the spacecraft and are essential for constructing the navigation measurements. In this section, based on the signal models and navigation system models given in the previous sections, the ML estimator used to estimate ∆p and ∆ f directly from the detected X-ray photon timestamps are presented, and then a modified unscented Kalman filter used to fuse the different measurements to update the spacecraft state are detailed.

ML Estimator
Based on the characteristics of NHPP, the joint pdf of the arrival timestamps {t i } K i=1 during a filtering period [t s , t e ] is derived as [2]: Utilizing Equation (22) as likelihood function, the unknown parameter (∆p, ∆ f ) can be estimated by: For the second term on the right hand side of Equation (23), since the term (∆v(t s ) · n)/c dependent of the unknown parameter ∆ f is much smaller than the other two terms independent of (∆p, ∆ f ), the overall second term can be approximately seen as independent of (∆p, ∆ f ). Furthermore, since the third term t e t s λ SC (s)ds is also independent of (∆p, ∆ f ), we end up with the following estimator: The 3σ search ranges of ∆p and ∆ f are respectively: in which σ( r(t s )) and σ( v(t s )) are the root mean squared error of r(t s ) and v(t s ) respectively. When the observation time is sufficiently long, the estimation accuracy of the presented ML estimator can achieve the Cramér-Rao lower bound (CRLB):  [20]. After ∆p and ∆ f have been estimated from the observed photon timestamps, the navigation measurements PP, DF and DFF can be obtained according to Equation (14), Equation (15) and Equation (21). Figure 1 shows the overall procedure of extracting navigation measurements from the observed photon timestamps.

Modified Unscented Kalman Filter
Since the orbit dynamic model and the measurement model of the proposed navigation method are both nonlinear, the Kalman filter cannot be used. In the XPNAV, the filtering period is usually several minutes. Under this condition, the EKF, which represents the nonlinear models by a first order Taylor series expansion and neglects the second and higher order terms, cannot perform very well. Since the UKF does not linearize the nonlinear equations, it does not cause the linearization error. For this reason, the UKF is a good choice for the XPNAV [24]. In this section, a modified unscented Kalman filter is deduced to upstate the spacecraft state. In order to reduce the timecorrelation of the DPP measurement equation, the previous state which is required in the DPP measurement equation is replaced by the sum of its estimate and the corresponding estimation error, and the mean square error matrix of the DPP measurement noise after the above transformation is derived.
The orbit dynamic model of the proposed navigation system can be rewritten as the following matrix form: The measurement equations of PP and DF can be expressed in the matrix form: (31) respectively, in which where k D is the DPP measurement vector:

Modified Unscented Kalman Filter
Since the orbit dynamic model and the measurement model of the proposed navigation method are both nonlinear, the Kalman filter cannot be used. In the XPNAV, the filtering period is usually several minutes. Under this condition, the EKF, which represents the nonlinear models by a first order Taylor series expansion and neglects the second and higher order terms, cannot perform very well. Since the UKF does not linearize the nonlinear equations, it does not cause the linearization error. For this reason, the UKF is a good choice for the XPNAV [24]. In this section, a modified unscented Kalman filter is deduced to upstate the spacecraft state. In order to reduce the time-correlation of the DPP measurement equation, the previous state which is required in the DPP measurement equation is replaced by the sum of its estimate and the corresponding estimation error, and the mean square error matrix of the DPP measurement noise after the above transformation is derived.
The orbit dynamic model of the proposed navigation system can be rewritten as the following matrix form: The measurement equations of PP and DF can be expressed in the matrix form: respectively, in which P k and F k are respectively the PP and DF measurement vectors, and V k and U k are the noise vectors, satisfying Formulating state of the previous epoch X k−1 as the sum of its estimateX k−1 and the corresponding estimation error ∆X k−1 , the DPP measurement equation can be rewritten as: where D k is the DPP measurement vector: and It can be deduced that: where K k−1 is the gain matrix and Q k = E V k V T k . After the above transformations, we can perform the procedures of standard unscented Kalman filter to update the spacecraft's position and velocity by recognizing Equation (10) as the orbit dynamic model and recognizing Equation (30), Equation (31) and Equation (32) as the measurement equations.

Photon-Level Simulations and Discussion
In this section, some simulations are carried out to evaluate the performance of the proposed navigation method, using photon event timestamps generated by the simulation method given in [33].
The pulsars PSR B1509-58, PSR B0540-69, PSR B0531+21 and PSR B0833-45 are selected as the navigation pulsars and their parameters are listed in Table 1. Figure 2 shows the normalized profiles of the selected four pulsars [20]. The initial orbit state is listed in Table 2. The simulation start time is chosen as MJD 52557.1155893. The Earth's position and velocity is predicted by the ephemeris JPL DE405. The detector area is 1 m 2 and the filtering period is 120 s. The background radiation flux for PSP B0531+21 is 1.015 ph/cm 2 /s [4]. For the other three pulsars, the background flux is assumed as 0.005 ph/cm 2 /s based on the research of the Naval Research Laboratory. The uncertainty of pulsar distance is set as 30% [32]. The initial clock error and the drift rate of clock frequency are respectively set as 1 µs and 10 −11 Hz [9].  Considering that the filtering period is longer than the observation time required for the selected four pulsars to achieve the Cramér-Rao lower bound (CRLB) and to account for other un-modeled error sources, the noise levels of ∆p and ∆f are set as two times the square-root of CRLB. The covariance matrix of the process noise is diag[q 2 1 , q 2 1 , q 2 1 , q 2 2 , q 2 2 , q 2 2 ], where q 1 = 2 × 10 −5 m and q 2 = 6 × 10 −4 m/s [10]. The measurement update period is set as 120 s. The initial position and velocity errors of spacecraft are respectively [100 km, 100 km, 100 km] and [100 m/s, 100 m/s, 200 m/s] and the covariance of the initial state is diag[q 2 3 , q 2 3 , q 2 3 , q 2 4 , q 2 4 , q 2 5 ], where q 3 = 100 km, q 4 = 100 m/s and q 5 = 200 m/s. error sources, the noise levels of p  andf  are set as two times the square-root of CRLB. The covariance matrix of the process noise is diag First, we compare our proposed method with the traditional method which only uses the PP measurement and the method which uses both the PP and DF measurements. For clarity, Table 3 gives the concepts of different measurements. When performing real-time photon level navigation simulations, the phase of the detected photons must be computed considering the dynamic motion of the spacecraft, otherwise the un-modeled Doppler effects caused by the motion of the spacecraft would result in the distortion of the folded pulse profiles of X-ray pulsars and severely degrade the pulse phase estimation accuracy which directly determines the navigation accuracy. Therefore, focusing on testifying the advantages of incorporating the DF and DPP measurements, in our experiments, we also use Equation (8) given in our article to compute the estimated phase of the detected photon timestamps for the traditional method which only uses the PP measurement. To focus on this purpose, we also use the UKF to update the spacecraft state for the latter two methods. Figure 3 shows the performance comparison of the three methods without systematic biases, and Figure 4 shows the results with systematic biases. For clarity, Table 4 lists the steady-state position and velocity estimation errors corresponding to the results shown in Figure 3 and Figure 4. It can be seen that no matter whether there are systematic biases, the method using the PP and DF measurements both obviously outperform the conventional XPNAV method which only uses the PP measurement, illustrating that the incorporation of the DF measurement in the navigation filter can yield better navigation accuracy. It is because that the modelling and estimation of the DF can eliminate the effect of the spacecraft orbital motion on the PP estimation accuracy to a large extent. First, we compare our proposed method with the traditional method which only uses the PP measurement and the method which uses both the PP and DF measurements. For clarity, Table 3 gives the concepts of different measurements. When performing real-time photon level navigation simulations, the phase of the detected photons must be computed considering the dynamic motion of the spacecraft, otherwise the un-modeled Doppler effects caused by the motion of the spacecraft would result in the distortion of the folded pulse profiles of X-ray pulsars and severely degrade the pulse phase estimation accuracy which directly determines the navigation accuracy. Therefore, focusing on testifying the advantages of incorporating the DF and DPP measurements, in our experiments, we also use Equation (8) given in our article to compute the estimated phase of the detected photon timestamps for the traditional method which only uses the PP measurement. To focus on this purpose, we also use the UKF to update the spacecraft state for the latter two methods. φ(t e ) DPP difference ofφ(t e ) andφ(t s ) Figure 3 shows the performance comparison of the three methods without systematic biases, and Figure 4 shows the results with systematic biases. For clarity, Table 4 lists the steady-state position and velocity estimation errors corresponding to the results shown in Figures 3 and 4. It can be seen that no matter whether there are systematic biases, the method using the PP and DF measurements both obviously outperform the conventional XPNAV method which only uses the PP measurement, illustrating that the incorporation of the DF measurement in the navigation filter can yield better navigation accuracy. It is because that the modelling and estimation of the DF can eliminate the effect of the spacecraft orbital motion on the PP estimation accuracy to a large extent.   Figure 5 shows the performance comparison between the proposed method, the method with augmenting state proposed in Ref. [24] and the conventional method which only uses the PP measurement. It can be seen that compared with the conventional method, the other two methods both can reduce the impact of systematic biases. Besides, the proposed method can achieve a more accurate steady state compared to the method with augmenting state. We attribute this to the reason that when focusing on compensating the composite impacts of various systematic biases, our   Figure 5 shows the performance comparison between the proposed method, the method with augmenting state proposed in Ref. [24] and the conventional method which only uses the PP measurement. It can be seen that compared with the conventional method, the other two methods both can reduce the impact of systematic biases. Besides, the proposed method can achieve a more accurate steady state compared to the method with augmenting state. We attribute this to the reason that when focusing on compensating the composite impacts of various systematic biases, our  It can also be observed that our proposed method which incorporates the DPP measurement as well as the PP and DF measurements obviously outperforms the method which only uses the PP measurement no matter whether there are systematic biases. Furthermore, when there are systematic biases, the navigation accuracy of our proposed method is also obviously higher than the one of method using the PP and DF measurements, which verifies that the incorporation of DPP measurement can reduce the effect of the systematic biases on navigation accuracy. However, when there are no systematic biases, it is observed that the performance of our proposed method is slightly worse than the method which uses the PP and DF measurements. We attribute this to the reason that the DPP measurement introduces more noise than useful information in this situation. Figure 5 shows the performance comparison between the proposed method, the method with augmenting state proposed in Ref. [24] and the conventional method which only uses the PP measurement. It can be seen that compared with the conventional method, the other two methods both can reduce the impact of systematic biases. Besides, the proposed method can achieve a more accurate steady state compared to the method with augmenting state. We attribute this to the reason that when focusing on compensating the composite impacts of various systematic biases, our proposed method has a better robustness to measurement noise compared to the method with augmenting state.  The performance of the proposed method is also studied under different initial velocity errors. Table 5 lists the steady-state position estimation errors of our proposed method and the method using PP and DPP measurements under different initial velocity error, with the other parameters being set as same as described above. It can be seen that the position error of the method using only PP and DPP measurements rapidly grows as the initial velocity error increases; on the other hand, position error of our proposed method only slightly grows as the initial velocity error increases and the advantage of our proposed method compared to the method without DF measurements becomes more and more notable. Finally, we investigate the effectiveness of the proposed method under different levels of systematic biases. Due to space limitations, here we only depict the impacts of pulsar directional errors on navigation performance. The impacts of pulsar distance error and spacecraft-borne clock error on navigation performance show the same regularity as that of pulsar directional errors. Figure  6 shows the performance comparison of our proposed method and the method using the PP and DF measurements under different angular position errors. The declination errors and right ascension errors are respectively set as n   and n   , of which  and  are the actual values of the angular position errors listed in Table 1 and n is the multiplication coefficient. We can see that as the pulsar The performance of the proposed method is also studied under different initial velocity errors. Table 5 lists the steady-state position estimation errors of our proposed method and the method using PP and DPP measurements under different initial velocity error, with the other parameters being set as same as described above. It can be seen that the position error of the method using only PP and DPP measurements rapidly grows as the initial velocity error increases; on the other hand, position error of our proposed method only slightly grows as the initial velocity error increases and the advantage of our proposed method compared to the method without DF measurements becomes more and more notable. Finally, we investigate the effectiveness of the proposed method under different levels of systematic biases. Due to space limitations, here we only depict the impacts of pulsar directional errors on navigation performance. The impacts of pulsar distance error and spacecraft-borne clock error on navigation performance show the same regularity as that of pulsar directional errors. Figure 6 shows the performance comparison of our proposed method and the method using the PP and DF measurements under different angular position errors. The declination errors and right ascension errors are respectively set as n * δθ and n * δϕ, of which δθ and δϕ are the actual values of the angular position errors listed in Table 1 and n is the multiplication coefficient. We can see that as the pulsar angular position error increases, our proposed method consistently outperforms the method which uses the PP and DF measurements and the advantage becomes more and more notable. It is noteworthy that the position estimation errors of the proposed method also increase as the pulsar angular position error increases. We expect this phenomenon since though the time differential method can eliminate the major part impact of the systematic biases on the DPP measurement, the PP and DF measurements are still affected by the systematic biases, which can lead to a decrease in the navigation performance.

PP PP + DF PP + DF + DPP
uses the PP and DF measurements and the advantage becomes more and more notable. It is noteworthy that the position estimation errors of the proposed method also increase as the pulsar angular position error increases. We expect this phenomenon since though the time differential method can eliminate the major part impact of the systematic biases on the DPP measurement, the PP and DF measurements are still affected by the systematic biases, which can lead to a decrease in the navigation performance.

Conclusion
This paper develops an innovative navigation method to eliminate the impact of Doppler effect caused by the spacecraft orbital motion and the systematic biases introduced by the pulsar directional error, spacecraft-borne clock error, etc. The method simultaneously employs the PP, the DF and the DPP of X-ray pulsars as navigation measurements. The PP and DF relative to the spacecraft's state estimation error is estimated by using the joint probability density function of the arrival photon timestamps as the likelihood function. Since the major part of systematic biases introduced by pulsar directional error, pulsar distance error and spacecraft-borne clock error is eliminated in the DPP, the DPP is also employed as navigation measurement to weaken the impact of systematic biases. Results of photon-level simulations testify to the advantages of the proposed method and show that simultaneously incorporating the PP, the DF and the DPP measurements can yield much better navigation accuracy than the method only using PP, the method using both PP and DF and the method using both PP and DPP for Earth orbit, and its advantage becomes more and more notable when the systematic biases increase.
It is noteworthy that the correlation between the PP and DPP measurements will inappropriately bias the filter, which is a shortcoming of the proposed method. If an effective method can be found to eliminate the correlation between the PP and DPP measurements, the navigation accuracy can be further improved, which is our next work. The batch least squares method can also be used to provide the DPP measurement, which is also one of our next work. APPENDIX Important parameters' notation and meaning.

Conclusions
This paper develops an innovative navigation method to eliminate the impact of Doppler effect caused by the spacecraft orbital motion and the systematic biases introduced by the pulsar directional error, spacecraft-borne clock error, etc. The method simultaneously employs the PP, the DF and the DPP of X-ray pulsars as navigation measurements. The PP and DF relative to the spacecraft's state estimation error is estimated by using the joint probability density function of the arrival photon timestamps as the likelihood function. Since the major part of systematic biases introduced by pulsar directional error, pulsar distance error and spacecraft-borne clock error is eliminated in the DPP, the DPP is also employed as navigation measurement to weaken the impact of systematic biases. Results of photon-level simulations testify to the advantages of the proposed method and show that simultaneously incorporating the PP, the DF and the DPP measurements can yield much better navigation accuracy than the method only using PP, the method using both PP and DF and the method using both PP and DPP for Earth orbit, and its advantage becomes more and more notable when the systematic biases increase.
It is noteworthy that the correlation between the PP and DPP measurements will inappropriately bias the filter, which is a shortcoming of the proposed method. If an effective method can be found to eliminate the correlation between the PP and DPP measurements, the navigation accuracy can be further improved, which is our next work. The batch least squares method can also be used to provide the DPP measurement, which is also one of our next work. E(∆t) corresponding systematic biases caused by ∆t v max spacecraft's maximum velocity in the SSB coordinate frame ∆φ(t e ) the difference ofφ(t e ) and the estimated phase of the previous filtering period,φ(t s ) P PP measurement vector of the navigation filter F DF measurement vector of the navigation filter V PP noise vector of the navigation filter U DF noise vector of the navigation filter D DPP measurement vector of the navigation filter Γ DPP noise vector of the navigation filter