2. Modified measurement model without considering the pulsar template

This paper studies a novel method of pulsar-based navigation, where the point is to avoid dependency on the template of the Crab pulsar. Schematics of three pulse profiles commonly used in dynamic signal processing are shown in Figure 1. The pulsar template is the basic input to the existing methods for calculating the pulse TOA, which usually needs two or three days of accumulated observation data to build. The modified observation profile in this paper is the pulse profile recovered after period search (Zhou et al., Reference Zhou, Ji and Ren2013), the Doppler effects from the motion of the vehicle are basically eliminated. The observation profile is that directly recovered with the prior information of pulsar frequency.

Figure 1. Schematic of observation profiles used in dynamic signal processing

For the existing theoretical framework of pulsar-based navigation, the pulse TOA can be calculated by the phase difference between the modified observation profile and the pulsar template. The general measurement model is

(1)\begin{equation}\Delta \textrm{TO}{\textrm{A}_1} = \frac{1}{c}\boldsymbol{n}\Delta \boldsymbol{r}\end{equation}

where $\boldsymbol{n}$ denotes the direction of the pulsar, c is the speed of light, and $\varDelta \boldsymbol{r}$ usually represents the distance between the investigated vehicle and SSB (solar system barycentre).

For the Crab pulsar, however, the aim here is to find a different method of TOA process without using the template. We can easily find that only the frequency parameters of the pulsar used for building the profile are different between the other two recovered profiles. So, there will remain a distort phase difference. As the recovered profiles are integrated with all pulses during the observation period, this distort phase difference should equal the mean phase propagation error over the entire period, which can be converted to a relative expression of velocity error of the vehicle at a given moment, in theory. This distort phase difference can also be regarded as a difference model between two kinds of phase propagation error based on the pulsar template with the advantage that most system errors can be ignored. There is also no need to consider the problem of ambiguity of whole cycles.

With reference to the beginning of the observation ${t_0}$, the detected phase propagation model at the onboard sensor is (Emadzadeh and Speyer, Reference Emadzadeh and Speyer2011a)

(2)\begin{equation}{\phi _{\det }}(t )= {\phi _0} + \int_{{t_0}}^t {f(\tau )} d\tau\end{equation}

where ${\phi _0}$ is the initial phase at ${t_0}$, $f(t )$ is the observed signal frequency, which can be decomposed into two different components: X-ray source frequency ${f_s}$ (the spinning frequency of pulsar), and the Doppler frequency shift ${f_d}(t )= {f_s}\frac{{v(t )}}{c}$.

Taking the orbiting motion of the vehicle into account, ${\phi _{\det }}(t )$ can be written as

(3)\begin{equation}{\phi _{\det }}(t )= {\phi _0} + {f_s}({t - {t_0}} )+ \frac{{{f_s}}}{c}\int_{{t_0}}^t {v(\tau )} d\tau\end{equation}

where v is the vehicle's velocity projected on the direction of the pulsar.

Given that the pulsar is far from the solar system, the velocity of the vehicle projected along the direction of the pulsar is

(4)\begin{equation}v = \boldsymbol{n} \cdot ({\dot{\boldsymbol{r}} - \dot{\boldsymbol{D}}} )\end{equation}

where $\boldsymbol{r}$ is the position of the vehicle relative to SSB, $\boldsymbol{D}$ is the position of the pulsar in the same reference system.

Substituting Equation (4) into Equation (3) yields

(5)\begin{equation}{\phi _{\det }}(t )= {\phi _0} + {f_s}({t - {t_0}} )+ \frac{{{f_s}}}{c}\int_{{t_0}}^t {\boldsymbol{n} \cdot \dot{\boldsymbol{r}}(\tau )d\tau } - \frac{{{f_s}}}{c}\int_{{t_0}}^t {\boldsymbol{n} \cdot \dot{\boldsymbol{D}}(\tau )d\tau }\end{equation}

Then, if the proper motion of the pulsar during the observation period can be neglected, that is,$\boldsymbol{D}(t )\approx \boldsymbol{D}({{t_0}} )$, the detected phase propagation model can be simplified as (Wang, Reference Wang2016)

(6)\begin{equation}{\phi _{\det }}(t) = {\phi _0} + {f_s}(t - {t_0}) + \frac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\boldsymbol{r}(t) - \boldsymbol{r}({t_0})} )\end{equation}

The true state of the vehicle at an arbitrary epoch t can be expressed as the term of predicted state and propagation error.

(7)\begin{equation}\left\{ {\begin{array}{*{20}{c}} {\boldsymbol{r}(t )= \tilde{\boldsymbol{r}}(t )+ \delta \boldsymbol{r}(t )}\\ {\boldsymbol{v}(t )= \tilde{\boldsymbol{v}}(t )+ \delta \boldsymbol{v}(t )} \end{array}} \right.\end{equation}

where $\tilde{\boldsymbol{r}}({\cdot} )$ and $\tilde{\boldsymbol{v}}({\cdot} )$ denote the predicted position and velocity, $\delta \boldsymbol{r}({\cdot} )$ and $\delta \boldsymbol{v}({\cdot} )$ denote the corresponding propagation errors.

Substituting Equation (7) into Equation (6) we have

(8)\begin{equation}{\phi _{\det }}(t) = {\phi _0} + {f_s}(t - {t_0}) + \frac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\tilde{\boldsymbol{r}}(t) - \tilde{\boldsymbol{r}}({t_0})} )+ \frac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\delta \boldsymbol{r}(t) - \delta \boldsymbol{r}({t_0})} )\end{equation}

To subtract $\frac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\tilde{\boldsymbol{r}}(t) - \tilde{\boldsymbol{r}}({t_0})} )$, which is called the orbital motion effect term, all the photon TOAs should be transformed to the position of the vehicle at ${t_0}$ with the aid of predicted state information. Equation (7) can then be rewritten as

(9)\begin{equation}{\phi _{\det }}(t) = {\phi _0} + {f_s}(t - {t_0}) + \frac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\delta \boldsymbol{r}(t) - \delta \boldsymbol{r}({t_0})} )\end{equation}

Assume that the observation period is ${T_{obs}} = T - {t_0}$ and comprising ${N_p}$ pulsar periods, namely, ${T_{obs}} = {N_p}P$. In the ith bin, the profile function after epoch folding can be represented as (Emadzadeh and Speyer, Reference Emadzadeh and Speyer2011b)

(10)\begin{align}\begin{aligned} \hat{\lambda }({t_i}) & =\dfrac{1}{{{N_p}}}\sum\limits_{i = 1}^{{N_p}} {\lambda \left( {{t_i};{\phi_0} + \dfrac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\delta \boldsymbol{r}(iP) - \delta \boldsymbol{r}({t_0})} )} \right)} + n({t_i})\\ & \approx \lambda ({{t_i};{\phi_0}} )+ \dfrac{1}{{{N_p}}}\sum\limits_{i = 1}^{{N_p}} {\lambda \left( {\dfrac{{{f_s}}}{c}\boldsymbol{n} \cdot ({\delta \boldsymbol{r}(iP) - \delta \boldsymbol{r}({t_0})} )} \right)} \end{aligned}\end{align}

where P denotes the period of the pulsar, $\lambda ({\cdot} )$ is the photon rate function of the pulsar signal, ${\lambda ^{\prime}}({\cdot} )$ expresses the derivative of $\lambda ({\cdot} )$ and ${t_i}$ denotes the middle time of the $i$th bin. $n({t_i})$ is referred to the epoch folding noise, which is uncorrelated over a long observation time.

Using the orbital transition matrix, $\delta \boldsymbol{r}(iP)$ can be expressed as a term of $\delta \boldsymbol{r}({t_0})$ and $\delta \boldsymbol{v}({t_0})$, i.e.,

(11)\begin{equation}\delta \boldsymbol{r}({iP} )= {\boldsymbol{\Phi }_{rr}}(iP,{t_0})\delta \boldsymbol{r}({{t_0}} )+ {\boldsymbol{\Phi }_{rv}}(iP,{t_0})\delta \boldsymbol{v}({{t_0}} )\end{equation}

where ${\boldsymbol{\Phi }_{rr}}(t,{t_0})$ and ${\boldsymbol{\Phi }_{rv}}(t,{t_0})$ are the partition matrices of orbital transition matrix (Zhang, Reference Zhang2015), and the concrete expressions can be found in the appendix.

Substituting Equation (11) into Equation (10) yields

(12)\begin{equation}\hat{\lambda }({t_i}) \approx \lambda ({{t_i};{\phi_0}} )+ \lambda ^{\prime}({{t_i};{\phi_0}} )\left\{ {\begin{array}{@{}c@{}} {\dfrac{{{f_s}}}{c}\boldsymbol{n} \cdot \sum\limits_{i = 1}^{{N_p}} {\left[ {\left( {{\boldsymbol{\varPhi }_{rr}}\left( {i\dfrac{{{T_{obs}}}}{{{N_p}}},{t_0}} \right) - {\boldsymbol{I}_{3 \times 3}}} \right)\dfrac{1}{{{N_p}}}} \right]\delta \boldsymbol{r}({t_0})} }\\ {+ \dfrac{{{f_s}}}{c}\boldsymbol{n} \cdot \sum\limits_{i = 1}^{{N_p}} {\left[ {{\boldsymbol{\varPhi }_{rv}}\left( {i\dfrac{{{T_{obs}}}}{{{N_p}}},{t_0}} \right)\dfrac{1}{{{N_p}}}} \right]\delta \boldsymbol{v}({t_0})} } \end{array}} \right\}\end{equation}

The period of the pulsar is much smaller than the whole observation time. It can therefore be considered that ${N_p} \to \infty$, then we can obtain

(13)\begin{gather} \begin{gathered} \dfrac{{{f_s}}}{c}\boldsymbol{n} \cdot \sum\limits_{i = 1}^{{N_p}} {\left[ {\left( {{\boldsymbol{\varPhi }_{rr}}\left( {i\dfrac{{{T_{obs}}}}{{{N_p}}},{t_0}} \right) - {\boldsymbol{I}_{3 \times 3}}} \right)\dfrac{1}{{{N_p}}}} \right]\delta \boldsymbol{r}({t_0})} \\ \approx \dfrac{{{f_s}}}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rr}}({\tau ,{t_0}} )d\tau } - {T_{obs}}{\boldsymbol{I}_{3 \times 3}}} \right)\delta \boldsymbol{r}({t_0}) \end{gathered} \end{gather}

(14)\begin{gather}\frac{{{f_s}}}{c}\boldsymbol{n} \cdot \sum\limits_{i = 1}^{{N_p}} {\left[ {{\boldsymbol{\varPhi }_{rv}}\left( {i\frac{{{T_{obs}}}}{{{N_p}}},{t_0}} \right)\frac{1}{{{N_p}}}} \right]\delta \boldsymbol{v}({t_0}) \approx } \frac{{{f_s}}}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rv}}({\tau ,{t_0}} )d\tau } } \right)\delta \boldsymbol{v}({t_0})\end{gather}

Make $\varDelta \boldsymbol{\lambda }$ represent the propagation effects of initial error.

(15)\begin{equation}\varDelta \boldsymbol{\lambda } = \frac{{{f_s}}}{{c{T_{obs}}}}\left\{ \begin{array}{@{}l} \boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rr}}({\tau ,{t_0}} )d\tau } - {T_{obs}}{\boldsymbol{I}_{3 \times 3}}} \right)\delta \boldsymbol{r}({t_0})\\ + \boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rv}}({\tau ,{t_0}} )d\tau } } \right)\delta \boldsymbol{v}({t_0}) \end{array} \right\}\left[ {\begin{array}{@{}c@{}} {\lambda^{\prime}({{t_1};{\phi_0}} )}\\ {\lambda^{\prime}({{t_2};{\phi_0}} )}\\ \vdots \\ {\lambda^{\prime}({{t_{{N_b}}};{\phi_0}} )} \end{array}} \right]\end{equation}

The additional phase estimation error caused by $\varDelta \boldsymbol{\lambda }$ can be written as

(16)\begin{equation}\varDelta \hat{\phi } = {({{\boldsymbol{H}^T}\boldsymbol{H}} )^{ - 1}}{\boldsymbol{H}^T}\varDelta \boldsymbol{\lambda }\end{equation}

where

(17)\begin{equation}\boldsymbol{H} = \boldsymbol{\lambda ^{\prime}}({\phi _0})\end{equation}

Substitute $\varDelta \boldsymbol{\lambda }$ into Equation (16), then we can derive the expression of additional phase estimation error caused by the propagation effects of initial error.

(18)\begin{align}\begin{aligned} \varDelta \hat{\phi } & =\dfrac{{{f_s}}}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rr}}({\tau ,{t_0}} )d\tau } - {T_{obs}}{\boldsymbol{I}_{3 \times 3}}} \right)\delta \boldsymbol{r}({t_0})\\ & \quad + \dfrac{{{f_s}}}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rv}}({\tau ,{t_0}} )d\tau } } \right)\delta \boldsymbol{v}({t_0}) \end{aligned}\end{align}

where $\varDelta \hat{\phi }$ represents the distort phase difference between the modified observation profile and the directly recovered one.

From simple simulation, the value of the term $\frac{{{f_s}}}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rr}}({\tau ,{t_0}} )d\tau } - {T_{obs}}{\boldsymbol{I}_{3 \times 3}}} \right)\delta \boldsymbol{r}({t_0})$ is extremely small and can be neglected in Equation (18). For the vehicle in Earth orbit, for every 10 km increase in the position error, the value only increases by about 1 × 10⁻⁵, which is much lower than the phase difference that can be resolved by signal processing (Wang, Reference Wang2016).

Finally, the constructed template-independent measurement model can be simplified and expressed as follows, which denotes the estimation error of TOA attributed to the velocity error.

(19)\begin{equation}\varDelta TO{A_2} \approx \frac{1}{{c{T_{obs}}}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^T {{\boldsymbol{\varPhi }_{rv}}({\tau ,{t_0}} )d\tau } } \right)\delta \boldsymbol{v}({t_0})\end{equation}

3. Template-independent navigation method

According to Equation (19), measurements with respect to the velocity error of the vehicle can be obtained without using the pulsar template. However, navigation is not available with access to these measurements of velocity alone, the system is unobservable, the position information of vehicle cannot be calculated.

Extending Equation (11), the state error of the vehicle at an arbitrary epoch t can be predicted by propagating the initial state error by using the orbital transition matrix (Liu, Reference Liu1992).

(20)\begin{equation}\left\{ {\begin{array}{@{}c} {\delta \boldsymbol{r}(t )= {\boldsymbol{\Phi }_{rr}}(t,{t_0})\delta \boldsymbol{r}({{t_0}} )+ {\boldsymbol{\Phi }_{rv}}(t,{t_0})\delta \boldsymbol{v}({{t_0}} )}\\ {\delta \boldsymbol{v}(t )= {\boldsymbol{\Phi }_{vr}}(t,{t_0})\delta \boldsymbol{r}({{t_0}} )+ {\boldsymbol{\Phi }_{vv}}(t,{t_0})\delta \boldsymbol{v}({{t_0}} )} \end{array}} \right.\end{equation}

where ${\boldsymbol{\Phi }_{vr}}(t,{t_0})$ and ${\boldsymbol{\Phi }_{vv}}(t,{t_0})$ are the other two partition matrices of the orbital transition matrix.

As shown in Equation (21), in order to achieve system observability, the initial state errors of vehicle can be completely transformed into velocity errors at arbitrary epochs by performing multi-segment observations.

(21)\begin{equation}\left\{ {\begin{array}{@{}c} {\delta {\boldsymbol{v}_1} = {\boldsymbol{\varPhi }_{vr}}({{T_1},{t_0}} )\delta {\boldsymbol{r}_0} + {\boldsymbol{\varPhi }_{vv}}({{T_1},{t_0}} )\delta {\boldsymbol{v}_0}}\\ {\delta {\boldsymbol{v}_2} = {\boldsymbol{\varPhi }_{vr}}({{T_2},{t_0}} )\delta {\boldsymbol{r}_0} + {\boldsymbol{\varPhi }_{vv}}({{T_2},{t_0}} )\delta {\boldsymbol{v}_0}}\\ {\ldots \ldots }\\ {\delta {\boldsymbol{v}_j} = {\boldsymbol{\varPhi }_{vr}}({{T_j},{t_0}} )\delta {\boldsymbol{r}_0} + {\boldsymbol{\varPhi }_{vv}}({{T_j},{t_0}} )\delta {\boldsymbol{v}_0}} \end{array}} \right.\end{equation}

where ${T_j},j = 1,2, \ldots ,n$ represent the end moments of multi-segment observations.

Then we can establish the navigation equations with the term of

(22)\begin{align}\left\{ \begin{aligned} \varDelta {\phi_1} & = \dfrac{{{f_s}}}{{c({{T_1} - {t_0}} )}}\boldsymbol{n} \cdot \left( {\int_{{t_0}}^{{T_1}} {{\boldsymbol{\varPhi }_{rv}}({\tau ,{t_0}} )d\tau } } \right)\delta {\boldsymbol{v}_0}\begin{array}{*{20}{c}} {}& {} \end{array}\\ \varDelta {\phi_2} & = \frac{{{f_s}}}{{c({{T_2} - {T_1}} )}}\boldsymbol{n} \cdot \left( {\int_{{T_1}}^{{T_2}} {{\boldsymbol{\varPhi }_{rv}}({\tau ,{T_1}} )d\tau } } \right)\delta {\boldsymbol{v}_1}\begin{array}{*{20}{c}} {}& {} \end{array}\\ & \hspace{6pc}{\ldots\ldots}\\ \varDelta {\phi_n} & = \frac{{{f_s}}}{{c({{T_n} - {T_{n - 1}}} )}}\boldsymbol{n} \cdot \left( {\int_{{T_{n - 1}}}^{{T_n}} {{\boldsymbol{\varPhi }_{rv}}({\tau ,{T_{n - 1}}} )d\tau } } \right)\delta {\boldsymbol{v}_{n - 1}} \end{aligned}\right.\end{align}

To meet the number of state elements, at least six observation segments are required.

Combining Equations (21) and (22), the simultaneous equations can be formed as

(23)\begin{equation}\varDelta \varPhi = A \cdot \left[\begin{array}{@{}c@{}} {\delta {\boldsymbol{r}_0}}\\ {\delta {\boldsymbol{v}_0}} \end{array} \right]\end{equation}

where $\varDelta \varPhi$ expresses the matrix of distort phase differences obtained by dynamic signal processing, A denotes the calculated orbital transition coefficient matrix of the navigation system.

To solve the simultaneous equations, we can build a function optimisation problem, which can be solved easily by the common optimisation algorithms (i.e., genetic algorithm, differential evolution or particle swarm optimisation) (Kennedy and Eberhart, Reference Kennedy and Eberhart1995; Hang and Hang, Reference Hang and Hang2003; Storn and Price, Reference Storn and Price2006). The strategies of crossover and mutation are improved adaptively to keep diversity of samples and improve the computation efficiency.

In summary, the proposed pulsar-based template-independent navigation method can be complemented by the following procedure (Figure 2):

1. Calculate the distort phase difference between the modified observation profile and the directly recovered profile for each data segment by dynamic signal processing.

2. Set the elements of $\delta {\boldsymbol{v}_0},\delta {\boldsymbol{r}_0}$ as variables, substitute them into Equations (21) and (22) to calculate the orbital transition coefficient matrix A with the aid of predicted state information of the vehicle, and then establish the navigation equations.

3. Solve the constructed equations, the information of orbit determination or navigation can be achieved directly.

Figure 2. Flow chart of the proposed method

4. Simulation results

In this section, the simulation data and real data from HXMT and NICER of the Crab pulsar are employed to verify the performance of the proposed method. For simulation data, when the propagation of orbital dynamic is performed, the two-body model is utilised, and the data are generated by adding different initial state errors. For real data, the predicted orbit information used is propagated by adding different initial errors based on the accurate orbit model (Wang, Reference Wang2016).

4.1 Analysis of theoretical model

4.1.1 Problem of Crab template

In this part, the real data from NICER are utilised to simulate the situation of building the high precision templates in orbit. The phase propagation errors between the templates with different time intervals are analysed.

In Figure 3, the two templates are built separately by the series of photon TOAs detected during September and November 2018. The estimation value of the phase propagation with the time interval is calculated according to the time and phase model, and the first template is translated to the folding moment of the second template. We can see that there exists a certain phase difference between the two templates.

Figure 3. Phase propagation error between the templates with intervals of two months

Figure 4 shows the variation of the phase propagation error between the constructed templates with different time intervals. The phase propagation error which may result in an extra position error increases sharply with the time interval. Therefore, it can be concluded that, for the Crab pulsar, the pre-calibrated template cannot be used in the long term and the timing parameters must be updated quite frequently (generally within 10 days).

Figure 4. Phase propagation errors between the templates with different intervals

4.1.2 Feasibility of the derived measurement model

In Section 2, a template-independent measurement model is derived as a term of the velocity error of the vehicle. To verify that the derived model is correct, it is investigated via simulation.

Figure 5 shows the comparison between the signal processing results and the calculation results obtained by the derived theoretical model. The minimum phase resolution of signal processing is 1 × 10⁻⁴. We can see that the set of theoretical results is consistent with the signal processing results. For the comparison result, the mean error is about 1⋅504 × 10⁻⁴, and the standard variance is about 1⋅833 × 10⁻⁴.

Figure 5. Comparison between the theoretical and the signal processing results with different velocity of vehicle

Therefore, it can be concluded that the measurement model derived in this paper is correct and can be applied to orbit determination or navigation.

4.2 Analysis of the proposed method

4.2.1 Experiments with simulation data

The parameters of the simulated Crab pulsar are listed in Table 1. The pulsar is assumed to be stationary during the observation period and the effect of higher order derivatives of frequency is neglected.

Table 1. Simulation parameters of Crab pulsar

The relationship between the accuracy of distort pulse TOA and the observation time is analysed first, and the result is shown in Figure 6. We can see that the accuracy of the distort pulse TOA increases with observation time, and when the observation time reaches 600–1,000 s, this variation will slow down.

Figure 6. Accuracy of distort pulse TOA in relation to the observation time

In this situation, therefore, 20 groups of simulation data were generated with the whole observation period of 6,000 s (each segment about 1,000 s), and the added different initial state errors were randomly produced by uniform distribution.

Figure 7 shows that the navigation result of simulation data, the average estimation error of position and velocity are about 7⋅23 km and 7⋅62 m/s. This result demonstrates that the proposed template-independent method is feasible, and can provide the service of orbit determination or navigation with relatively high accuracy.

Figure 7. Estimation error of the proposed method for the simulation data

The superiority of the proposed method over the existing pulsar-based navigation methods is further investigated. Table 2 shows the comparison of navigation performance between the traditional filtering algorithm and the proposed method under the same simulation conditions. It can be seen that the performance of the methods is basically the same. The computational cost of the proposed method is slightly higher but still very low. It should be noted that we can obtain the navigation information directly by this proposed method. For the traditional filtering algorithm, however, it usually takes time to converge the current navigation accuracy.

Table 2. Comparison between traditional filtering algorithm and the proposed method

For the existing pulsar-based navigation methods, pulsar template is the basic input. The phase propagation error existing in the template will directly cause a corresponding extra error in the processing result of the pulse TOA. Considering that the true state of the vehicle cannot be accurately obtained, the template constructed in orbit will also contain a certain phase error. Based on Figure 4, if the template and the timing parameters of the Crab pulsar were not updated in time, the resulting extra position error, which is in proportion, will exceed 10 km. The proposed template-independent method in this paper is not affected by this.

4.2.2 Experiments with real data

In this part, real data from the Crab pulsar are employed to verify the performance of the proposed method. The initial frequency parameters can be found in the database of Jodrell Bank, UK.

The real observation data from HXMT on 31 December 2018 (Obs_ID: 0101299008) was utilised. Ten groups of data detected by the onboard high energy sensor were chosen with an observation period of more than 3,000 s and Table 3 shows the specifics. Figure 8 shows the navigation results, and the average estimation error of position and velocity are about 9⋅53 km and 9⋅89 m/s.

Figure 8. Estimation error of the proposed method for the real data from HXMT

Table 3. Specifics of real observation data from HXMT

Real observation data from NICER on 26 December 2018 (Obs_ID: 1013010147) were then utilised. Twelve groups of data with observation period more than 2,000 s were chosen and Table 4 shows the specifics.

Table 4. Specifics of real observation data from NICER

For the proposed method, it is necessary to divide the chosen data into six segments, which will result in less data per segment. This will have an impact on the navigation performance. Figure 9 shows the navigation results; the average estimation error of position and velocity are about 8⋅94 km and 9⋅37 m/s, respectively.

Figure 9. Estimation error of the proposed method for the real data from NICER

For the real observation data, most groups of the estimation results of position and velocity are better than 10 km and 10 m/s. This also verifies the feasibility of the proposed method. Compared with the navigation results by processing simulation data, the navigation accuracy is lower. This is because, for the model derived in Section 2, the variation of ${f_s}$ is not considered, which will have an impact on calculating the distort pulse TOA and come to affect the navigation performance. Figure 10 shows the variation of the estimation error of $\varDelta \varPhi$ by processing simulation and real data. We can see that the phase estimation error increases significantly with time for the real data, and basically remains unchanged for the simulation data.

Figure 10. Variation of phase estimation error for using different kinds of data