3.1. X-ray pulsar signal model for orbiting vehicles

The process for a vehicle to record a series of X-ray photon TOAs can be modelled as a nonhomogeneous Poisson process with a periodic rate function $\lambda \lpar t\rpar \ge 0$ . λ (t) indicates the aggregate rate of photons from pulsar and background, and is of the form (Emadzadeh and Speyer, Reference Emadzadeh and Speyer2010a)

(1) $$\lambda \lpar t \rpar =\alpha h\lpar {\phi_{\rm det} \lpar t \rpar } \rpar +\beta$$

where $h\lpar \cdot\rpar $ is the periodic pulsar profile, $\phi_{\rm det} \lpar t\rpar $ is the detected phase, and α and β are the detected rate constant of the pulsar and of the background respectively.

Taking the orbital motion of the vehicle into consideration, $\phi_{\rm det} \lpar t \rpar $ can be modelled as

(2) $$\phi_{\rm det} \lpar t \rpar =\phi_0 +f_s \lpar {t-t_0} \rpar +\displaystyle{{f_s}\over{c}}\int_{t_0}^t {v\lpar \tau \rpar } d\tau$$

where ϕ₀ is the detected phase at initial time t ₀, f _s is the intrinsic frequency of the pulsar, c is the speed of light, and v is the velocity of the vehicle towards the pulsar.

When the vehicle has a velocity of ${\bi v}$ and the direction vector of the pulsar is ${\bi n}$ , Equation (2) can be rewritten as

(3) $$\phi_{\rm det} \lpar t \rpar =\phi_0 +f_s \lpar {t-t_0} \rpar +\displaystyle{{f_s}\over{c}}\int_{t_0}^t {{\bi n}\cdot {\bi v}\lpar \tau \rpar } d\tau$$

3.2. Epoch folding

Epoch folding that could recover the pulse profile from a series of recorded photon TOAs has been widely applied in the astronomy field.

Suppose an observation period lasting for T _obs seconds consists of N _p periods of the received signal, P. P can be further divided into N _b bins, the size of which is T _b seconds. Epoch folding proceeds as follows. Firstly, all the photons are folded back into the first period according to their TOAs. Then, the photons collected in each bin are counted. Finally, an empirical profile can be recovered by normalising the computed photon counts.

The empirical profile in the ith bin i ∈ [1, N _b ], $\tilde {\lambda}\lpar T_{i}\rpar $ , can be described by (Emadzadeh and Speyer Reference Emadzadeh and Speyer2010b)

(4) $$\tilde {\lambda}(T_i)=\lambda (T_i,\phi_0)+n(T_i)$$

where $\lambda \lpar T_{i}\comma \; \phi_{0}\rpar $ is the template with an initial phase ϕ₀ and n(T _i ) is the epoch folding noise.

When the empirical profile is built, ϕ₀ can be estimated by comparing the template and the empirical profile by some computationally efficient methods such as the cross-correlation method or the fast near-maximum likelihood estimator (Rinauro et al., Reference Rinauro, Colonnese and Scarano2013). Figure 3 describes the whole process of the epoch folding. It is obvious that the epoch folding method relies on the assumption that the period of the signal is constant. Thus, the classical epoch folding method cannot be directly employed by orbiting vehicles, as the period of the received signal of those vehicles would vary with time.

Figure 3. Diagram of the epoch folding.

3.3. Linearized method

In Wang and Zheng (Reference Wang and Zheng2016), we rewrote the integral term on the right side of Equation (3) as

(5) $$\int_{t_0}^t {{\bi n}\cdot {\bi v}\lpar \tau \rpar } d\tau ={\bi n}\cdot {\bi r}\lpar t \rpar -{\bi n}\cdot {\bi r}\lpar {t_0} \rpar $$

where ${\bi r}\lpar \cdot \rpar $ is the position of the vehicle.

As ${\bi r}\lpar \cdot \rpar $ can be expressed as the sum of the position predicted by the orbital dynamics of vehicle $\tilde {{\bi r}}\lpar \cdot \rpar $ and the error within $\tilde {{\bi r}}\lpar \cdot \rpar $ , $\delta {\bi r}\lpar \cdot \rpar $ , and the transition matrix that connects the $\delta \,{\bi r}\lpar \cdot \rpar $ s can be written in a polynomial form, we have the linearized phase propagation model of

(6) $$\phi_{\rm det} \lpar t \rpar =\phi_0 +f_s \lpar {t-t_0} \rpar +\sum\limits_{k=1}^\infty {f_{dy}^{(k)} \lpar {t-t_0} \rpar ^k} +\displaystyle{{f_s}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}\lpar t \rpar -\tilde {{\bi r}}\lpar {t_0} \rpar } \rsqb $$

where $f_{dy}^{\lpar k\rpar }$ is the coefficients of the polynomial that need to be estimated.

For Earth-orbiting vehicles, Equation (6) can be rewritten as

(7) $$\eqalign{ \phi _{{\rm det}}(t) =& \phi _0 + f_s(t - t_0) + \sum\limits_{k = 1}^\infty {f_{dy}^{(k)} {(t - t_0)}^k} + \displaystyle{{f_s} \over c}{\bf n}\cdot [\widetilde{{\bf r}}_{SC/E}(t) - \widetilde{{\bf r}}_{SC/E}(t_{\rm 0})] \cr & + \displaystyle{{f_s} \over c}{\bf n}\cdot [\widetilde{{\bf r}}_{E}(t{\rm )} - \widetilde{{\bf r}}_{E}(t_{0}{\rm )}]} $$

where $\tilde {{\bi r}}_{SC/E} \lpar \cdot \rpar $ is the position of the vehicle with respect to the Earth and $\tilde {{\bi r}}_{E} \lpar \cdot \rpar $ is the position of the Earth with respect to the SSB that could be provided by the Jet Propulsion Laboratory (JPL) ephemeris such as DE405 or DE421.

4.1. Velocity acquisition via starlight Doppler measurements

Assuming a spectrometer measures the Doppler frequency shift of the ith star (i = 1, 2, 3,.....) and obtains the velocity of the vehicle towards the ith star, v _i , the measurement model for velocity acquisition is (Liu et al., Reference Liu, Fang, Ma, Kang and Wu2015)

(8) $$v_i ={\bi s}_i ^T{\bi v}_{SC/E} +{\bi s}_i ^T{\bi v}_E -v_{s,i} +\zeta_i$$

where $s_{i} \in {\open R}^{{3} \times {1}}$ is the direction vector of the ith star, v _{s, i} is the radial velocity of the ith star, ${\bi v}_{SC/E} $ is the velocity of the vehicle relative to the Earth, ${\bi v}_{E} $ is the velocity of the Earth that could also be provided by DE405, and ζ _i is the measurement noise.

Assuming three stars are observed, Equation (8) can be rewritten in a matrix form of

(9) $${\bi Z}={\bi S}\lpar {{\bi v}_{SC/E} +{\bi v}_E} \rpar +{\bi Z}_s +{\bf{\zeta}}$$

where ${\bi Z}=[ {v_{1}} \quad {v_{2}} \quad {v_{3}}] ^{T}$ , ${\bi S}=[ {{\bi s}_{1}^{T}} \quad {{\bi s}_{2}^{T}} \quad {{\bi s}_{3}^{T}}]^{T}$ , ${\bi Z}_{s} =[ {v_{s\comma 1}} \quad {v_{s\comma 2}} \quad {v_{s\comma 3}}]^{T}$ , and ${{\bf{\zeta}}=[{\zeta_{1}} \quad {\zeta_{2}} \quad {\zeta_{3}}\rsqb ^{T}}$ .

Thus, ${\bi v}_{SC/E} $ can be estimated by

(10) $${\bf{\hat {v}}}_{SC/E} =\lpar {{\bi S}^T{\bi R}^{-1}{\bi S}} \rpar ^{-1}{\bi S}^T{\bi R}^{-1}\lpar {{\bi Z}-{\bi Z}_s} \rpar -{\bi v}_E$$

where ${\bi R}$ is the covariance of ${\bf{\zeta}}$ . The covariance of ${\bf{\hat {v}}}_{SC/E} $ is ${cov({\hat v}_{SC/E}) = ({\bf S}^{T}{\bf R}^{ - {1}}{\bf S})^{-1}}$ .

4.2. Starlight Doppler measurement-aided phase propagation model

4.2.1. Derivation

In contrast to the way the linearization method modifies the integral term in Equation (3), in this subsection, we solve the integral term by numerically integrating rather than analytically solving it. For Earth-orbiting vehicles, we have

(11) $$\int_{t_0}^t {{\bi n}^T{\bi v}\lpar \tau \rpar } d\tau =\int_{t_0}^t {{\bi n}^T{\bi v}_{SC/E} \lpar \tau \rpar } d\tau +{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb $$

As the second term on the right hand of Equation (11) is known, we focus on solving the first term. It is well known that $\int_{t_{0}}^{t} {{\bi n}^{T}{\bi v}_{SC/E} \lpar \tau \rpar } d\tau $ can be solved using the composite trapezoid formula (Hildebrand, Reference Hildebrand1974). Assuming the observation period $[ {t_{0}} \quad t]$ consists of N sampling periods of the starlight Doppler measurement, T _s , we have

(12) $$\int_{t_0}^t {{\bi n}^T{\bi v}_{SC/E} \lpar \tau \rpar } d\tau \approx \displaystyle{{T_s}\over{2}}\sum\limits_{j=1}^{N-1} {\lpar {{\bi n}^T{\bi v}_{SC/E,j} +{\bi n}^T{\bi v}_{SC/E,j+1}} \rpar }$$

where ${\bi v}_{SC/E\comma j} $ is the true velocity of vehicle at the end epoch of the jth sampling period of the starlight Doppler measurement.

Given that the true velocity of the vehicle is unknown, we substitute the estimated velocity ${\bf{\hat {v}}}_{SC/E} $ into Equation (12), yielding

(13) $$\int_{t_0}^t {{\bi n}^T{\bi v}_{SC/E} \lpar \tau \rpar } d\tau \approx \displaystyle{{T_s}\over{2}}\sum\limits_{j=1}^{N-1} {\lpar {{\bi n}^T{\bf{\hat {v}}}_{SC/E,j} +{\bi n}^T{\bf{\hat {v}}}_{SC/E,j+1} +{\bi n}^T{\bf{\xi}}_j +{\bi n}^T{\bf{\xi}}_{j+1}} \rpar }$$

where ${\bf{\xi}}_{j} $ is the error within ${\bf{\hat {v}}}_{SC/E\comma j} $ and can be modelled as a zero-mean Gaussian noise.

When N approaches infinity, we have

(14) $$\displaystyle{{T_s}\over{2}}\sum\limits_{j=1}^{N-1} {\lpar {{\bi n}^T{\bf{\xi}}_j +{\bi n}^T{\bf{\xi}}_{j+1}} \rpar } \approx 0$$

In practice, the sampling period of starlight Doppler measurement is in the order of a second and the observation period of a pulsar is in the order of hours. Thus, the assumption that N is approximated to be infinite readily holds in real space missions.

Substituting Equations (13) and (14) into Equation (3) yields

(15) $$\phi_{\rm det} \lpar t \rpar =\phi_0 +f_s \lpar {t-t_0} \rpar +\displaystyle{{f_s T_s}\over{2c}}\sum\limits_{j=1}^{N-1} {\lpar {{\bi n}^T{\bf{\hat {v}}}_{SC/E,j} +{\bi n}^T{\bf{\hat {v}}}_{SC/E,j+1}} \rpar } +\displaystyle{{f_s}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb $$

Equation (15) could be used as the starlight Doppler measurement-aided phase propagation model, if the sampling period of starlight Doppler measurement is close to the difference between the recorded TOAs of photons next to each other. However, in practice, the sampling period of the starlight Doppler measurement is in the order of a second, but the difference between the TOAs of photons next to each other is usually shorter than the intrinsic period of a pulsar, which is usually in the order of a millisecond. This means there would be many photon TOAs that are not employed to estimate the pulse phase if Equation (15) were used as it would reduce the phase estimation accuracy. In order to introduce those “lost” photon TOAs into the pulse phase estimation, we use linear interpolation to modify Equation (15).

When $t_{p} \in [ {jT_{s}} \quad {\lpar {j+1} \rpar T_{s}}\rsqb $ the velocity of the vehicle at t _p can be approximated as a term of ${\bf{\hat {v}}}_{SC/E\comma j} $ and ${\bf{\hat {v}}}_{SC/E\comma j+1}$ ,

(16) $${\bi v}_p =\left({1-\displaystyle{{t}\over{T_s}}} \right){\bf{\hat {v}}}_{SC/E,j} +\displaystyle{{t}\over{T_s}}{\bf{\hat {v}}}_{SC/E,j+1}$$

Then, we have

(17) $$\eqalign{\phi_{\rm det} \lpar {t_p} \rpar =&\phi_0 +f_s \lpar {t_p -t_0} \rpar +\displaystyle{{f_s T_s}\over{2c}}\sum\limits_{k=1}^{j-1} {\lpar {{\bi n}^T{\bf{\hat {v}}}_{SC/E,k} +{\bi n}^T{\bf{\hat {v}}}_{SC/E,k+1}} \rpar} \cr & +\displaystyle{{f_s}\over{2c}}\lpar {t-jT_s} \rpar \lpar {{\bi n}^T{\bf{\hat {v}}}_{SC/E,j} +{\bi n}^T{\bi v}_p} \rpar +\displaystyle{{f_s}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb}$$

In summary, the starlight Doppler measurement-aided phase propagation model is

(18) $$\phi_{\rm det} \lpar t \rpar =\phi_0 +f_s \lpar {t-t_0} \rpar +\displaystyle{{f_s}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb +\phi_{star} \lpar t \rpar $$

where

(19) $$\phi _{star}{\rm }(t{\rm )} = \left\{ {\matrix{{\displaystyle{{f_sT_s} \over {2c}}\sum\limits_{k = 1}^{j - 1} {({\bi n}^{T}{{\hat {\bi v}}}_{{SC}/{E},{k}} + {\bi n}^{T}{{\hat {\bi v}}}_{{SC}/{E},{k} + {1}})} } & {t = jT_s} \cr {\displaystyle{{f_sT_s} \over {2c}}\sum\limits_{k = 1}^{j - 1} {({\bi n}^{T}{{\hat {\bi v}}}_{{SC}/{E},{k}} + {\bi n}^{T}{{\hat {\bi v}}}_{{SC}/{E},{k} + {1}})} } & {t \in \left[ {\matrix{ {jT_s} & {(j + 1)T_s} \cr } } \right]} \cr {\quad + \displaystyle{{f_s} \over {2c}}(t - jT_s){\rm }({\bi n}^{T}{{\hat {\bi v}}}_{{SC}/{E},{j}} + {\bi n}^{T}{\bi v}_{p})}}} \right.$$

and $j=0\comma \; 1\comma \; 2\comma \; \cdots \cdots\comma \; N-1$ .

4.2.2. Removal of the orbital effect within the recorded photon TOAs

As illustrated in Section 3.2, the epoch folding operates on the premise that the investigated signal is of constant frequency. This subsection will show how to remove the time-varying orbital effect within the recorded photon TOAs and to have photon TOAs with a constant frequency.

Equation (18) can be rewritten as

(20) $$\displaystyle{{1}\over{f_s}}\phi_{\rm det} \lpar t \rpar -\displaystyle{{\phi_0}\over{f_s}}=t-t_0 +\displaystyle{{1}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb +\displaystyle{{1}\over{f_s}}\phi_{star} \lpar t \rpar $$

Setting $\bar {t}=\phi_{\rm det} \lpar t \rpar /f_{s} $ and $\bar {t}_{0} =\phi_{0} /f_{s}$ , we have

(21) $$\bar {t}-\bar {t}_0 =t-t_0 +\displaystyle{{1}\over{c}}{\bi n}\cdot [ {\tilde {{\bi r}}_E \lpar t \rpar -\tilde {{\bi r}}_E \lpar {t_0} \rpar } \rsqb +\displaystyle{{1}\over{f_s}}\phi_{star} \lpar t \rpar $$

Assuming there are M photon TOAs recorded within an observation period of a pulsar, the photon TOAs series is denoted as $\lcub {t_{i}} \rcub _{i\, =\, 1}^{M} $ . When Equation (21) is applied to all the photon TOAs, the orbital effect within $\lcub {t_{i}} \rcub _{i\, =\, 1}^{M} $ can be roughly removed, bringing in new photon TOAs, ${\lcub {\bar {t}_{i}} \rcub}_{i=1}^{M} $ . Compared with ${\lcub {t_{i}} \rcub}_{i=1}^{M}$ , the period of ${\lcub {\bar {t}_{i}} \rcub}_{i=1}^{M} $ is the intrinsic period of the pulsar, i. e., 1/f _s . Thus, the initial phase of ${\lcub {\bar {t}_{i}} \rcub}_{i=1}^{M} $ can be readily estimated by the classical epoch folding methods.

4.3. Summary of the integrated pulse phase estimation method

The proposed integrated pulse phase estimation method can be complemented by the following procedure:

• Step 1 Estimate the velocity of the vehicle at each sampling period of the starlight Doppler measurements according to Equation (10);

• Step 2 Remove the orbital effect in the photon TOAs $\lcub {t_{i}} \rcub _{i=1}^{M} $ using Equation (21), and obtain ${\lcub {\bar {t}_{i}} \rcub}_{i=1}^{M} $ ;

• Step 3 Calculate the initial phase of ${\lcub {\bar {t}_{i}} \rcub}_{i=1}^{M} $ using epoch folding methods.

5.1. Simulation conditions

Three Earth orbiting satellites are investigated, and the initial elements of the three satellites are shown in Table 1. The pulsar PSR B1821-24 with parameters listed in Table 2 is employed as the observed pulsar. The template of PSR B1821-24 is given in Figure 5 (Rutledge et al., Reference Rutledge, Fox, Kullkarni, Jacoby, Cognard, Backer and Murray2004). The ephemeris DE405 is employed to predict the position and velocity of Earth within the observation period of the pulsar. The fast-near maximum likelihood estimator is used to estimate the initial phase.

Figure 5. PSR B1821-24 profile.

Table 1. Initial orbital elements of satellites.

Table 2. Parameters of Simulated PSR B1821-24.

The sampling period of starlight Doppler measurement is 0·1 s and the measurement accuracy is 0·1 m/s. The navigation stars, the data of which is taken from Liu et al. (Reference Liu, Fang, Ma, Kang and Wu2015), are listed in Table 3.

Table 3. Navigation stars.

Root Mean Square (RMS) error is selected to evaluate the initial phase estimation result. The phase estimation results are obtained from 1000 Monte Carlo trials.

5.2. Hardware-in-loop system validation

The simulations are performed by the newly-developed hardware-in-loop system. Figure 6 provides a block diagram of the hardware-in-loop system. As a vehicle-borne experiment is difficult to perform, it is an effective and economical way to develop a hardware-in-loop system to simulate the actual space mission and to analyse the performance of the proposed integrated navigation system.

Figure 6. Block diagram of the hardware-in-loop system.

The system is mainly constituted of the following parts:

• (1)  X-ray signal simulation portion: Generate the X-ray photons according to the nominal orbital information of the investigated vehicle, and detect the photon TOAs by employing a Silicon Drift Detector (SDD).   It is worth noting that the generation and detection of the X-ray photons are performed in a vacuum pipeline with a length of 3 m, as the X-ray radiation would rapidly attenuate in the atmosphere.

• (2)  Starlight Doppler measurement simulation computer: Simulate the starlight Doppler measurement according to the nominal orbital information of the investigated vehicle as well as the ephemeris of the observed stars.

• (3)  Fusion unit: perform the proposed integrated pulse phase estimation method to estimate the pulse TOA.

• (4)  High-fidelity orbital dynamics simulation computer: Simulate a high-fidelity space dynamics environment to provide the nominal orbital information to the starlight Doppler measurement simulation computer and the X-ray signal simulation portion.

• (5)  Integrated navigation simulation computer: Calculate the position and velocity of the investigated vehicle using the pulse TOA provided by the fusion unit.

As the key portion of the integrated navigation system is the fusion unit, we will first analyse the performance of the integrated pulse phase estimation method.

Figure 7 shows the phase estimation result of the proposed method for the three satellites. For the three satellites, the RMS errors of estimated phase all reduce as the observation period of a pulsar increases. When the observation period is less than a threshold of around 100 s, the RMS errors experience an unreliable region, and the phase estimation result should not be used. When the observation period exceeds the threshold, the RMS errors become linearly convergent to the Cramér-Rao Bound (CRB) and finally reach an accuracy of about 0·001. This means the proposed method could effectively reduce the orbital effect. In addition, the estimation accuracy grows as the orbital altitude of the satellite increases. Thus, we will first select satellite 1 to analyse the factors that might affect the performance of the proposed method.

Figure 7. Initial phase estimation results for three satellites.

Figure 8 shows the pulse phase estimation results with different sampling periods of starlight Doppler measurements. As the observation period shorter than 100 s is the unreliable region, Figure 8 merely provides the results when the observation period is greater than 100 s. When the sampling period is greater than 1s, the RMS error curves all diverge from the CRB. This is because the sampling period growth is accompanied with an increasing approximation error of the satellite orbit. As Equation (17) showed, the velocity of the satellite within the sampling period is calculated by linearly interpolating. In other words, the orbit of the satellite is assumed as a line within a sampling period of starlight Doppler measurement. The sampling period continuing to grow would violate the assumption, bringing in an error impossible to ignore. Thus, in order to have a reliable pulse phase estimation result, the sampling period of starlight Doppler measurement should not be longer than 1 s.

Figure 8. RMS error of initial phase with different sampling periods of the starlight Doppler measurement.

Assuming the sampling period of starlight Doppler measurement is 1 s, Figure 9 displays the pulse phase estimation results with different starlight Doppler measurement accuracies ranging from 0·1 m/s to 1 m/s. All the RMS error curves converge to the CRB as the observation period increases, and the difference among the RMS error curves is slight. This means the accuracy of the starlight Doppler measurement puts little impact on the pulse phase estimation. Thus, the proposed method drops the requirement of the accuracy of starlight Doppler measurement and a device with an accuracy of 1 m/s is sufficient.

Figure 9. RMS error of initial phase varying with the starlight Doppler measurement accuracy.

Next, we compare the performance of the proposed integrated phase estimation method to the linearized method. As shown, for the satellites with orbit altitude of higher than 28,000 km (e.g. satellite 3 in this section), the order of the polynomial component of the linearized phase propagation model should at least be 2. In order to make a simple performance comparison, we use satellite 3 to perform the comparison.

Figure 10 compares the phase estimated by the proposed method and that estimated by the linearized method. Both methods converge to the CRB as the observation period increases. The phase estimation results are close to each other. In the simulation environment that contains Intel I7-4710MQ@2·5 GHz and Visual studio 6·0, Figure 11 compares the CPU time cost by the two methods. The linearized method is performed using the 10-threads parallel computation technique and the proposed method is performed using only one thread. Even if the linearized method employs the parallel computation technique, the integrated phase estimation method expends much less CPU time than the linearized method. Thus, the integrated pulse phase estimation method could achieve a phase estimation result similar to the linearized method but requires much less computational complexity.

Figure 10. Initial phase estimation result of the integrated phase estimation method and of the linearized method.

Figure 11. The CPU time cost by the integrated phase estimation method and by the linearized method.

Finally, we investigate the navigation performance of the proposed integrated navigation system. Besides PSR B1821 − 24, PSR B0540 − 69, and PSR B0531 + 21 are adopted. Satellites 1, 2, and 3 are all investigated, and the initial position and velocity error of the satellites are set to be 1 km and 1 m/s along each axis in the Earth-centred Inertial Coordinate System J2000·0. The investigated satellites are assumed to sequentially observe the three pulsars for 1,800 s to perform navigation.

Figure 12 shows the positioning results of the three satellites. The solid lines represent the positioning errors and the dashed lines denote the 3σ outlines. The positioning errors of the three satellites are all within the 3σ outlines, indicating that the proposed integrated navigation method could guarantee a reliable positioning result. As the integrated pulse phase estimation method could save computation resource, the corresponding integrated navigation method is promising.

Figure 12. Positioning performance of the integrated navigation method.