X-ray Pulsar-Based Navigation Using Pulse Phase Delay between Spacecraft and Verification with Real Data

: Pulsars are neutron stars with high rotation speeds and have extraordinary long-term rotational stability. X-ray pulsar-based navigation (XNAV) is a navigation method that estimates the position and velocity of a spacecraft using the X-ray radiation from pulsars. Flight experiments on Insight-Hard X-ray Modulation Telescope (Insight-HXMT) and Neutron Star Interior Composition Explorer (NICER) have successfully verified the feasibility of using XNAV for a single spacecraft. For spacecraft in formation, a pulsar-based navigation method that uses the pulse phase delay between spacecraft is derived. Moreover, a direct estimation method for pulse phase delay, which is independent from the pulsar template, is proposed. The proposed method is verified with simulation data of the Crab pulsar and real data of the same pulsar obtained from Insight-HXMT and NICER.


Introduction
With the development of aerospace technology, on-orbit spacecraft missions have become more complex [1].Formation spacecraft have the advantages of low system cost and high flexibility [2,3].Compared to a single spacecraft, spacecraft formation flying is more flexible for complex space missions [4].Thus, spacecraft formation flying has attracted a lot of attention [5,6].Furthermore, navigation services are important for spacecraft formation flying missions [7].
XNAV is a developing autonomous navigation method providing navigation services for Earth-orbit and deep-space spacecraft [8,9].The basic idea of the XNAV for a single spacecraft is similar with that of the Global Navigation Satellite System.A spacecraft estimates its position and velocity using received photons from different pulsars [10].Since its introduction in 1981, a significant growth in the navigation algorithm for XNAV has occurred in the past four decades [11,12].In addition, some flight experiments, such as NICER and Insight-HXMT, were performed [13][14][15][16].
Aside from providing navigation service for a single spacecraft, XNAV can also be used for spacecraft in formation [17][18][19].In [17], Xiong proposed the use of time difference of arrival (TDOA) measurements from pulsars and pseudo-range observations from intersatellite links to determine the position of a satellite in the constellation.In [18], Zhang proposed using TDOA, phase increment, the baseline of double satellites, and the plane formed by three satellites for spacecraft formation navigation.In [19], Emadzadeh proposed the use of pulse delay between two spacecraft to estimate the relative position and studied the effect of velocity errors on position estimation.However, in [19], the pulse delay is calculated after pulse phases of spacecraft are estimated, respectively.This can lead to unnecessary redundant calculations.Moreover, in [19], the pulse phase of spacecraft is estimated using the template of pulsar and empirical profiles recovered form pulsar photons.However, this can affect the estimation performance of the pulse phase based on the accuracy of the template.In [17][18][19], the proposed methods were only verified through simulation.
In this paper, we first described a pulse phase estimation method and the measurement model of XNAV for a single spacecraft.Then, taking the formation of two spacecraft as an example, a pulsar-based navigation method is used to derive the pulse phase delay between spacecraft.To improve computational efficiency and eliminate dependence on the template of pulsars, a direct estimation method for pulse delay without using a template is proposed.The proposed method is verified with simulation data of the Crab pulsar and real data of the Crab pulsar derived from NICER and Insight-HXMT.For convenience, we only consider the Romer delay and Shapiro delay when determining the transform photon arrival time from spacecraft to solar system barycenter (SSB).To provide an accurate navigation service for spacecraft near gravitational sources, an accurate transform model should be considered in further studies [20].
The remainder of the paper is organized as follows: The related method of XNAV for a single spacecraft is shown in Section 2; Section 3 derives a pulsar-based navigation method using the pulse phase delay between spacecraft, and a direct estimation method is proposed for pulse phase delay.In Section 4, the simulation data of the Crab pulsar and the real data of the Crab pulsar derived from NICER and Insight-HXMT are used to verify the proposed method.The conclusions are presented in Section 5.

Measurement Model
The measurement of pulsar-based navigation is the pulse phase.The pulse phase at t i can be expressed as [16] follows: where F 0 and F 1 are the frequency of the pulsar and the time derivative of the frequency at T 0 , respectively, and F 0 and F 1 can be found in the public pulsar ephemeris [21].In addition, g(t i ) can be described by [20,22] the following: where c is the speed of light, r i is position vector of the spacecraft relative to the Earth at t i , n is the direction vector of the pulsar, r E,i is the position vector of Earth relative to the SSB at t i , r k,i is the position of the spacecraft relative to the kth celestial body(unit m), µ k is the gravitational coefficient of the kth celestial body, and ∥ ∥ is the magnitude of the vector.The second and third terms on the right side of Equation ( 2) are the Romer delay and the Shapiro delay, respectively.The Romer delay represents a vacuum delay between the spacecraft and the SSB, while the Shapiro delay represents a time delay caused by the pulse passing through curved space-time.
Assuming that is the velocity of the spacecraft relative to the Earth at t i .Thus, Equation (1) can be rewritten as [16] follows:

Pulse Phase Estimation Method
The pulsar is very weak such that a spacecraft can only receive a series of events rather than a continuous pulsed signal.Thus, a spacecraft must estimate the pulse phase from received events [23].
Let us assume a spacecraft observes a pulsar from t i to t i+1 , i = 1, 2, . ..I and receives N events τ j } N j=1 .At first, τ j } N j=1 should be corrected to SSB by (2); thus, we have 16: Therefore, the pulse phase at τ j,SSB is expressed as .
where ϕ i = ϕ(t i ) is the pulse phase at g(t i ), F i is the frequency at g(t i ), and .
F i is the time derivative of the frequency at g(t i ), respectively.
However, the true position of the spacecraft, r τ j in (2), is inaccessible in an autonomous navigation task.r τ j is expressed as r τ j = ∼ r τ j + δr τ j to estimate the pulse phase, where ∼ r τ j is the predicted position of the spacecraft and δr τ j is the error of ∼ r τ j .Moreover, δr τ j can be derived as functions of δv(τ i ) and δr(τ i ) [24-26]: where δv(τ i ) is the error of the predicted velocity of the spacecraft ∼ v τ j .Inaddition, Φ rr τ j , τ i and Φ rv τ j , τ i can be expressed as [24,25] follows: where φ m and γ m are constant matrixes that can be calculated using analytical, numerical, or finite difference techniques [26].Therefore, (5) can be rewritten as [27]: where In the equations above, ∼ ϕ τ j,SSB and ∼ g(t i ) are the pulse phase and photon arrival time calculated by ∼ r (t i ), respectively.Applying (8) to every event τ j } N j=1 , we have ϕ τ j,SSB } N j=1 .Thus, we have to estimate ϕ i 0 and F orbit i from ϕ τ j,SSB } N j=1 .Here, F orbit i and ϕ i 0 can be estimated using the epoch folding method.The observation period is assumed to consist of M pulsar periods in the epoch folding method.The period can be divided into M b equal-length bins.The length of each bin is assumed as t b , and c m l is the number of detected photons in the l − th bin of the m − th period.Therefore, the empirical profile can be recovered.Here, we define the empirical profile function at φ l as [19] the following: where φ l is the pulse phase at the middle of l − th bin.We adopt the Chi-squared test to calculate the statistical characteristics of the empirical profile as follows [28]: where Using Equations ( 14)-( 16), υ 1 can be estimated with Inaddition, ϕ i can be estimated by solving a nonlinear least squares estimation problem [19]: where λ Te is the template of the pulsar.At the same time, ϕ i can be estimated with φi = argminJ(ϕ i )

Measurement Model
The diagram of XNAV for spacecraft in formation is shown in Figure 1.Assume that spacecraft A and spacecraft B observe the same pulsar; thus, according to (1), the pulse phases are given by where where r A i is the position of spacecraft A relative to the Earth at t i andr B i is the position of spacecraft B relative to the Earth at t i .Given that the distance between the two spacecraft is far less than the distances between the celestial bodies and the spacecraft, we assume that r A k,i ≈ r B k,i , and we denote it as r k,i .We use pulse phase delay between spacecraft A and spacecraft B as a measurement: Substituting ( 20) and ( 21) into ( 24), we have ) ) Since c is about 3 × 10 8 m/s,  is smaller than 10 −10 HZ/s for most pulsars.Assume that the distance between spacecraft A and B is shorter than 100 km.Therefore, is close to 0, and we have We also define Substituting ( 27) into (25), we have Based on Equation ( 28), the difference in pulse phases of spacecraft A and spacecraft B reflects their distance in the direction of the pulsar.Assuming that Δ , =  ( ) −  ( ), , = Δ , and ℎ , (Δ , ) =  , ⋅(  ) , then Equation ( 28) can be rewritten as a standard measurement model of  , = ℎ , ( , ). (29)

Pulse Phase Delay Estimation
In this section, we explain how the pulse phase delay can be abstained.We use pulse phase delay between spacecraft A and spacecraft B as a measurement: Substituting ( 20) and ( 21) into ( 24), we have Since c is about 3 × 10 8 m/s, F 1 is smaller than 10 −10 HZ/s for most pulsars.Assume that the distance between spacecraft A and B is shorter than 100 km.Therefore, is close to 0, and we have We also define F A,B i as follows: Substituting ( 27) into (25), we have Based on Equation ( 28), the difference in pulse phases of spacecraft A and spacecraft B reflects their distance in the direction of the pulsar.Assuming that ∆x A,B i , then Equation ( 28) can be rewritten as a standard measurement model of

Pulse Phase Delay Estimation
In this section, we explain how the pulse phase delay can be abstained.Assume that spacecraft A and spacecraft B observe the same pulsar from t i to t i+1 , i = 1, 2,. ..I and receives τ A j } N A j=1 and τ B j } N B j=1 , respectively.Applying (8) to τ A j } N A j=1 and τ B j } N B j=1 , we have ϕ τ A j,SSB } N j=1 and ϕ τ B j,SSB } N j=1 , where As shown in Section 2.2, F orbit,A i and F orbit,B i can be estimated by (17).In Equation ( 19), ϕ A i and ϕ B i can be estimated by (18), and ∆ϕ A,B i can be calculated by (24).However, when estimating ϕ A i and ϕ B i , respectively, we have to solve the optimization problem shown in (18) twice, and the template is indispensable.To improve computational efficiency and eliminate dependence on the template, we proposed a direct estimation method for ∆ϕ A,B i .Similar to (19), we define the cost function for estimating ∆ϕ A,B i as follows: Therefore, ∆ϕ A,B i can be estimated by As shown in (32) and (33), ∆ϕ A,B i can be estimated by solving (33) once and is independent from the template of pulsars.

Procedure of Pulsar-Based Navigation Using Pulse Phase Delay
In this section, we discuss the procedure of pulsar-based navigation using pulse phase delay between two spacecraft.We assume that the true position and velocity of spacecraft A is known, and the position and velocity of spacecraft B can be estimated by measuring pulse phase delay between spacecraft A and B. By measuring ∆ϕ A,B i , we can obtain the distance of the direction of the pulsar between spacecraft A and spacecraft B. Using one measurement, we can determine a ring of coordinate triples for spacecraft B.Moreover, given that the dynamic model of spacecraft B can be modeled in advance, we use an unscented Kalman filter (UKF) to estimate the state of spacecraft B using the measurement and dynamic model.The procedure of pulsar-based navigation using pulse phase delay is shown in Algorithm 1.At first, we apply (30) and (31) to τ A j } N A j=1 and τ B j } N B j=1 and obtain ϕ τ A j,SSB } N j=1 and ϕ τ B j,SSB } N j=1 , respectively.Then, we estimate Forbit,A i , Forbit,B i by (17) and estimate ∆ϕ A,B i by (33).The state of spacecraft B is then estimated using UKF.
Algorithm 1. Procedure of pulsar-based navigation using pulse phase delay Predict the state of spacecraft B at {τ B j } N B j=1 by propagating ∼ x B (t i ).

Experiments
In this section, we used the simulation data of Crab pulsar and real data of Crab pulsar derived from NICER and Insight-HXMT to verify the proposed pulsar-based navigation method using pulse phase delay.

Experiments with Simulation Data
Here, the Crab pulsar with simulation parameters shown in Table 1 was taken as the observed pulsar.Figure 2 shows the template of Crab pulsar.The initial orbit elements of spacecraft A and spacecraft B are shown in Table 2.Both spacecraft were assumed to observe Crab pulsar with observation t obs = 100 s.In this section, we assumed that the true position and velocity of spacecraft A was known, and the observation data of Crab pulsar were used to estimate state of spacecraft B with Algorithm 1.The initial state errors of spacecraft B were set as [5 km, 5 km, 5 km] and [5 m/s, 5 m/s, 5 m/s].Here, the Crab pulsar with simulation parameters shown in Table 1 was taken as the observed pulsar.Figure 2 shows the template of Crab pulsar.The initial orbit elements of spacecraft A and spacecraft B are shown in Table 2.Both spacecraft were assumed to observe Crab pulsar with observation  = 100 s.In this section, we assumed that the true position and velocity of spacecraft A was known, and the observation data of Crab pulsar were used to estimate state of spacecraft B with Algorithm 1.The initial state errors of spacecraft B were set as [5 km, 5 km, 5 km] and [5 m/s, 5 m/s, 5 m/s]., computed via (29).Therefore, the pulse phase delay can reflect the distance between spacecraft A and B.  The initial state of spacecraft A and B was used to recover empirical profiles with ϕ τ HXMT j,SSB and ϕ τ N ICER j,SSB via (13).The recovered empirical profiles are shown in Figure 3. Due to the distance between spacecraft A and B, pulse phase of Spacecraft A and B in the same exposure were different.Using (33) to estimate ∆ φA,B i , the value of estimated pulse phase delay was 0.2650, and ∆ϕ A,B i , calculated by ( 29), was 0.2652.It can be seen that ∆ φA,B i was close to ∆ϕ A,B i , computed via (29).Therefore, the pulse phase delay can reflect the distance between spacecraft A and B. Figures 4 and 5 show the navigation result of spacecraft B. As can be seen, the position and velocity gradually reduced to below 1 km and 1 m/s, respectively, thus verifying the feasibility of proposed navigation method using pulse phase delay.

Experiments with Real Data
We use real data of Crab pulsar from LE (Low Energy X-ray Telescope) on Insight-HXMT and NICER on 13 March 2018 and 14 March 2018 (Obs_ID: P0111605054 (Insight-HXMT), 1013010125 (NICER), and 1013010126 (NICER)).The start and finish time of real data are shown in Table 3. Figures 4 and 5 show the navigation result of spacecraft B. As can be seen, the position and velocity gradually reduced to below 1 km and 1 m/s, respectively, thus verifying the feasibility of proposed navigation method using pulse phase delay.Figures 4 and 5 show the navigation result of spacecraft B. As can be seen, the position and velocity gradually reduced to below 1 km and 1 m/s, respectively, thus verifying the feasibility of proposed navigation method using pulse phase delay.

Experiments with Real Data
We use real data of Crab pulsar from LE (Low Energy X-ray Telescope) on Insight-HXMT and NICER on 13 March 2018 and 14 March 2018 (Obs_ID: P0111605054 (Insight-HXMT), 1013010125 (NICER), and 1013010126 (NICER)).The start and finish time of real data are shown in Table 3.  Figures 4 and 5 show the navigation result of spacecraft B. As can be seen, the position and velocity gradually reduced to below 1 km and 1 m/s, respectively, thus verifying the feasibility of proposed navigation method using pulse phase delay.

Experiments with Real Data
We use real data of Crab pulsar from LE (Low Energy X-ray Telescope) on Insight-HXMT and NICER on 13 March 2018 and 14 March 2018 (Obs_ID: P0111605054 (Insight-HXMT), 1013010125 (NICER), and 1013010126 (NICER)).The start and finish time of real data are shown in Table 3.

Experiments with Real Data
We use real data of Crab pulsar from LE (Low Energy X-ray Telescope) on Insight-HXMT and NICER on 13 March 2018 and 14 March 2018 (Obs_ID: P0111605054 (Insight-HXMT), 1013010125 (NICER), and 1013010126 (NICER)).The start and finish time of real data are shown in Table 3.To verify the proposed pulsar-based navigation method using pulse phase delay, we chose the real data when Insight-HXMT and NICER simultaneously observe the Crab pulsar.The exposures on Crab pulsar from NICER and Insight-HXMT are shown in Figure 6.As shown in the figure, 13 March 2018 (MJD:58190) Insight-HXMT and NICER did not observe Crab pulsar simultaneously.On 14 March 2018 (MJD:58191), there are nine exposures that meet the criteria.We select these nine exposures and align the start time of each exposure.The start times of the selected exposures are shown in Table 4.To verify the proposed pulsar-based navigation method using pulse phase delay, we chose the real data when Insight-HXMT and NICER simultaneously observe the Crab pulsar.The exposures on Crab pulsar from NICER and Insight-HXMT are shown in Figure 6.As shown in the figure, 13 March 2018 (MJD:58190) Insight-HXMT and NICER did not observe Crab pulsar simultaneously.On 14 March 2018 (MJD:58191), there are nine exposures that meet the criteria.We select these nine exposures and align the start time of each exposure.The start times of the selected exposures are shown in Table 4.Given that there are only nine exposures in which Insight-HXMT and NICER simultaneously observe Crab pulsar, we verified the proposed pulsar-based navigation method by comparing pulse phase delay calculated by (28) and that estimated by (33).With the pulse phases of events received by Insight-HXMT and NICER, {( , )} and {( , )} were calculated by (20) and ( 21) using the real position of NICER and Insight-HXMT.The empirical profiles of NICER and Insight-HXMT in each exposure were then recovered using {( , )} and {( , )} via (13).The empirical profiles in each exposure are shown in Figure 7.Given the distance between NICER and Insight-HXMT in the direction of pulsar, the pulse phases of NICER and Insight-HXMT in the same exposure were different.In accordance with (24), the pulse phase delay Δ , reflects the distant between NICER and Insight-HXMT in the direction of the pulsar.Given that there are only nine exposures in which Insight-HXMT and NICER simultaneously observe Crab pulsar, we verified the proposed pulsar-based navigation method by comparing pulse phase delay calculated by (28) and that estimated by (33).With the pulse phases of events received by Insight-HXMT and NICER, ϕ τ HXMT j,SSB and ϕ τ N ICER j,SSB were calculated by (20) and ( 21) using the real position of NICER and Insight-HXMT.The empirical profiles of NICER and Insight-HXMT in each exposure were then recovered using ϕ τ HXMT j,SSB and ϕ τ N ICER j,SSB via (13).The empirical profiles in each exposure are shown in Figure 7.Given the distance between NICER and Insight-HXMT in the direction of pulsar, the pulse phases of NICER and Insight-HXMT in the same exposure were different.In accordance with (24), the pulse phase delay ∆ φHXMT,NICER i reflects the distant between NICER and Insight-HXMT in the direction of the pulsar.
Next, we used the estimation error of ∆ φHXMT,NICER i to reflect the measurement error of the distance between Insight-HXMT and NICER.The pulse phase delays ∆ φHXMT,NICER i , i = 1, . . ., 9 were estimated by (33) using the empirical profiles of NICER and Insight-HXMT shown in Figure 6.The real pulse phase delays were calculated by (28) using the real positions of NICER and Insight-HXMT.Furthermore, the estimation error of pulse phase delay was defined as the difference between estimated pulse phase delay and real pulse phase delay.Table 5 shows the estimation errors of the pulse phase delay of nine exposures.The results revealed that the estimated pulse phase delay can accurately reflect the distance between Insight-HXMT and NICER, thus verifying the feasibility of the proposed pulsar-based navigation using the pulse phase delay.Next, we used the estimation error of Δ , to reflect the measurement error of the distance between Insight-HXMT and NICER.The pulse phase delays Δ , ,  = 1, . . .,9 were estimated by (33) using the empirical profiles of NICER and Insight-HXMT shown in Figure 6.The real pulse phase delays were calculated by (28) using the real positions of NICER and Insight-HXMT.Furthermore, the estimation error of pulse phase delay was defined as the difference between estimated pulse phase delay and real pulse phase delay.Table 5 shows the estimation errors of the pulse phase delay of nine exposures.The results revealed that the estimated pulse phase delay can accurately reflect the distance between Insight-HXMT and NICER, thus verifying the feasibility of the proposed pulsar-based navigation using the pulse phase delay.Finally, we compared the performance of the separate estimation method for pulse phase delay used in (19) and the proposed direct estimation method.In the proposed direct estimation method, pulse phase delays ∆ φHXMT,NICER i , i = 1, . . ., 9 were directly estimated by (33) using the empirical profiles of Insight-HXMT and NICER in each exposure.In the respective estimation method, ∆ φHXMT i and ∆ φNICER i were estimated, respectively, by (19) using empirical profiles and templates, after which ∆ φHXMT,NICER i , i = 1, . . ., 9 were then calculated by (24).The evaluation criterion of the estimation performance is defined as RMS (root mean square) error of the pulse phase delay.We defined CPU time cost as evaluation criterion of computational cost.The computation environment contained Intel Core i5-1240P@1.7 GHZ (Intel, Santa Clara, CA, USA) and python 3.8.Table 6 shows the RMS error and mean CPU time cost of the proposed direct estimation method and the separate estimation method for pulse phase delay.The results showed that the estimation performance of the proposed direct estimation method was close to that of respective estimation method.Furthermore, the CPU time cost of the proposed direct estimation method was about 36.99% lower than that of respective estimation method.

Conclusions
This study derived a pulsar-based navigation method using pulse phase delay between spacecraft.A direct estimation method for pulse phase delay is proposed to eliminate reliance on templates of pulsars for pulse phase delay estimation.The proposed method is verified with the simulation data of the Crab pulsar and the real data of the Crab pulsar form Insight-HXMT and NICER.Real data experiment results show that the pulse phase delay estimated by the real data of the Crab pulsar can accurately reflect the distance between Insight-HXMT and NICER.Moreover, compared with the estimation method that estimates the pulse phase of two spacecraft, the proposed direct estimation method can estimate pulse phase delay between two spacecraft by running the pulse phase estimation algorithm just once.The experimental results show that the estimation performance of the proposed direct estimation method is close to that of the respective estimation method and that the CPU time cost of the proposed method is about 36.99% lower.Notably, the transform model shown as (2) only considers the Romer delay and Shapiro delay.An accurate transform model that considers the Einstein delay and high-order terms of relativistic effects should be considered in further studies to provide an accurate navigation service for spacecraft near gravitational sources.

Figure 1 .
Figure 1.Diagram of XNAV for spacecraft in formation.

Figure 1 .
Figure 1.Diagram of XNAV for spacecraft in formation.

Figure 2 .
Figure 2. Template of the Crab pulsar.

Figure 2 .
Figure 2. Template of the Crab pulsar.

Figure 3 .
Figure 3. Empirical profiles of spacecraft A and B.

Figure 4 .
Figure 4. Position estimation error of spacecraft B.

Figure 5 .
Figure 5. Velocity estimation error of spacecraft B.

Figure 3 .
Figure 3. Empirical profiles of spacecraft A and B.

Figure 3 .
Figure 3. Empirical profiles of spacecraft A and B.

Figure 4 .
Figure 4. Position estimation error of spacecraft B.

Figure 5 .
Figure 5. Velocity estimation error of spacecraft B.

Figure 4 .
Figure 4. Position estimation error of spacecraft B.

Figure 3 .
Figure 3. Empirical profiles of spacecraft A and B.

Figure 4 .
Figure 4. Position estimation error of spacecraft B.

Figure 5 .
Figure 5. Velocity estimation error of spacecraft B.

Figure 5 .
Figure 5. Velocity estimation error of spacecraft B.

Figure 7 .
Figure 7. Empirical profiles of each exposure.

Table 2 .
Orbital elements of spacecraft A and B.

Table 2 .
Orbital elements of spacecraft A and B.

Table 3 .
Start and finish time of real data.

Table 3 .
Start and finish time of real data.

Table 4 .
The start times of the selected nine exposures.

Table 4 .
The start times of the selected nine exposures.

Table 5 .
Estimation errors of the pulse phase delay.

Table 5 .
Estimation errors of the pulse phase delay.

Table 6 .
RMS error and mean CPU time cost for different pulse phase delay estimation methods.