Stellar Angle-Aided Pulse Phase Estimation and Its Navigation Application

: X-ray pulsar-based navigation (XNAV) is a promising autonomous navigation method, and the pulse phase is the basic measurement of XNAV. However, the current methods for estimating the pulse phase for orbiting spacecraft have a high computational cost. This paper proposes a stellar angle measurement-aided pulse phase estimation method for high Earth orbit (HEO) spacecraft, with the aim of reducing the computational cost of pulse phase estimation in XNAV. In this pulse phase estimation method, the effect caused by the orbital motion of the spacecraft is roughly removed by stellar angle measurement. Furthermore, a deeply integrated navigation method using the X-ray pulsar and the stellar angle is proposed. The performances of the stellar angle measure-ment-aided pulse phase estimation method and the integrated navigation method were verified by simulation. The simulation results show that the proposed pulse phase estimation method can handle the signals of millisecond pulsars and achieve pulse phase estimation with lower computational cost than the current methods. In addition, for HEO spacecraft, the position error of the proposed integrated navigation method is lower than that of the stellar angle navigation method.


Introduction
Due to the significant increase in the number of spacecraft in recent years, the burden of ground-based tracking systems has also increased, reducing the survivability of spacecraft [1]. Thus, autonomous navigation is necessary [2]. For Earth-orbit spacecraft, the Global Navigation Satellite System (GNSS) is employed for autonomous navigation and has excellent performance. However, the signals of GNSS are subject to interference, which significantly affects the survivability of spacecraft [3]. Moreover, because the intensity of the GNSS signals decays with the increase in the client spacecraft's orbital altitude, the navigation performance of GNSS is not well enough for spacecraft at orbits higher than 3000 km [4]. Thus, developing an autonomous navigation technique with anti-interference capability for these spacecraft is necessary. Furthermore, due to the increase in the number of spacecraft launches, the required improvement in navigation systems can be understood from the perspective of the new generation of instrumentation for Space Traffic Management (STM) [5][6][7].
XNAV is a promising autonomous navigation technique that has anti-interference characteristics and can be applied in near-Earth and deep space. In this paper, we focus on XNAV for near-Earth spacecraft, which are within the gravitational field of the Earth. Pulsars can provide radiation with stable periodicities because they spin with highly stable periodicities [8]. Although pulsars can provide radiation ranging from optical to gamma frequencies, Emadzadeh and Speyer recommended X-ray radiation for navigation because it can be utilized by small-sized detectors [9]. Since being introduced by Chester and Butman [10], XNAV has experienced significant growth [11,12]. Moreover, flight experiments have been performed. In 2018, the United States announced their SEXTANT (Station Explorer for X-ray Timing and Navigation Technology) flight experiments with a position error less than 10 km [13]. In addition, China also verified the feasibility of XNAV, in principle, using real data from TG-2 (Tiangong-2) Spacelab and the Insight-HXMT (Insight-Hard X-ray Modulation Telescope) [14,15].
The difference between the pulse Time of Arrival (TOA) predicted at the SSB (Solar System Barycenter) with that estimated at the spacecraft can reflect the position of the spacecraft in the direction of pulsar [8]. However, because the signal of the pulsar is very weak, the spacecraft cannot record a continuous pulsed signal but a series of photon TOAs. Thus, one key technique of XNAV is to estimate the pulse TOA using the recorded photon TOAs. For a spacecraft that is stationary or performs a uniform linear motion towards the pulsar, two types of methods are effective for estimating the pulse TOA, namely, the epoch folding method [16] and the maximum likelihood estimator (MLE) method [17]. The epoch folding method compares the template with the empirical profile recovered from photon TOAs and estimates the pulse TOA [16]. Based on epoch folding, a large number of methods for pulse phase estimation have been proposed in recent years. In the weighted average method, the final accurate pulse phase is estimated from the generated estimations, which are derived from a universal set of photon TOAs [18]. Rao et al. proposed an FrFT-based GCC method, which can estimate the pulse phase with higher accuracy than that of the GCC method [19]. In addition, methods have been proposed to reduce the computational cost of the epoch method [20,21]. The MLE method is also well developed. The photon TOAs recorded by the spacecraft follow the non-homogeneous Poisson process [17]. Thus, the MLE method is used to derive a likelihood function and can estimate the pulse TOA by maximizing the likelihood function.
However, in practice, spacecraft perform orbital motions, causing the pulsar signal to have a time-varying Doppler frequency. The velocity and position of the spacecraft are unknown in an autonomous navigation task; thus, it is difficult to reduce the impact of the introduced Doppler frequency. The epoch folding and MLE methods fail when used in this situation. To solve this problem, a phase tracking algorithm was proposed by Golshan and Sheikh [22]. This method divides the pulsar observation period into several intervals within each of which the spacecraft can be assumed to perform a linear uniform motion. Therefore, the Doppler frequency can be assumed to be a constant and can be tracked by a Digital Phase-Locked Loop (DPLL). The performance of the phase tracking method depends on the duration of the interval. In order to assume the spacecraft performs a uniform linear motion within an interval, the interval's duration should be sufficiently short. However, when the duration is shorter than a certain threshold, the pulse phase method cannot obtain a reliable result [23]. The flux of the pulsar determines the threshold. For a faint millisecond pulsar, PSR B1821-24, the threshold is about 100 s [24]. However, the assumption that an Earth-orbiting spacecraft performs a uniform linear motion within 100 s is usually violated. Consequently, the phase tracking method fails to handle the signal of a faint millisecond pulsar.
To overcome the problem of the phase tracking method, a phase estimation method with the aid of the orbital dynamics of the spacecraft was proposed by Wang and Zheng [25,26]. In the remainder of this paper, this method is referred to as the linearized method. The linearized method derives a pulse phase propagation model. In this model, the pulse phase consists of a polynomial term and a predicted term. The linearized method estimates the pulse phase by estimating the parameters of the polynomial term. Because the observation period is not divided into intervals, the linearized method is appropriate for handling the signal of millisecond pulsars. In addition, the Doppler frequency and pulse phase can be simultaneously estimated by a maximum likelihood estimator [27]. Moreover, Xue et al. proposed to estimate the Doppler frequency and pulse phase using the H-test and the fast maximum likelihood method, respectively, thereby partly reducing the computational cost [28]. However, an accurate pulse phase estimation result still needs a large number of photons, which leads to high computational costs. Thus, the above phase estimation methods lead to a rigorous demand for onboard computers.
In addition to the orbit dynamics of spacecraft, the information from other measurements can be used to estimate the pulse TOA. A deeply integrated navigation method that utilizes starlight Doppler measurement was proposed by Wang et al. [24]. Using the starlight Doppler measurement, the effect of orbital motion is roughly removed. In addition, Liu et al. proposed a deeply integrated navigation method that utilizes direction, velocity, and distance measurements [29]. This navigation method uses direction and velocity measurements to remove the effect of orbital motion. However, the spectrometer, the equipment which measures the starlight Doppler, remains underdeveloped, and the measurement accuracies used in the above papers are higher than the current status.
Compared with the spectrometer, the star tracker and the horizon sensor, which can be used to derive the stellar angle measurement, are well developed [30]. Thus, we proposed a deeply integrated navigation method using the X-ray pulsar and stellar angle for high-Earth orbit (HEO) spacecraft.
The integrated navigation method using the X-ray pulsar and stellar angle has been widely studied for Earth-orbit, Mars-orbit, and Jupiter-orbit spacecraft [31][32][33]. In the classical integrated navigation methods using the X-ray pulsar and stellar angle, the pulse TOA is estimated independently, and the pulse TOA and stellar angle measurement are fused to estimate the state of the spacecraft. In contrast, the proposed navigation method employs the stellar angle measurement to estimate the pulse TOA, and thus the integration is more sophisticated than the classical integrated navigation method. We derived a stellar angle-aided pulse propagation model which roughly removes the effect of the orbital motion of the spacecraft. Thus, the pulse phase can be estimated by the epoch folding method. The epoch folding method is used because it is more computationally efficient than MLE. Then, the position of the spacecraft is estimated by the pulse TOA. Compared with the current method for estimating the pulse TOA for orbiting spacecraft, the proposed method can greatly reduce the computational burden. In addition, the proposed integrated navigation approach can also be used for Lunar/Mars satellite missions. Limited by its technical level, the current X-ray detector is extremely large; thus, the proposed method is only applicable to large spacecraft. In the future, with the improvement in detection efficiency, this method can hopefully be applied to nano satellites. This paper is organized as follows. In Section 2, we describe the overall design of the deeply integrated navigation system. We review the relevant works of pulse phase estimation in Section 3. Section 4 outlines the stellar angle measurement-aided pulse phase estimation method. In Section 5, simulations are performed to verify the performance of the stellar angle measurement-aided pulse phase estimation method and integrated navigation method. Finally, the conclusions are given in Section 6.

Overall Design of the Deeply Integrated Navigation System
The composition of the proposed deeply integrated navigation system using the Xray pulsar and stellar angle is shown in Figure 1. It can be seen that the integrated navigation system comprises five units: (1) X-ray detectors, which detect the X-ray radiation signal of the pulsars; (2) atomic clocks, which provide a timestamp for the recorded photons; (3) star trackers and horizon sensors, which provide the stellar angle measurement; (4) a fusion unit, which provides the pulse TOA by processing the photon TOAs, the stellar angle measurement, and the orbital dynamics information of the spacecraft; (5) a navigation unit, which estimates the position of the spacecraft by processing the orbital dynamics information of the spacecraft and the pulse TOA.  Figure 1. Composition of the proposed deeply integrated navigation system using the X-ray pulsar and stellar angle.
As shown in Figure 2, the pulse TOA is estimated independently in the current integrated navigation system. The position of the spacecraft is estimated in the navigation unit by processing the pulse TOA, the stellar angle, and the orbital dynamics information of the spacecraft.  Figure 2. Composition of the current integrated navigation system using the X-ray pulsar and stellar angle.
Compared with the current integrated navigation system, the navigation unit in the proposed deeply integrated navigation system only needs to process the orbital dynamics information and the pulse TOA, which has been previously well studied. According to Ref. [34], the navigation unit is shown in Figure 3. In the navigation unit, the orbital dynamics information and the pulse TOA are processed using a Unscented Kalman Filter (UKF) and the estimated position of the spacecraft is obtained [34].  In the proposed navigation method, the fusion unit is the key unit. Thus, in the following sections, the design of the fusion unit is introduced.

Principle of Pulse Phase Estimation
The recorded pulsar photons obey a non-homogeneous Poisson process with a periodic rate function. ( ) t λ is the rate function and is of the form [35]: The parameters α and β represent the average signal and total background count rates in units of counts per second, respectively. det ( ) t φ represents the detected pulsar signal phase and ( ) h ⋅ is the normalized profile function of the pulsar.
Considering the orbital motion of the spacecraft, det ( ) t φ can be modelled as [24]: where 0 φ is the initial phase at 0 t , c represents the speed of light, s f represents the frequency of the pulsar, n represents the direction vector of the pulsar, and v represents the velocity of the spacecraft.

Epoch Folding Method
The epoch folding method is used to recover the pulse profile from the recorded photons and is a widely applied method. Assuming that the observation period lasts can be described by [16]: where j c is the number of photons in the ith bin folded by the jth period, and i T is the center of the ith bin. After recovering the empirical profile, we can estimate the pulse phase by comparing the empirical profile and the template. The representative methods for epoch folding are the nonlinear least square method, the cross-correlation method, and the fast near-maximum likelihood estimator [36,37]. The epoch folding method is effective under the premise that the pulse period is constant. However, for an orbiting spacecraft, the recorded pulse period varies with time and the epoch folding cannot be directly used.

Stellar Angle Measurement-Aided Pulse Phase Estimation Method
In this section, we take the Earth-orbit spacecraft as an example and show the estimation of the position of the spacecraft using stellar angle measurements and the derivation of the stellar angle measurement-aided pulse phase estimation method.

Position Estimation by Stellar Angle Measurement
As shown in Figure 4, the stellar angle measurement model has the form [31]:

Stellar Angle Measurement-Aided Phase Propagation Model
For a spacecraft in an Earth-orbit, we have [24]: n r r n r r (5) where E r is the position of the Earth and can be obtained from DE405 or DE421. Because E r is known, our task is to solve the first term on the right-hand side of Equation (5).
is unknown, and the estimated /S E r is substituted into Equation (5), giving: Substituting Equations (5) and (6) into Equation (2), yielding: Equation (7) is the stellar angle measurement-aided phase propagation model based on the assumption that the difference between each recorded photon TOA is close to the sampling period of the stellar angle measurement. However, the difference between each recorded photon TOA is considerably shorter than the sampling period of the stellar angle. Thus, if Equation (7) is used to estimate the pulse phase, numerous photons are not applied, thus significantly reducing the accuracy of the pulse phase estimation. To employ all of the photon TOAs, we use a cubic spline interpolation to obtain the estimated position ) , p r can be estimated as: where j S is the cubic polynomial in ,( 1) Then, the stellar angle measurement-aided phase propagation model is: where: n r r n r r (10) We fold the photon TOAs to the last pulse period to ensure the real-time navigation performance, yielding: (11) where: n r r n r r (12) In Equation (11)

Summary of the Stellar Angle Measurement-Aided Pulse Phase Estimation Method
The procedure of the proposed stellar angle measurement-aided pulse phase estimation method is described below: Step 1 Estimate the position of the spacecraft by the stellar angle measurement and orbital dynamics information of the spacecraft; Step 2 Use Equation (11) to remove the orbit effect of the spacecraft in the recorded photon TOAs; Step 3 Estimate the pulse phase end φ in Equation (11) by the epoch folding method.

Accuracy of the Stellar Angle Measurement-Aided Phase Estimation
As shown in Equation (11), the accuracy of det ( ) t φ is depends on the accuracy of ( ) stellar angle end t φ − . We define , and the estimation error of ( ) stellar angle end t φ − has the form: ,   (13) where ⋅ Δ n d and ˆp ⋅ Δ n d are the displacement of the spacecraft in the pulsar direction.
It is obvious that the estimated accuracy of ⋅ Δ n d and ˆp ⋅ Δ n d determines the accuracy of ( ) stellar end t φ − and det ( ) t φ . Thus, the performance of the stellar angle measurementaided pulse phase estimation depends on the estimated accuracy of the displacement of the spacecraft in the direction of the pulsar.

Simulations
We investigated three Earth-orbit spacecraft with the orbital elements shown in Table  1 [24]. We take the PSR B1821-24 pulsar as the observed pulsar, for which the simulation parameters are listed in Table 2 [39]. The template profile of this pulsar is given in Figure  5. We use DE405 to predict the position of Earth. The pulse phase end φ is estimated by the nonlinear least square method using the epoch folding method.
The data of navigation stars are shown in Table 3 [24]. We define the evaluation criterion of the estimation performance as the Root Mean Square (RMS) error of pulse phase end φ from 1000 Monte Carlo trials.  Table 2. Simulation parameters of PSR B1821-24 (Reprinted from [39]).  First, we analyze the performance of the proposed stellar angle measurement-aided pulse phase estimation method. Because the accuracy of the horizon sensors is significantly lower than that of star trackers [30], the accuracy of the stellar angle is mainly determined by the accuracy of the horizon sensors. We assume that the accuracy of the horizon sensors is 0.02°, the accuracy of the star trackers is 3'' [30], the sampling period of stellar angle is 10 s, and the observation period of the pulsar is 3000 s. Figure 6 shows the displacement error of spacecraft 3 in the direction of the pulsar 3 Δd , which is estimated by orbit propagation and the stellar angle measurement. It can be seen that the stellar angle measurement significantly improves the accuracy of 3 Δd . The pulse phase estimation results of the stellar angle measurement-aided pulse phase estimation method for the three spacecraft are shown in Figure 7. The RMS errors decrease with the increase in the observation period of the pulsar. In addition, the RMS errors first experience an unreliable region, and the RMS errors decrease and gradually converge to the Cramer-Rao Lower Bound (CRLB) when the observation period exceeds 100 s. This indicates that the stellar angle measurement-aided pulse phase estimation method reduces the orbit effect. Moreover, the estimation accuracies of the pulse phase reduce with the increase in the orbital altitude of the spacecraft. Because it has the highest orbital altitude, we select spacecraft 3 to analyze the factors that influence the accuracy of the proposed stellar angle measurement-aided pulse phase estimation method.  Assuming that the accuracy of the horizon sensors is 0.02°, the estimation results of the pulse phase with different sampling periods of the stellar angle ranging from 5 to 100 s are shown in Figure 9. When the observation period is shorter than 100 s, the pulse phase estimation result is unreliable. Consequently, only results when the observation period exceeds 100 s are shown. The RMS errors decrease as the observation period increases. In addition, the RMS errors increase with the increase in the sampling period of the stellar angle. When the sampling period of stellar angle exceeds 50 s, there is little improvement in the pulse phase estimation with the increase in the observation time. Consequently, to achieve a suitable estimation performance, we recommend that the sampling period of the stellar angle should be shorter than 50 s. Next, we compare the pulse phase estimation accuracy and computational cost of the linearized method with the proposed stellar angle measurement-aided pulse phase estimation method. We assume that the initial state errors are [5 km, 5 km, 5 km] and [5 m/s, 5 m/s, 5 m/s], the sampling period of the stellar angle is 10 s, the accuracy of the horizon sensor is 0.02°, and the accuracy of the star tracker is 3''. In the linearized method, the third term on the right-hand side of Equation (2) is approximated as a linear function [23]. The pulse phase estimation results of the linearized method and the stellar angle measurement-aided phase estimation method are shown in Figure 10. The phase estimation accuracy of the proposed stellar angle measurement-aided phase estimation method is slightly lower than that of the linearized method. For the PSR B1821-24 pulsar, the disparity in the phase estimation is about 1 × 10 −3 cycle, which would cause an error of position at about 900 m. Moreover, the pulse phase error is a random error in practice. Consequently, the disparity of phase estimation has little effect on the final navigation performance. The computation environment comprises an Intel Core i5-7500@3.4GHZ (Intel, Santa Clara, CA, USA) and Python 3.8.1. We access the CPU time cost as the computational cost. As shown in Figure 11, as the observation period increases, the CPU time for the linearized method increases from 0.45 to 20.9 s, whereas that for the proposed method increases from 0.35 to 0.4 s. Thus, the computational cost of the stellar angle measurementaided phase estimation method is significantly lower than that of the linearized method. Finally, for spacecraft 3, the navigation performance of the stellar angle navigation method is compared with that of the proposed deeply integrated navigation method. We chose PSR B1937+21, PSR B1821-24, and PSR J0218+4232 as the navigation pulsars. We assume that the initial state errors are [5 km, 5 km, 5 km] and [5 m/s, 5 m/s, 5 m/s], the accuracy of the star trackers is 3'', the accuracy of the horizon sensors is 0.05° [28], the sampling period of the stellar angle measurement is 10 s, and the observation period of the pulsars is 1000 s. Figure 12 shows the navigation results of the proposed deeply integrated navigation method and the stellar angle navigation method. Initially, the position error of the stellar angle navigation method is lower than that of the deeply integrated navigation method. With the increase in the navigation period, the position error of the deeply integrated navigation method gradually decreases. When the navigation period is about 9 days, the position error of the stellar angle navigation method exceeds that of the deeply integrated navigation method. When the navigation period is 20 days, the position error is reduced by about 42.3%.

Conclusions
In this paper, a stellar angle measurement-aided pulse phase estimation method is proposed that can significantly reduce the computational cost of the estimation of the pulse phase. In this pulse phase estimation method, the time-varying Doppler frequency caused by orbital motion is roughly removed by introducing stellar angle measurement into the pulse phase estimation. Consequently, the pulse phase can be estimated by the epoch folding method. In addition, we proposed a deeply integrated navigation method using the X-ray pulsar and the stellar angle, which is suitable for HEO spacecraft. Some simulations were conducted for three Earth-orbiting spacecraft. The simulation results show that the stellar angle measurement-aided pulse phase estimation method can process millisecond pulsar signals with a significantly lower computational cost than the current pulse phase estimation method. Compared with the stellar angle navigation method, the position error for the HEO spacecraft of the proposed navigation method is reduced by about 42.3%.