Novel X-ray Communication Based XNAV Augmentation Method Using X-ray Detectors

The further development of X-ray pulsar-based NAVigation (XNAV) is hindered by its lack of accuracy, so accuracy improvement has become a critical issue for XNAV. In this paper, an XNAV augmentation method which utilizes both pulsar observation and X-ray ranging observation for navigation filtering is proposed to deal with this issue. As a newly emerged concept, X-ray communication (XCOM) shows great potential in space exploration. X-ray ranging, derived from XCOM, could achieve high accuracy in range measurement, which could provide accurate information for XNAV. For the proposed method, the measurement models of pulsar observation and range measurement observation are established, and a Kalman filtering algorithm based on the observations and orbit dynamics is proposed to estimate the position and velocity of a spacecraft. A performance comparison of the proposed method with the traditional pulsar observation method is conducted by numerical experiments. Besides, the parameters that influence the performance of the proposed method, such as the pulsar observation time, the SNR of the ranging signal, etc., are analyzed and evaluated by numerical experiments.


Introduction
Being an autonomous navigation method applied for solar system and beyond, XNAV (X-ray pulsar-based NAVigation) has been a hot research area during the past ten years. Multiple aspects of XNAV, including Time of Arrival (TOA) measurement [1][2][3], time transfer [4][5][6], ambiguity resolution [7][8][9], filtering algorithms [10,11] etc., have been studied to verify the feasibility of XNAV. The theoretical positioning accuracy of XNAV is 10 m. Currently, the available accuracy of XNAV is about several hundred meters, which is far from the designed accuracy [7]. The gap is caused by various factors, including the inaccuracy of the noise model, the relativistic effects, the limitations of current X-ray detectors, the ephemeris error, etc. The limited available accuracy greatly hinders the application of XNAV, therefore, how to improve the positioning accuracy of XNAV under current conditions has become a critical and urgent issue.
The main solution of improving XNAV accuracy is to obtain extra information from the available navigation systems or multiple space sources. Reference [12] provided an integration of XNAV with a strap-down inertial system (SINS), which provides higher accuracy and requires shorter filtering period, but the SINS error accumulation would affect the performance of the integrated system, especially for long-duration missions. Another main method to augment XNAV is by observing stars, planets, or asteroids. In the solar system, the Sun is an obvious and readily available source for observation, which is used in XNAV augmentation. In Reference [13], a Sun sensor was utilized to observe the line-of-sight vector towards the Sun and the line-of-sight vector was integrated with X-ray pulsar observation to achieve autonomous navigation in Halo orbit. Reference [14] also utilized the line-of-sight vector of the Sun as the observation variable and adopted a residual orthogonal unscented Kalman filter (ROUKF) for navigation filtering. Planets are also common observation targets in the navigation of satellites that orbit around a central body. Accurate navigation is realized by combining the pulsar observation and the central body observation that is achieved by a star sensor [15,16]. This method of observing a central body has the obvious disadvantage of a limited application range. The method can only be used near the central body. For specific missions, asteroids are used to improve XNAV performance. For example, in Reference [17], the images of asteroids, captured by a navigation camera, were combined with pulsar observation to realize navigation in an interplanetary cruise. This method would have a great dependence on the quality of the asteroid images. Though utilizing the observation of space sources to improve XNAV has shown some potential, both the integration method using a Sun sensor and the augmentation method with central body observation or asteroid observation share the same flaw in that the star sensors have a limited accuracy and easily suffer from interference by noise. Thus, the improvements of XNAV based on the abovementioned methods are limited.
In this paper, we propose an augmentation method for XNAV based on X-ray communication (XCOM), which utilizes accurate X-ray ranging as the extra observation to improve the XNAV accuracy. XCOM is a newly emerged concept for deep space communication [18], which utilizes X-rays as the transmission medium. Being high energy light and not being affected by electrical and magnetic field [19], X-rays show fine space propagation performance. The idea of X-ray ranging is derived from XCOM. In our previous research, we have introduced X-ray ranging and presented a detailed performance analysis [20]. Compared with available augmentation methods, X-rays could provide accurate range measurements, which could serve as the observation variable in XNAV. In other words, by both X-ray pulsar observation and X-ray ranging measurement, the position and velocity of a spacecraft could be estimated by a Kalman filtering algorithm. Actually, there is much more potential in combining XNAV with XCOM. Besides accurate range measurement, XCOM-based simultaneous communication and ranging could also provide communication services [21]. As is known to all, the phase revolution model is of great importance in estimating TOA in XNAV and should be updated periodically [6]. With XCOM-based simultaneous communication and ranging, the phase revolution model could be updated in a timely fashion, which would benefit deep space missions far from the Earth. In a word, XCOM could not only provide the accurate observation variable needed for XNAV augmentation, but also a communication service for the transmission of key information. XCOM would be a great supplement for XNAV and provide the possibility of integrated X-ray navigation and communication for future deep space explorations.

Coordinate Systems
To calculate the accurate TOA, the time measurement should be conducted in an inertial coordinate frame. For pulsar observations, the solar system barycenter (SSB) frame is the preferable coordinate. The SSB frame is also called the International Celestial Reference Frame (ICRF), which is an inertial reference system that uses the solar system barycenter as the center of the frame [22].
When in the vicinity of Mars, the Mars-centered coordinate system should be utilized. The fundamental Mars-centered inertial coordinate system is the Mars-centered Earth Mean Equator and Equinox of Epoch [23]. This coordinate is formed by moving the Earth Mean Equator and Equinox of Epoch reference system from the Earth center to the Mars center. The Earth Mean Equator and Equinox of Epoch reference system is defined as follows: (1) the reference plane is set to be the Earth mean equator; (2) the reference direction is set to be the vernal equinox of the Earth; (3) the J2000.0 epoch is selected as the reference epoch. Another commonly used Mars-centered inertial reference system is the Mars-centered Mars Mean Equator and Equinox of Epoch, which is derived from the Mars-centered Earth Mean Equator and Equinox of Epoch. The Mars-centered Mars Mean Equator and Equinox of Epoch is defined as follows: (1) the reference plane is the Mars mean equator of J2000.0; (2) the reference direction is the ascending node of the Mars mean equator relative to the Earth Mean Equator; (3) the reference is set to the J2000.0 epoch. Currently, the Mars-centered Mars Mean Equator and Equinox of Epoch is the most popular reference system in Mars research.

Orbit Dynamics
We define the system state vector as Let ( ) f ⋅ denote the description function of the spacecraft dynamics. Then, the spacecraft state can be expressed as: is the process noise. Let X(t)  be the optimal estimation of the system state and define If the process noise is ignored, ( ( )) f t X can be expressed by the velocity, v , and the acceleration, a : Then, ( ) F t can be expressed as: where I is the unit matrix. The decentralized form of Equation (3) can be written as: where k Φ is the state transition matrix that satisfies: Generally, k Φ can be calculated by the Runge-Kutta (RK) method [24]. Considering the perturbations, the acceleration of a spacecraft orbiting the Mars, a , can be written as: where 0 a is the gravitational acceleration of the Mars, p a is the acceleration caused by various perturbation sources, M μ is the gravitational constant of the Mars, and || || r = r . The perturbation sources of Mars orbit include the non-spherical perturbation, the third-body perturbation, the solar radiation pressure perturbation, and the atmosphere drag perturbation, etc. For the high Mars orbits, the main perturbation sources are the non-spherical perturbation and the third-body perturbation. Other perturbations, such as the atmosphere effect, the solar radiation pressure, etc., are ignored for the high orbit case.
The non-spherical perturbation is caused by the irregularity of the mass distribution and the shape of Mars. The gravitational potential function of Mars can be expressed as [24]: (9) where ( ) nm P ⋅ is the n-order m-degree associated Legendre function, R is the average radius of Mars, r is the distance of the spacecraft relative to the Mars center, λ is the longitude of the spacecraft, ϕ is the reduced latitude. lm C and lm S are the gravitational potential parameters that reflect the mass distribution of Mars, whose values can refer to the Goddard Mars Model: GMM-2B [25]. The acceleration can be expressed as: where ∇ is the gradient computation. Define the zonal harmonic terms, l J as: Other gravitational potential parameters are named as the sectorial harmonic terms ( m l < ) and the tesseral harmonic terms ( m l = ). 2 J perturbation is the main part of the non-spherical perturbation. Thus, the acceleration caused by the non-spherical perturbation can be expressed as: M non The third-body perturbation is caused by the gravitational force from other celestial bodies, mostly from the Sun and the Earth. The third-body acceleration, denoted by third body − a , could be expressed as:  (13) where s μ is the gravitational constant of the Sun, e μ is the gravitational constant of the Earth, s r is the position vector of the Sun in the Mars inertia coordinate, sp r is the position vector of the Sun relative to the spacecraft, e r denotes the position vector of the Earth in the Mars inertia coordinate, and ep r is the position vector of the Earth relative to the spacecraft. Based on the aforementioned analysis, the perturbation acceleration can be calculated as:

Measurement Models
As for the proposed XNAV augmentation method, both the pulsar signal and the range measurement obtained based on XCOM are observed. In this section, the observation equations are established.

Pulsar Timing Observation
Pulsar signals are pulsed period signals with high stability. At the Solar System Barycenter (SSB), a pulsar phase revolution model is established based on the long-duration pulsar observation, which can be utilized to predict the Time of Arrival (TOA) of pulsar pulses at SSB [26]. By observing the TOA of the pulsar signals and comparing the TOA with the time predicted by the phase prediction model [27,28], one can obtain the position of the spacecraft [7]. As the TOA is observed at the spacecraft, it should transferred to Solar System Barycenter (SSB) [7]. During the time transfer, both the geometric effect and other effects, for example, the relativistic effects, should be taken into account. The principle of the time transfer is indicated in Figure 1, where b is the vector pointing from the Sun center to SSB, where ( , ) h t X is a nonlinear function of X , i denotes the i th pulsar that is observed and i v is the observation noise, which is generally modeled as a white Gaussian noise.
As defined previously, r satisfies SSB be the estimation of the system state. Then, the measurement residual can be defined as: In Equation (17), the first order approximation of ( , ) h t X is utilized, i.e.: where ( ) t ε are the high-order terms that are ignored. As the estimation of the state is T T T

X = [r , v ]
   , we have: where SSB MSSB = + r r r   . Then, the measurement residual can be written as: where p H is the pulsar observation matrix defined as: and v is the pulsar observation noise. The decentralized form of Equation (20) can be expressed as: As to the pulsar timing observation, the observation noise is related to the pulsar parameters:  Figure 1. Pulsar timing observation scheme.

Range Observation Based on XCOM
XCOM is a revolutionary concept in deep space exploration, which utilizes X-ray to convey information [29]. As X-rays are high-energy light and have little attenuation in vacuum, they are suitable for deep space exploration. Besides, X-ray detectors can be quite small and energy-saving, which makes them suitable for long distance and long duration missions.
In this paper, the range between the spacecraft and a known location near the Earth is utilized as an observation to improve the performance of XNAV. For the range measurement, an X-ray circularly polarized ranging method based on XCOM has been proposed in our previous research, which will be illustrated in what follows.
Based on the idea of XCOM, we propose a circularly polarized ranging method (XCPolR) [20]. The proposed method utilizes circularly polarized X-ray as the ranging signal. The generation and detection of X-ray polarization signal can refer to the related references [30][31][32][33]. The polarized ranging signal can be expressed as: where r( ) t is the ranging signal, M is the length of the ranging signal, ( ) k Ψ determines the polarization states to be either the left-hand circular polarization state or the right-hand circular polarization state. ( ) k Ψ is defined as:  x y x y x y is the intensity of the light vector components and δ is the phase error between the two orthogonal light vector components. These three parameters can be found in Figure 2. As shown in Figure 2, the two orthogonal light vector components, x E and y E , are defined as: are the phases of the components, and ω is the angular frequency.
In essence, the X-ray range measurement is a two-way measurement, whose procedure can be summarized as: (1) At the emitter, modulate the X-ray signal with the binary ranging code using the circular polarization modulation and send the modulated X-ray signal to the receiver. (2) Through the propagation in space channel, the X-ray ranging signal is received by the receiver and the signal is demodulated to recover the binary ranging sequence.
The estimated range measurement can be defined as: Then, the observation residual can be written as: where || || r = r   , and r H is the range observation matrix defined as: and r v is the ranging jitter. According to the analysis in Reference [20], the ranging jitter mainly comes from the acquisition of the ranging signal. The ranging jitter can be evaluated by the standard variance of the delay measurement error expressed by:

Kalman Filtering Algorithm for State Estimation
As for the proposed method, both the pulsar timing observation and the range measurement are considered in the position determination of a spacecraft. In previous sections, the pulsar timing observation and the range measurement have been analyzed. Here, we would like to establish the observation equation based on the pulsar observation and the range measurement.
Assume two pulsars are observed and one range measurement is conducted. The observation equation can be established as: where ( 1,2) i i σ = is the pulsar timing observation noise related to pulsar parameters and r σ is the ranging jitter determined by Equation (34).

Based on the state transition equation and the observation equation, the Kalman filtering process could be established:
(1) State prediction As indicated by the state transition equation, the one-step state prediction can be expressed as: where k Φ is the state transition matrix, k X  is the optimal estimation of the system state at k , and 1 k + X is the one-step prediction of the system state at +1 k . k w is the process noise whose covariance matrix is defined as The covariance matrix of the prediction error, 1k + P , can be written as: where k P  is the covariance matrix of the prediction error at k .
(2) Filtering gain calculation The filtering gain can be calculated as

Simulation Conditions
The area of the X-ray detector is set to 10,000 cm 2 and the background flux intensity is set to be 0.0050 ph/cm 2 /s. The parameters of the pulsars used in the simulation are listed in Table 1.
The initial error is set to: The covariance matrix of the process noise is set to: The initial covariance matrix of the prediction error, 0 P , is set to: The covariance matrix of the observation noise is calculated based on Equations (24) and (34), respectively.
The accurate orbits are simulated by the Satellite Tool Kit (STK), which is a professional software for simulating the design, test, and operation of space missions. The High-Precision Orbit Propagator (HPOP) is utilized to generate the Mars orbit in this paper. The parameters of the designed orbit are set as follows: (1) the orbit is a circular orbit with a radius of 46,792.48 km; (2) the orbit inclination is 45°. For simplicity, this orbit is named as "Orbit 1".
The position error and velocity error can be evaluated by the following equations:

State Estimation Based on Kalman Filtering
In this section, the Kalman filtering results obtained based on the proposed method are presented and comparisons are made between the proposed method and the traditional pulsar observation method that estimates the system state by observing three pulsars.
The Mars orbit used in the simulation is Orbit 1, whose parameters has been given in Section 3.1. The filtering interval is set to 100 s and the pulsar timing observation interval is set to 800 s. The first two pulsars listed in Table 1 are utilized to achieve pulsar timing observation. The SNR of the ranging signal is set to be 0 dB and the slot duration, rc T , is set to be 1e−7 s. Figure 3 presents the state estimation results of the proposed method. It can be seen from the figures that the position and velocity error can converge to a low level, which demonstrates the feasibility of the proposed method.
The comparison between the proposed method and the traditional method with only pulsar observation are conducted. For the proposed method, the range measurement, together with the timing observation of the first two pulsars listed in Table 1 are utilized for the state estimation. For the traditional method, all three pulsars listed in Table 1 are utilized for the state estimation. Figure 4 presents the comparison of the proposed method with the traditional method. Though the proposed method shows little superiority in velocity estimation compared with the traditional method, the proposed method does show better performance in position estimation under the same condition. The mean position estimation error of the traditional method remains at a level of 200 m, while the proposed method obtains a mean position error of 124 m. Besides, the standard variance of the position estimation error of the traditional method is 102.7 m, which is much larger than that of the proposed method whose standard variance is 65 m. This indicates that the proposed method shows better stability in the estimation of the position.

Impact of Pulsar Observation Interval on Filtering
To analyze the influence of the pulsar observation time on the filtering, a numerical experiment is conducted. The orbit used in the simulation is Orbit 1. The first two pulsars listed in Table 1 are utilized to achieve pulsar timing observation. The SNR of the ranging signal is set to be 0 dB and the slot duration is set to  Figure 5a, with the increasing pulsar observation time intervals, both the mean and standard variance of the position estimation error decrease. It can also be observed from Figure 5b that the mean value and standard variance of the velocity show a similar tendency with the increasing observation time intervals as that of the position estimation. The reason for the phenomenon is that the noise of the pulsar observation is closely related to the observation time. As pulsar signals are weak, the observed signal is accumulated over a certain duration to strength the pulsar signal. The longer the observation time is, the higher the SNR of the pulsar signal is. Thus, long pulsar observation time could mitigate the observation noise, and consequently improve the performance of the position and velocity estimation.

SNR of Ranging Signal on Filtering
As the range measurement is utilized as the observation variable in the proposed method, the accuracy of the range measurement would have an influence on the filtering performance. On the other hand, the range measurement accuracy is highly related to the SNR of the ranging signal. Thus, in this section, we analyze the impact of the ranging signal SNR on the position and velocity estimation.
In the simulation, the initial conditions are as follows: (1) the orbit utilizes Orbit 1; (2) the filtering interval is set to 100 s and the pulsar observation interval is set to 1000 s, respectively; (3) SNR of the ranging signal is set to be −8 dB, 0 dB, 8 dB, respectively. Based on the initial conditions, the filtering of the proposed method is conducted, and the results are plotted in Figure 6. As shown in the figure, the influence of SNR of ranging signal on position and velocity estimation is minor. This is because of the differences in level of noises of the pulsar observation and the ranging measurement. The noise level of the range measurement is much lower than that of the pulsar observation. Thus, when the noise of the range measurement is influenced by the variation of the SNR of the ranging signal, the influence of the change in the range measurement noise is minor because the change is small relative to the high level of the pulsar observation noises.

Filtering Algorithm on Different Orbits
Different Mars orbits are used in this section to analyze the performance of the proposed method. Besides of Orbit 1 introduced in Section 3.1, we introduce another orbit, which is a stationary orbit of Mars. Its parameters are as follows: (1) the orbit is a circular orbit with a radius of 20,429.893717 km; (2) the inclination of the orbit is 0.014°. For simplicity, this orbit is named as "Orbit 2". Other parameters are set as follows: (1) the filtering interval is set to 100 s and the pulsar observation time is set to 1000 s; (2) The SNR of the ranging signal is set to be 0 dB and the slot duration, rc T , is set to be 1e−7 s; (3) the first two pulsars in Table 1 are utilized for pulsar observation. Figure 7 presents the filtering results of the proposed method in different orbits. As shown in Figure 7, the performance of the position and velocity estimation on Orbit 1 and Orbit 2 shows much similarity.
This indicates that the proposed method applies not only to the high Mars orbits, but also to the geostationary Mars orbits.

Conclusions
The lack of accuracy of XNAV has become a huge obstacle for its further development and application. In this paper, a novel XNAV augmentation method is documented, aiming at improving the XNAV accuracy. The proposed method utilizes X-ray ranging as the extra observation, together with traditional pulsar timing observation, to achieve the position and velocity estimation. The performance of the proposed method was verified by numerical experiments. It is demonstrated that the proposed method shows better performance than the traditional pulsar observation method in position and velocity estimation. Numerical experiments also indicated that the performance of the proposed method was influenced by several parameters, such as the pulsar observation interval, the SNR of the ranging signal, etc. Two different orbits were utilized to verify the feasibility of the proposed method on different types of orbits.