A fast pulse phase estimation method for X-ray pulsar signals based on epoch folding

Xue Mengfan,Li Xiaoping,Sun Haifeng,Fang Haiyan

School of Aerospace Science and Technology,Xidian University,Xi'an 710126,China

A fast pulse phase estimation method for X-ray pulsar signals based on epoch folding

Xue Mengfan,Li Xiaoping*,Sun Haifeng,Fang Haiyan

School of Aerospace Science and Technology,Xidian University,Xi'an 710126,China

X-ray pulsar-based navigation(XPNAV)is an attractive method for autonomous deepspace navigation in the future.The pulse phase estimation is a key task in XPNAV and its accuracy directly determines the navigation accuracy.State-of-the-art pulse phase estimation techniques either suffer from poor estimation accuracy,or involve the maximization of generally nonconvex object function,thus resulting in a large computational cost.In this paper,a fast pulse phase estimation method based on epoch folding is presented.The statistical properties of the observed profile obtained through epoch folding are developed.Based on this,we recognize the joint probability distribution of the observed profile as the likelihood function and utilize a fast Fourier transform-based procedure to estimate the pulse phase.Computational complexity of the proposed estimator is analyzed as well.Experimental results show that the proposed estimator significantly outperforms the currently used cross-correlation(CC)and nonlinear least squares(NLS)estimators,while significantly reduces the computational complexity compared with NLS and maximum likelihood(ML)estimators.

1.Introduction

Pulsars are highly magnetized,rapidly rotating neutron stars emitting uniquely identifiable signals that are periodical and predictable,throughout the electromagnetic spectrum with periods ranging from milliseconds to thousands of seconds.

The repetition period of the radiation signals is simply the rotation period of the neutron star.For some pulsars,the stability of their rotation periods over long timescales is as precise as an atomic clock.1–3Of all pulsars,the ones which are visible in the X-ray band of the electromagnetic spectrum are called''X-ray pulsarquot;.3,4Compared with the other types of pulsars,the X-ray pulsars are more suitable for use in deep space navigation because of the existence of small size X-ray detectors that can be mounted on a spacecraft.5

X-ray pulsar-based navigation(XPNAV)is a developing celestial navigation method and receives increasing attention.It is promising to fulfill completely autonomous navigation to reduce the dependence of current navigation system to ground-based operations,or to argument the current systems by employing additional measurements to improve their performance.4–7The concept of employing X-ray pulsars to estimate the position of deep space explorer was first proposed in 1974 and has grew rapidly during the last 40 years.8United States and the European Space Agency have analyzed the feasibility of XPNAV and continuously studied on the subject.9–12

In the recent years,many researchers have investigated different applications of NPNAV,including both absolute and relative navigations.13–15It has been shown that one key issue of XPNAV is how to precisely estimate the phase delay between the observed profile and the predicted one.13–16To date,many pulse phase estimation algorithms have been developed for XPNAV.The maximum likelihood(ML)phase estimator presented in Ref.13directly utilizes the detected photon time of arrivals(TOAs);its estimation accuracy approaches the Cramér-Rao lower bound(CRLB)as the observation duration increases,but the direct employment of the photon TOAs makes the amount of calculation and storage rapidly increase with the growth of received photons and the direct search of the ML solution also leads to a high computational cost.The Fourier transform-based pulse delay estimation method given by Taylor has the advantage that its accuracy is independent of the bin size,but it inevitably involves a straightforward iterative procedure which is computationally intensive.17Emadzadeh has proposed two different pulse delay estimators based on the observed profile obtained through epoch folding.4,13One uses the cross-correlation(CC)function between the photon intensities and greatly reduces the computational cost compared with the ML estimator.The other is a nonlinear least squares(NLS)estimator.However,both estimators suffer from poor estimation accuracy and are not asymptotically efficient,and the NLS estimator still involves the maximization of a non-convex objective function,thus resulting in a direct search-based procedure.In Ref.18,the authors describe an approximate ML estimator at low signal to noise ratio(SNR)values.The appropriate approximation of the likelihood function reduces the computational complexity,but results in a degradation of the estimation accuracy as the SNR values increase due to the deviation between the established statistical model of the observed profile and the exact situation.In Ref.19,the authors recast the problem of pulsar phase estimation as a cyclic shift parameter estimation problem under multinomial distributed observations and develop a fast near-maximum likelihood phase estimation method based on this.This strategy is essentially based on the conditional probability density function of the photon TOA and actually complicates the phase estimation problem of X-ray pulsar,which can be directly mathematically formulated by the epoch folding procedure as presented in the subsequent content of this paper.

In this paper,the X-ray pulsar radiation is characterized as a non-homogeneous Poisson process(NHPP)and the observed profile obtained through the epoch folding procedure is statistically modeled as a heteroscedastic Poisson sequence.Upon this model,we recognize the joint probability distribution of the observed profile as the likelihood function and employ a computationally efficient fast Fourier transform(FFT)based procedure to estimate the pulse phase of X-ray pulsar signals.The rest of this paper is organized as follows.In Section 2,the mathematical models describing the X-ray pulsar signals are presented.Based on this,the heteroscedastic Poisson model formulating the observed profile is established.

Section 3 explains how the pulse phase can be estimated by a FFT based procedure by employing the joint probability distribution of the observed profile as the likelihood function.Computational complexity of the proposed estimator is studied in Section 4.In Section 5,experiments are carried out to evaluate this new technique's performance,using both simulated data and real data.Finally,Section 6 concludes the study.

2.Heteroscedastic Poisson model of X-ray pulsar observed profile

The original measurement of X-ray pulsar is the TOAs of all the X-ray photons from the pulsar source as well as the background.TOA of a photon is recorded by the X-ray detector when the photon hits the detecting material.20–22In XPNAV,a low power X-ray detection system capable of measuring the photon TOAs with submicro second accuracy is required.The detector must have a low background noise,a large detection area and a light weight.The Naval Research Laboratory(NRL),as part of the Defense Advanced Research Projects Agency(DARPA),has developed a new X-ray silicon-based detector that satisfies all the above-mentioned requirements.19

To obtain the observed profile,the measured photon TOAs are first transformed to the solar system barycenter(SSB)and then assembled into a single pulse cycle through the procedure of epoch folding.3–5,13In this paper,to focus our discussion,we assume that the photon TOAs have been transformed to SSB.In what follows,according to the presented X-ray detection method,mathematical equations are used to describe the photon TOAs at the SSB;then upon this,statistical properties of the observed profile are given.

2.1.Mathematical model of X-ray pulsar signals

Let Ntbe the number of arrival photons at the time interval(0,t).The counting process{N(t),t≥0}can be modeled by a NHPP with a time-varying intensity λ(t) ≥ 0.23,24For a fixed time interval(ts,te),the number of arrival photons Nts,teis a Poisson random variable with parameterIts distribution low is

and its mean and variance are

Furthermore,since{N(t),t≥0}has independent increments,the numbers of detected photons in any non overlapping time intervals are independent from each other.The intensity function λ(t)whose unit is ph/s includes all the arriving photons from the X-ray pulsar and the background.It is expressed as

where h(φ)is the normalized pulsar profile,φ(t)represents the evolution of the pulse phase with respect to the time t as seen at the SSB,and α and β are the known effective background and source photon arrival rates,respectively.The pulsar profile h(φ),which is unique to a particular pulsar and defines the non-homogeneous nature of the source photons,is defined at the phase interval [0,1) and normalized to satisfyandand then its definition domain is extended beyond[0,1)by letting h(φ+n)=h(φ) ?n∈ Z.23,24Fig.1 shows the profile of the Crab pulsar(PSR B0531+21)in the X-ray band(2–16 keV),created by using about 1000 s real data observed by the rossi X-ray timing explorer(RXTE).

The pulse phase φ(t)is the sum of the initial phase φ0at a reference time t0and the accumulation of phase since t0.It can be expressed as the Taylor series

where f and˙f are the X-ray source frequency and its first derivative valid at time t0,respectively,and O(m)is the highorder item and is expressed as

with f(m)the mth derivation of the source frequency at time t0.23The parametersand f(m)of a pulsar are obtained through pulsar timing techniques and they need to be constantly revised to minimize the pulsar timing residuals which would result in additional uncertainties of the pulse phase.For example,the above mentioned Crab pulsar is the brightest rotation-powered pulsar in the X-ray band,yielding about 9.9X10-9ergs/(cm2/s)of X-ray energy flux in the 2–10 keV band,which is a significant advantage for its application in spacecraft navigation,3but it is also a young pulsar and its timing residuals exhibit significant timing noise,so its phase model needs to be frequently revised to minimize the timing noise and optimal extrapolation of its timing residuals is essential for it to be used in navigation.25,26

2.2.Epoch folding

The process of epoch folding is to recover the observed profile from the photon TOAs.3–5,13,23,24It is implemented as follows:All the time tags within an observation are folded back into one phase cycle of the source by calculating their corresponding phase values according to the phase model of Eq.(4)and then only retaining the fractional part of each phase value,which is referred to as the normalized phase.Afterwards,divide the phase cycle(0,1)into some equal-length bins and drop each of the normalized phase values into the appropriate phase bin.The resulting histogram over one phase cycle is the so-called observed profile of X-ray pulsar.The overall epoch folding procedure is illustrated in Fig.2.

Fig.1 Profile of Crab pulsar obtained by RXTE.

Fig.2 Procedure of epoch folding.

Assume that the observation duration is Tobs,one phase cycle contains Nbbins,and cidenotes the photon count of the ith bin after the epoch folding procedure.According to Eqs.(1)and(2),ciis a Poisson random variable with the following mean and variance:

where

φiis the center phase of the ith bin andIt is shown that the variance of cichanges with the phase.Therefore,the observed profile,denoted bycan be mathematically formulated by a heteroscedastic Poisson model.Moreover,according to the independent increment property of the Poisson process,ciis statistically uncorrelated,namely,

3.Fast maximum likelihood phase estimator

Additionally,since ciis statistically uncorrelated,the joint probability distribution of the observed profile is multinomial and is given by

Recognizing the joint probability distribution given in Eq.(10)as the likelihood function,an ML estimator based on the observed profile is provided by maximizing Eq.(10)with respect to the unknown parameter φ0.Equivalently,the loglikeli hood function(LLF)can be maximized.It is expressed as

where

It is obvious that the term A does not depend on the unknown parameter φ0.Furthermore,since λp(φ)is a periodic function,its integral sum over one period is not a function of the initial phase φ0,which means that the term B is also independent of φ0.Therefore,the term A and B can be dropped from the LLF.We end up with the following LLF:

The initial phase can be estimated by solving the following optimization problem:

Since the statistic λ(φ0)is a periodic function,the inner multiplication on the right-hand side of Eq.(16)can be recognized as the cyclic CC function between the observed profile^λ and the logarithmic transformation of the actual radiation intensity as a function of φ0.

In order to exploit the cyclic convolution property of the FFT for the evaluation of the global maximum of Eq.(16),we discretize the range(0,1]of the parameter φ0into Nbintervals of width 1/Nb.The index kφ,corresponding to a coarse estimateof φ0with a resolution of 1/Nb,namely,can be obtained from(16)rewritten as27

where D is the following orthogonal unit matrix27:

Since the resolution of the coarse estimateprovided by the discrete cross correlation is definitely limited by the value of Nb,a fine estimate^φ0is attained by means of a parabolic interpolation of the objective sequencefaround its maximum f(kφ).

Moreover,it is interesting to notice that the estimation rule Eq.(15)is close to the CC estimator,4from which it differs mainly by the employment of the logarithm of λp(φ)instead of λp(φ).Actually,since the logarithmic transformation is a kind of variance stabilizing transformations,28it can,to some extent,eliminate the influence of the heteroscedastic property of the observed profile described in Section 2.2 on the phase estimation accuracy.Therefore,our proposed estimator may have a higher precision compared with the CC estimator.This is also validated by the experimental results presented in Section 5.

4.Computational complexity analysis

In this section,the computational complexity of the proposed pulse phase estimator is studied and is compared with that of the currently used estimators.Since when considering the change in the source frequency,the procedure of transforming the photon TOAs into photon phases based on the phase model Eq.(4)is needed for all the phase estimators,in the following comparison analysis we assume the phase values of all the arrival photons have been known for the sake of simplicity.

To use the proposed estimator,the observed profile is needed.To obtain it,Nb(Nc-1)real additions must be performed,with Ncbeing the number of phase cycles within the observation time.Since the computed photon counts need not to be normalized in our epoch folding procedure,the 2Nbdivisions needed in the epoch folding procedure given in Refs.23,24can be eliminated.Since the Nbreal logarithms needed to construct the cost function can be pre-calculated,we do not contain this part of calculation into the computational cost of the proposed estimator.To calculate the cyclic cross-correlation sequence in the cost function Eq.(17),two Nb-point FFT and a Nb-point inverse FFT are involved,which requires a total number of 2Nblog2Nbreal multiplications and 3Nblog2Nbreal additions,assuming that one complex addition contains two real additions and one complex multiplication contains four real multiplications and two real additions.Hence,the total operations required by the proposed estimator are approximately 3Nblog2Nb+Nb(Nc-1)real additions and 2Nblog2Nbreal multiplications.We observe that the amount of necessary additions is a linear function of Ncand grows linearly as the observation time increases,and the required real multiplications only depend on the number of bins in one phase cycle,Nb.

We also give the computational complexities of state of the art phase estimators to perform a contrastive analysis.The ML estimator derived in Ref.13needs a total number of 5NphNgreal additions,2NphNgreal multiplications and NphNgreallogarithms,where Nphdenotes the total number of detected photons within the observation duration.Note that the amount of calculations is proportional to the number of received photons and it will significantly increase as the observation time becomes longer,which is a noticeable defect of the ML method.The NLS estimator presented in Refs.13,17,24requires Nb(Nc-1)real additions and 2Nbreal divisions for the epoch folding procedure and NbNgreal additions,NbNgreal multiplications and NbNgreal subtractions for the grid search,with Ngbeing the number of grid points considered.The CC estimator developed in Ref.4requires the same operations as the NLS estimator to implement the epoch folding procedure and the same operations as the proposed estimator to calculate the cyclic CC sequence.

Table 1 Computational costs of different estimators.

A comparison of the complexities of all the methods mentioned above is summarized in Table 1.It is worth noting that since Nph?Nbin most cases,the computational cost of the ML is the highest.Taking the Crab pulsar as an example,when the detector area is 200 cm2,approximately 9X105photons could be detected within an observation duration of 30 s,which corresponds to about898phase cycles.Setting Nb=Ng=1024,the total numbers of arithmetic operations needed by the ML,NLS and CC methods are about 7.373X109,4.066X106and 9.717X105,respectively,whereas the proposed estimator only requires about 9.696X105operations.

5.Experimental results and discussion

To evaluate the performance of the proposed pulse phase estimator,a series of experiments is carried out using both the simulated data obtained by the ground-based simulation system and the real data obtained by the RXTE satellite.We assess the performance of proposed estimator against the other three methods discussed in Section 4,including the ML estimator,the NLS estimator and the CC estimator.

5.1.Simulated photon TOAs

5.1.1.Ground-based simulation system of X-ray pulsar signals Fig.3 shows the principle diagram of the X-ray pulsar signal ground-based simulation system,which consists of three modules:The voltage source module,the optical attenuation module and the photon detection module.The voltage source module with high time–frequency stability is to synthesize a voltage signal of similar shape as the intensity function.It includes an atomic clock,a crystal oscillator,a frequency synthesizer and an arbitrary waveform generator.Under the premise of the physical process being consistent,we replace the X-ray source with the visible light source in order to reduce costs.The optical attenuation module generates the single photon sequence by simulating the practical attenuation process in the propagation of the X-ray signal.In this module,an electronic control visible light source and a variable optical attenuator are fixed inside an optical shielding cavity.The photon detection module,which includes a photomultiplier,an optical pulse discriminator and the electronic readout circuit,is designed to record the photon TOAs.

With the characteristics of high time resolution of 10 ns,high stability up to 10-9,low dark counts and small counting error,the simulation system can generate photon TOAs of any known pulsar.For detailed information,you can refer to our previous work.29It is noteworthy that what this system simulates is the photon TOAs transformed to the SSB,but not the photon TOAs directly detected by the detector onboard a moving spacecraft,and the system cannot simulate the change in the source frequency yet,meaning that in Eq.(4),the first and higher order derivatives of the source frequency are all set to zero.Since the simplified constant frequency model does not affect the validation of the pulse phase estimation algorithms,in the following numerical simulations,we use this simulation system to generate the photon TOAs of the selected pulsars.

5.1.2.Results and analyses

Two pulsars are selected:B0531+21(Crab)and B1509-58,whose simulation parameters are listed in Table 2.4,13The simulation results are obtained by using the Monte-Carlo technique,with 100 independent realizations of photon time tags generated by the ground-based simulation system.The initial phase φ0is set to be 0.3021.The number of bins in one pulse cycle,namely Nb,is set as 1024 for Crab and 512 for B1509-58.The number of grid points Ngof the NLS and ML estimators is set equal to Nbso as to ensure the same phase resolution as the CC and the proposed estimator.For a fair comparison,the procedure of parabolic interpolation of the cost function around its maximum is also employed in the ML,NLS and CC methods to obtain a fine estimate.We evaluate the estimation accuracy in terms of root mean square error(RMSE),which is defined as

Fig.3 Principle diagram of X-ray pulsar signal ground-based simulation system.

Table 2 Parameters of employed pulsars.

In the first scenario,giving a fixed SNR value,we test the estimators' performance for different observation times using Crab pulsar.In Fig.4,the RMSE of the four methods for each observation time is plotted against the square root of the CRLB calculated according to Ref.13,at a background arrival rate α=4.4 ph/(cm2/s),corresponding to SNR ≈ -10 dB.It can be seen that when the observation time is short,a relatively large deviation between the RMSE and the square root of the CRLB is introduced no matter what estimator is used.The reason is that for a short observation time,the realization of the photon TOAs is not enough to represent the intensity function precisely,leading to a cost function whose global optimum generally locates at a farther point to the true value,φ0.As time goes on,the RMSE of each method becomes more attached to the CRLB and approaches to zero,also indicating that the proposed estimator is consistent.Besides,simulation results show that the proposed estimator outperforms the NLS and CC estimators and keeps close to the ML estimator at almost all the observation times.This is due to the fact that the established heteroscedastic Poisson model can accurately formulate the observed profiles of any observation times.

Fig.4 RMSE of different estimators for Crab.

Fig.5 RMSE versus SNR for B1509-58 at Tobs=300 s.

As a second scenario,we evaluate the RMSE versus SNR.The pulsar B1509-58 is considered.We fix the observation time at 300 s and calculate the RMSE of each estimator at different SNR values varying from-23 dB to 32 dB.The observation time is set to be 300 s.Simulation results plotted in Fig.5 show that when the SNR is less than-5 dB,the proposed estimator performs only slightly better than the NLS and CC estimators,whereas as the SNR grows higher,the proposed estimator tends to have more and more obvious advantage than the NLS and CC estimators.When the SNR is higher than 5 dB,the NLS and CC estimators enter into the saturation region,with their RMSE values changeless with the increasing SNR values,while the proposed estimator still keeps tightly attached to the ML estimator and clearly outperforms the NLS and CC estimators.

Although the ML estimator has an advantage over the proposed estimator on the estimation accuracy,the operation time taken by the estimator should be taken into consideration.In Fig.6,the CPU time cost by MATLAB to implement the calculation for one Monte-Carlo realization is plotted as a function of the observation time.The processor utilized is an Intel 2.5 GHz dual core.It can be seen that the CPU time for the ML estimator is much bigger than the ones for the other estimators,especially the CC and the proposed estimators,and it linearly grows as the observation time becomes longer and grows significantly faster than the other three estimators.The NLS,CC and the proposed estimators all totally contain two procedures,namely the epoch folding and the construction of cost function,of which the CPU time for the former grows linearly with the observation time and the CPU time for the latter does not change with the observation time,leading to linearly growing in these three estimators' computational complexities with the observation time.Because the cost function construction time for the NLS estimator is bigger than the ones for the CC and the proposed estimators and the time costs of the epoch folding procedure are almost the same for the three estimators,both CC and the proposed estimators require much lower computational load than the NLS estimator,especially when the observation time is short,as is shown in Fig.6,which is consistent with the results discussed in Section 4.Furthermore,the proposed estimator costs slightly less time than the CC estimator due to the elimination of the 2Nbreal divisions in the epoch folding procedure.

Fig.6 CPU time of different estimators.

5.2.Real data from RXTE

We use real data of the Crab pulsar obtained with the proportional counter array(PCA)on board RXTE.The PCA,which is a low-energy(2–60 keV)detection instrument,consists of five xenon- filled proportional counter units(PCUs)with a total effective area of about 6500 cm2.Its time resolution is 1 μs.We perform the data extraction and the pro file generation using the RXTE standard data analysis software,FTOOLS 4.0.30The event mode data of the observation 96802-01-12-00 detected in MJD 55776 is selected.To identify periods of good data from the whole observation duration,the following selection criteria are applied.30

(1)The minutes since the peak of the last South Atlantic Abnormal passage are larger than 30.

(2)The pointing offset of the detector is less than 0.02°.

(3)The electron contamination of each PCU is smaller than 0.1.

(4)To avoid contamination due to Earth's X-ray bright limb,the elevation angle of the source above the spacecraft horizon must be greater than 10°.

After filtering out the band timespan,all the reserved time tags are first transformed to the SSB to eliminate the Doppler effect,in which process the various clock corrections,the Roemer delay and the relativistic effects including the Sun Shapiro and the Einstein delays are considered,30and then the transformed time tags are folded back into one phase cycle to generate the observed profile.The above two stepsare implemented by employing the 'Fasebin' and' Fbssum' operations,of which Nbis set to 1024.The observed profiles of the delayed version are generated by directly adding a delay of 9 ms,corresponding to an initial phase of 0.2674,into the corrected time tags before epoch folding.We take the observed profile obtained over the total left exposure time of about 1011 s as the reference profile to be compared with the delayed profile.The background and the source arrival rates are about 8.813X103ph/s and 1.0216X103ph/s,respectively,corresponding to SNR≈-10 dB.

Fig.7 shows the phase estimation errors of the four methods for different observation times.The proposed estimator outperforms the NLS and CC estimators at most observation times,especially when the observation time is long,and its estimation accuracies keep close to the ones of the ML method at most observation times that are longer than 350 s,which is consistent with the simulation results given in Section 5.1.It is noteworthy that since the source arrival rate is large,when the observation time is long,the phase estimation accuracy is limited by the time resolution of the detected time tags and does not present a notable improvement,as can be seen in Fig.7.The fluctuation occurring between 350 s and 400 s may be caused by a kind of unknown strong noise during this period of observation.

Fig.7 Phase estimation errors of different estimators derived by using real data of Crab.

6.Conclusions

This paper has presented a new computationally efficient phase estimator for the X-ray pulsar signals.According to the statistical properties of X-ray pulsar signals, a heteroscedastic Poisson model is established to formulate the observed pulsar profile obtained through epoch folding.Based on this model,a new pulse phase estimator which employs the joint probability distribution of the observed profile as the likelihood function and can be calculated by a FFT-based procedure is derived.The performance of the proposed estimator is evaluated by the simulated photon TOAs as well as the real data captured by the RXTE.It is shown that when the observation time is short,our estimator has an estimation accuracy slightly superior to the ones of the NLS and CC estimators,while when the observation time is long,it performs signif icantly better than the NLS and CC estimators.Furthermore,its computational cost is much smaller than the ones of the ML and NLS estimators,making the estimator more suitable to be used on the resource-limited spacecraft.

Acknowledgement

This study was supported by the National Natural Science Foundation of China(No.61301173).

1.Matsakis DN,Taylor JH,Eubanks TM,Marshall T.A statistic for describing pulsar and clock stabilities.Astron Astrophys 1997;326(3):924–8.

2.Rots AH,Jahoda K,Lyne AG.Absolute timing of the Crab pulsar with ROSSI X-ray timing explorer.Astrophys J 2004;605(2):L129–32.

3.Sheikh SI.The use of variable celestial X-ray sources for spacecraft navigation dissertation.Maryland:University of Maryland;2005.

4.Emadzadeh AA.Relative navigation between two spacecraft using X-ray pulsars dissertation.Los Angeles:University of California;2009.

5.Liu J,Kang ZW,White P.Doppler/XNAV-integrated navigation system using small-area X-ray sensor.IET Radar Sonar Navig 2011;5(9):1010–7.

6.Feng DZ,Xu LP,Zhang H,Song SB.Determination of intersatellite relative position using X-ray pulsars.J Zhejiang Univ Sci C 2013;14(2):133–42.

7.Xue MF,Li XP,Fu LZ,Fang HY,Sun HF,Shen LR.X-ray pulsar-based navigation using pulse phase and Doppler frequency measurements.Sci China Inf Sci 2015;58(12)122202:1-10.

8.Fei BJ,Yao GZ,Du J,Sun WJ.The pulse profile and united measurement equation in XNAV.Sci China Phys Mech Astron 2011;40(5):644–50.

9.Xue MF,Li XP,Fu LZ,Liu XP,Sun HF,Shen LR.Denoising of X-ray pulsar observed profile in the undecimated wavelet domain.Acta Astronaut 2016;118:1–10.

10.Sheikh SI,Pines DJ,Ray PS,Wood KS,Lovellette MN,Wolff MT.Spacecraft navigation using X-ray pulsars.J Guid Control Dyn 2006;29(1):49–63.

11.Huang LW,Liang B,Zhang T,Zhang CH.Navigation using binary pulsars.Sci China Phys Mech Astron 2012;55(3):527–39.

12.Qiao L,Liu JY,Zheng GL,Xiong Z.Augmentation of XNAV system to an ultraviolet sensor-based satellite navigation system.IEEE J Sel Top Signal Process 2009;3(5):777–85.

13.Emadzadeh AA,Speyer JL.On modeling and pulse phase estimation of X-ray pulsars.IEEE Trans Signal Process 2010;58(9):4484–95.

14.Zhang H,Xu LP.An improved phase measurement method of integrated pulse profile for pulsar.Sci China Technol Sci 2011;54(9):2263–70.

15.Wang YD,Zheng W,Sun SM,Li L.X-ray pulsar-based navigation using time-differenced measurement.Aerosp Sci Technol 2014;36:27–35.

16.Liu J,Fang JC,Kang ZW,Wu J,Ning XL.Novel algorithm for X-ray pulsar navigation against Doppler effects.IEEE Trans Aerosp Electron Syst 2015;51(1):228–41.

17.Taylor JH.Pulsar timing and relativistic gravity.Philos Trans Phys Sci Eng 1992;341(1660):117–34.

18.Li JX,Ke XZ.Maximum-likelihood TOA estimation of X-ray pulsar signals on the basis of Poison model.Chin Astron Astrophys 2011;35(1):19–28.

19.Rinauro S,Colonnese S,Scarano G.Fast near-maximum likelihood phase estimation of X-ray pulsars.Signal Process 2013;93(1):326–31.

20.Wang YD,Zheng W,Sun SM.X-ray pulsar-based navigation system/sun measurement integrated navigation method for deep space explorer.Proc IMechE Part G:J Aerosp Eng 2015;229(10):1842–53.

21.Zhang H,Xu LP,Xie Q.Modeling and Doppler measurement of X-ray pulsar.Sci China Phys Mech Astron 2011;54(6):1068–76.

22.Wang YD,Zheng W,Sun SM.X-ray pulsar-based navigation system with errors in the planetary ephemerides for Earth-orbiting satellite.Adv Space Res 2013;51(12):2394–404.

23.Emadzadeh AA,Speyer JL.X-ray pulsar-based relative navigation using epoch folding.IEEE Trans Aerosp Electron Syst 2011;47(4):2317–28.

24.Emadzadeh AA,Jason LS.Relative navigation between two spacecraft using X-ray pulsar.IEEE Trans Control Syst Technol 2011;19(5):1021–35.

25.Deng XP,Coles W,Hobbs G,Keith MJ,Manchester RN,Shannon RM,et al.Optimal interpolation and prediction in pulsar timing.Mon Not R Astron Soc 2012;424(1):244–51.

26.Deng XP,Hobbs G,You XP,Li MT,Keith MJ,Shannon RM,et al.Interplanetary spacecraft navigation using pulsars.Adv Space Res 2013;52(9):1602–21.

27.Colonnese S,Rinauro S,Scarano G.Generalized method of moments estimation of location parameters:application to blind phase acquisition.IEEE Trans Signal Process 2010;58(9):4735–49.

28.Guan Y.Variance stabilizing transformations of Poisson,binomial and negative binomial distributions.Stat Probab Lett 2009;79(14):1621–9.

29.Sun HF,Xie K,Li XP,Fang HY,Liu XP,Fu LZ,et al.A simulation technique of X-ray pulsar signals with high timing stability.Acta Phys Sin 2013;62(10):109701 Chinese.

30.The ABC of RXTE:Contents Internet.[cited 2015 July 18].Available from: http://heasarc.gsfc.nasa.gov/docs/xte/abc/contents.html/.

Xue Mengfanis a Ph.D.candidate at School of Aerospace Science and Technology,Xidian University.She received the B.S.degree in 2011 from Shijiazhuang Tiedao University.Her research focuses on signal processing theory and X-ray pulsar-based navigation.

Li Xiaopingwas born in 1961.She received her Ph.D.degree from Xidian University in 2004.She is now a professor and Dean of the School of Aerospace Science and Technology of Xidian University.Her research focuses on information fusion,intelligent information processing and information security.

Sun Haifengwas born in 1986.He received his M.S.and Ph.D.degrees from Xidian University respectively in 2011 and 2015.His researches focus on celestial navigation and signal processing theory.

26 July 2015;revised 3 December 2015;accepted 1 February 2016

Available online 10 May 2016

Epoch folding;

Maximum likelihood;

Phase estimation;

X-ray pulsar;

X-ray pulsar-based navigation(XPNAV)

?2016 Production and hosting by Elsevier Ltd.on behalf of Chinese Society of Aeronautics and Astronautics.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

*Corresponding author.Tel.:+86 29 81891036.

E-mail address:xpli@xidian.edu.cn(X.Li).

Peer review under responsibility of Editorial Committee of CJA.