On-Orbit Pulse Phase Estimation Based on CE-Adam Algorithm

: Pulse phase is the basic measurements of X-ray pulsar-based navigation, and thus how to estimate a pulse phase for an orbiting spacecraft is important. The current methods for on-orbit pulse phase estimation could provide an accurate estimation performance enhancing with the photon amount, but its central processing unit (CPU) time cost also increases sharply with the increase of photon amount. In this paper, an on-orbit pulse phase estimation method based on the cross-entropy adaptive moment estimation (CE-Adam) algorithm is proposed to reduce the CPU time cost while retaining decent estimation accuracy. This method combines the CE and Adam algorithms, and is able to obtain a global optimum with low CPU time cost. The performance of the proposed algorithm is veriﬁed by simulation data and real data from the Neutron Star Internal Composition Detector (NICER). The results show that the proposed algorithm could greatly reduce the CPU time cost, which is about 1.5% of the CE algorithm, and retain similar estimation accuracy of pulse phase with CE algorithm.


Introduction
X-ray pulsar-based navigation (XPNAV) is a developing spacecraft autonomous navigation technique [1]. After being first introduced in 1981 [2], the past 40 years have witnessed a significant growth in XPNAV, including the pulse phase estimation as well as the navigation algorithm. Moreover, some flight experiments were performed such as SEXTANT (Station Explorer for X-ray Timing and Navigation Technology) [3] by the United States and the TG-2 (Tiangong-2) spacelab [4] and the Insight-HXMT (Insight-Hard X-ray Modulation Telescope) Satellite [5] of China.
Since the pulsar signal is extremely weak, an orbiting spacecraft can only record a series of photon TOAs (times of arrival) rather than a continuous pulsed signal [1]. Assuming the spacecraft is stationary or performs a uniform linear motion towards the pulsar, there are two types of method to estimate the pulse phase: the epoch-folding method [6] and the maximum likelihood estimator (MLE) method [7]. The epoch-folding method estimates pulse TOA by comparing the template with the empirical profile recovered by photon TOAs. Many methods based on epoch folding have been proposed in recent years. A wavelet-bispectrum algorithm is proposed in reference [8] to reduce the central processing unit (CPU) time cost of pulse phase estimation. In this method, wavelet transform is used to decompose the empirical profile and obtain low-frequency components of the empirical profile. Then, the low-frequency components are used to estimate the pulse phase [8]. In 2020, reference [9] proposed a genetic algorithm (GA)-optimized EMD (empirical mode decomposition)-CS (compressed sensing) with high accuracy and small computational load. In this method, the template is decomposed by the EMD (empirical mode decomposition) to obtain IMFs (intrinsic mode functions), and the IMFs selected by GA constitute the optimized measurement matrix. Then the pulse phase is estimated using the optimized measurement matrix [9]. In order to solve the sensitivity for noise, reference [10] proposed a variable step size iteration method which estimates pulse phase using a one-dimensional slice of template and empirical profile. The epoch-folding method is simple to implement, but needs lots of photons for the recovery of empirical profile and has an estimation accuracy less than MLE. However, the CPU time cost of MLE increases sharply as photon amount increases [7].
For orbiting spacecraft, the phase estimation problem become more complicated in that the orbit motion of a spacecraft introduces an unknown Doppler frequency into the pulse signal. To cope with this, Golshan et al. proposed a phase-tracking algorithm that divides the whole observation period into several short intervals [11]. During each interval, the spacecraft is approximated to perform a linear uniform motion, and in this case, the Doppler frequency is assumed to be constant and is tracked by a digital phase-locked loop (DPLL) [11]. In 2013, Huang et al. modified the DPLL to be a two-order Kalman filter and improve the performance of phase tracking [12]. In 2015, the phase tracking method was verified by the lab equipment of the Neutron Star Internal Composition Detector (NICER) and SEXTANT teams [13]. The phase tracking method works well for young pulsars such as PSR B0531+21(Crab) but fails when applied to faint millisecond pulsars such as PSR B1821-24 and PSR B1937+21. (The flux of Crab pulsar is 3667 photons/m 2 /s, but that of PSR B1937+21 is 0.161 photons/m 2 /s [14].) This is because the threshold for faint pulsars to obtain a reliable result is over 100 s, during which the assumption that a spacecraft performing a uniform linear motion is usually violated when the spacecraft is orbiting in a low Earth orbit. To overcome the drawbacks of the phase-tracking method, Wang and Zheng proposed a phase estimation method with the aid of orbital dynamics of spacecraft [15,16]. This method derives a linearized pulse phase propagation model, which modifies the pulse phase propagation model as a term of the position of the spacecraft. In addition, this method estimates the parameters of the linearized model to eliminate the impact of Doppler frequency. Compared with the phase-tracking algorithm, the phase estimation method with the aid of orbital dynamics of spacecraft does not divide the observation into intervals and thus is feasible for millisecond pulsar. Wang and Zheng recommend estimating the parameters by the direct search method [17]. Unfortunately, when using the direct search method, the CPU time cost is still high when the photon amount is large. In addition, Xue proposed estimating the Doppler frequency and pulse phase by a maximum likelihood estimator simultaneously [18], and this method also incurs high CPU time cost. Also, the photon background is an important factor which will impact on the performance of pulse phase estimation. Reference [19] shows that the photon background constantly changes because of orbital motion of the spacecraft. Reference [20] argues the effect of orbital motion on photon background could be reduced by data screening and the photon background could be assumed as a constant after data screening. Thus, in the existing methods, the photon background is assumed as a constant in the observation period.
For an on-orbit pulse phase estimation problem, there is a contradiction between estimation accuracy and CPU time cost. In order to estimate an accurate pulse phase, large amount of photons is needed, which leads to a high CPU time cost [21]. In reference [3], the error of position of XNAV is less than 10 km. In order to achieve this goal, the accuracy of phase estimation for Crab pulsar and PSR B1821-24 should be about 1 × 10 −3 cycle. The linearized method could provide an accurate pulse phase estimation result with a high CPU time cost. According to the simulation results, the linearized method could obtain the pulse phase with enough accuracy when observation period is 1000 s. However, based on the computation environment of this paper, the CPU time cost of the linearized method is about 8000 s. This means that we can obtain the pulse phase of the first observation period after 8 observation periods, which is unacceptable for a navigation system. Thus, the CPU time of the pulse phase estimation method must shorter than the observation period. In addition, the position of spacecraft should be estimated using the estimated pulse phase. However, the CPU time cost by the pulse phase estimation might hinder the real-time navigation process. This is because the spacecraft keeps orbiting during the CPU time and the final estimated pulse phase reflects the position of the spacecraft before the CPU time. Thus, the CPU time cost by the pulse phase estimation should be as short as possible.
To cope with this problem, we aimed to balance the pulse phase estimation accuracy and CPU time cost, which demands a highly efficient computational method with a decent estimation performance.
We cast the pulse phase estimation problem into a multi-extremal optimization problem. The cross-entropy (CE) algorithm is a classical method to solve a multi-extremal optimization problem, but has a heavy computational burden [22][23][24]. The Adam (adaptive moment estimation) algorithm is of high computational efficiency, but could only obtain a local optimum [25,26]. To combine the advantages of the two methods, this paper proposes a CE-Adam algorithm which features a low CPU time cost and global optimum. We compare the performance of the CE-Adam algorithm with the CE algorithm, DE (differential evolution) algorithm [27] and PSO (particle swarm optimization) [28] and find that the CE-Adam algorithm could provide accurate estimation results of pulse phase and significantly reduce the CPU time cost. In addition, we consider that the proposed algorithm could also be used to detect period signal in photon sequences over larger observational periods, which might indicate new physical phenomena.
The organization of this paper is as follows. Section 2 describes the principle of pulse phase estimation. Section 3 briefly reviews the CE algorithm and Adam algorithm. Section 4 shows the proposed CE-Adam algorithm. Simulation data and the real data are processed in Section 5 to verify the proposed algorithm. Finally, the summary of the conclusions is given in Section 6.

Principle of Pulse Phase Estimation
The arrival event of X-ray photons {t k } N k=1 follows the non-homogeneous Poisson process. Within the interval (t a , t b ), the probability of k photons arrival is [29]: where Λ satisfies Λ(t) = t 0 λ(s)ds. λ is the rate function and can be given by: The parameters α, β represent average signal and total background count rates in units of counts per second, respectively. h is the normalized profile function of pulsar, φ(t) represents the pulsar signal phase at the detector. φ(t) can be expressed as [14]: where φ 0 is the phase evolution at the reference observatory (we usually choose the Solar System Barycenter (SSB) as the reference observatory), and τ(t) represents the light propagation time from the detector to the reference observatory.
For an orbiting spacecraft, φ(t) can be expressed as [14]: where x(t) is the spacecraft position vector, n is the pulsar direction vector, x (t) is the predicted position vector of spacecraft, δx(t) = x(t) − x(t), e(t) can be linearized as e(t) = q + f (t − t a ), q is the initial phase at time t a , and f = .
φ is the frequency of the pulsar signal.
According to Equation (4), φ(t k ) can be expressed as: The joint probability density function of the pulsar photon arrival event is: We define Equation (5) as the likelihood function. The natural logarithm of the likelihood function is: When the observation period is long enough, the second term in Equation (6) is a constant [29]. Thus, Equation (6) can be transformed into: Parameters q and f can be estimated by solving the following optimization problem [14]: The CRLB (Cramer-Rao lower bound) for estimation for q is given by [21]: where T obs is the observation period and I p is a constant for each pulsar. It can be seen that the estimation accuracy of q is increased with observation period. As a result, large amounts of photons are needed for an accurate estimation of pulse phase. Equation (9) shows CPU time cost enhancing with the amount of photons. Thus, there is a contradiction between estimation accuracy and CPU time cost, and a computationally efficient algorithm with a decent estimation accuracy is needed. As shown in Figure 1, the objective function shown as (8) has many extrema. Then, the gradient-descent technique could only arrive at a local optimum. Thus, a global optimization algorithm is needed.

Cross-Entropy (CE) Algorithm
Assume S(x) is an objective function over all elements x in χ, and the maximum of it is γ * . We have [22]: The CE method defines a family of probability density functions { f (·; v), v ∈ V } on χ. Then, the CE method is associated with (11) the estimation of [22]: If γ is close to γ * , {S(x) γ} is a rare event. The CE method solves the maximization problem by making adaptive changes to the pdf according to the Kullback-Leibler and creating a sequence f (·; u), f (·; v 1 ), f (·; v 2 ), . . .. By finding the theoretically optimal density f (·; v * ), the tuples {(γ t , v t )}, which are generated by the CE method gradually converge to the small neighborhood of the optimal tuple (γ * , v * ) [22]. The steps of the CE method solving the continuous multi-extremal optimization problems are as follows [22].

1.
Choseμ 0 ,σ 2 0 , the elite sampling rate ρ, and the sample size N. Let the number of iterations t and the elite sample size N e be: 2.

3.
Choose the best N e samples as the elite samples.

Adam Algorithm
The Adam algorithm is a stochastic gradient descent (SGD) algorithm, which only calculates the first-order gradients. The expression of the Adam algorithm is as follows [25]: where θ = (q, f ), η represents the step size, g n is the gradient of the objective function at θ n , β 1 and β 2 are exponential decay rates for the moment estimates. m n and v n are the moment vectors, is a very small constant.

Recursive CE Algorithm
The CE algorithm has a good global optimization performance. However, the objective function values of all photons need to be calculated during each iteration, which leads to a large CPU time cost. A recursive CE algorithm is proposed to reduce the CPU time cost. The recursive CE algorithm divides all N photons into subsets [N 1 , N 2 , . . . , N n ] according to the observation period and selects a subset in each iteration to calculate the objective function and update the parameters. In other words, the improvement of the recursive CE algorithm to the CE algorithm is that the objective function is changed to S t (•). The expression of S t (•) is: where t is the number of iterations, m l is the capacity of N l , and N represents all photons, N = [N 1 , N 2 , . . . , N n ], N 1 , N 2 , . . . , N n are the subsets of photons, and N l is the currently used subset.
The µ k and σ k are calculated from Equation (14) and are used to update theμ t andσ t by Equation (15).

Proposal of CE-Adam Algorithm
The recursive CE algorithm uses few photons in each iteration, which reduce the CPU time cost as well as the estimation accuracy of q and f. In order to overcome this disadvantage, the parameter updated by the recursive CE algorithm is used as the initial value of the Adam algorithm. Then, the final value of the Adam algorithm is used to update the distribution parameters of the recursive CE algorithm.
where θ Adam−0 and θ Adam−end respectively represent the initial value and final value of Adam algorithm. The above algorithm is called the CE-Adam algorithm, and its procedure can be described as follows. Firstly, the recursive CE algorithm is adopted. After each recursion, the Adam algorithm is applied with the updated parameter's mean value as the initial value to obtain better than expected values of the parameters and update the distribution parameters of the recursive CE algorithm. The flow chart of the CE-Adam algorithm is shown in Figure 2.

Experiments
The simulation data of PSR B1821-24 and Crab pulsar and real data from NICER are used to verify the performance of the proposed CE-Adam algorithm. The algorithm is utilized to estimate q (initial phase) and f (frequency) in Equation (5). For simulation data, we define the evaluation criterion of the estimation performance is the RMS (root mean square) error of q. For the real data, the estimation performance is assessed by the error of estimated q relative to that obtained by the Global Positioning System (GPS). The computational cost is accessed by the CPU time cost. The computation environment contains Intel Core i5-7500@3.4GHZ (Intel, Santa Clara, CA, US), memory of 8G and python 3.8.

Experiments with Simulation Data
We chose the PSR B1821-24 pulsar, with the simulation parameters listed in Table 1, as the observed pulsar [15], and the profile of it is given in Figure 3. The results are all calculated from 1000 Monte Carlo simulations.  Assuming that the observation period is 1000 s, the photon TOAs recorded by an orbiting spacecraft are simulated. The objective function for estimating q and f is Equation (5). We compare the estimation performance of the CE algorithm, recursive CE algorithm, and CE-Adam algorithm. The average estimation error of q from 1000 Monte Carlo simulations is shown in Figure 4. It can be seen that the estimation accuracy of q for the three algorithms increases with the iterations. The estimation error of q for the recursive CE algorithm is higher than the other algorithms. This is because that the recursive CE algorithm uses a few photons in each iteration, which reduces the estimation accuracy of q. The estimation result for the CE-Adam and the estimation result for the CE algorithm are close to each other. This is because the CE-Adam algorithm overcomes the problem of the recursive CE algorithm and greatly improves the estimation accuracy.  Figure 5 shows the estimation error of q obtained by the four algorithms. It can be seen that the estimation accuracy of q for the four algorithms increases with the increase of the observation period. When the observation period exceeds 1000 s, the estimation error of q of the DE algorithm is higher than PSO, CE and CE-Adam. The estimation error of q for PSO, CE and CE-Adam are close to each other.   Table 2 shows the estimation error of q and f. It can be seen that the PSO, CE and CE-Adam perform a lower estimation error than DE. Since the pulse phase estimation is mainly focused on the pulse phase, the estimation of f will not be analyzed in the rest of this paper.   Figure 7 shows the CPU time cost by the four algorithms. Although the CPU time of CE-Adam algorithm is not the shortest when the observation period is less than 1000 s, its advantage gradually manifests when the observation period overtakes 1000 s. When the observation period reaches 3000 s, the CPU time cost of the CE-Adam algorithm is much lower than the others. For faint millisecond pulsar, in order to obtain an accurate estimation result, the observation period of more than 1000 s is needed. Therefore, compared with other algorithms, the algorithm proposed in this paper can reduce the CPU time cost while ensuring the accuracy of pulse phase estimation. We assume that the observation period is 1000 s and compare the linearized method with the proposed method. Table 3 shows the estimation error of q and CPU time cost. It can be seen that the proposed method could obtain similar estimation error of q with shorter CPU time. In addition, we investigated the estimation performance of the CE-Adam algorithm for Crab pulsar. Compared with PSR B1821-24, Crab pulsar has a lager flow rate, resulting in much higher CPU time cost. Table 4 shows the simulation parameters of the pulsar [16], and Figure 8 shows its profile. Assuming the observation period is 500 s, the photon TOAs recorded by the orbiting spacecraft are simulated. The objective function of estimating q and f is Equation (5). We compare the estimation performance of CE, DE, PSO, and the CE-Adam algorithm proposed in this paper.  Figure 9 shows the estimation error of q and CPU time for the four algorithms. It can be seen that the estimation accuracy of q for the four algorithms is not very different, but the CPU time cost of the CE-Adam algorithm is much lower than the other algorithms. Compared with the CE algorithm, the CPU time cost is about 1.5%; and compared with DE and PSO, the CPU time cost is about 2.8%.

Experiments with Real Data
NICER is the International Space Station (ISS) payload, which was launched in June 2017. It is mainly used to observe pulsars and neutron stars. NICER has a large effective area (>1800 cm 2 ), wide band-pass (0.2-12.0 keV) and high time resolution (better than 300 ns) [30].
We use the real data of the Crab pulsar on 26 December 2018 (Obs (observation)_ID: 1013010147). The predicted orbit is used to correct the barycenter of photon TOAs, and the CE-Adam algorithm proposed in this paper is used for on-orbit pulse phase estimation. The parameters of the predicted orbit are listed in Table 5.
We selected eight groups of data with observation period longer than 2000 s, which are shown in Table 6. The estimation performance was assessed by the error of estimated q relative to that obtained by GPS. Figure 10 shows CPU time and the estimation error of q for the CE-Adam algorithm. The average CPU time is 63.1 s, and the average error of estimation of q is 1.06 × 10 −3 .
Finding that the CPU time and q estimation error of dataset 1 are both exceptionally low, we used CE-Adam algorithm to process dataset 1 for 100 times. The results of estimation error of q and CPU time are shown in Figure 11. It can be seen that the estimation error of q and CPU time are uncertain. Most results of estimation error of q are about 1 × 10 −3 , and most results of CPU time are about 30 s. Although the RMS error of q is satisfactory, the uncertainty of the CE-Adam method could impact the estimation of q. Thus, further study is needed to overcome the uncertainty of the CE-Adam method.

Conclusions
In order to solve the problem in on-orbit pulse phase estimation that the CPU time cost increases sharply with the observation period, a pulse phase algorithm based on the CE-Adam algorithm, which has low CPU time cost and retains decent estimation accuracy, was proposed. Simulation results showed that the proposed algorithm had a lower CPU time cost than other algorithms. The on-orbit pulse phase estimation results obtained with real data of NICER showed that the algorithm proposed in this paper can quickly and accurately process the real data, which has practical engineering significance for pulsar navigation. In addition, we found that the uncertainty of the CE-Adam algorithm would impact the CPU time and accuracy of pulse phase estimation, which should be overcome in a further study.