1. INTRODUCTION
High-orbit satellites, with an orbit altitude of more than 20 000 km, have important applications in meteorological detection, disaster warning, and data relay (Michael et al., Reference Michael, Davis and Carpenter2002), and many have been launched in recent years. Currently, orbit determination of high-orbit satellites is mainly supported by ground tracking systems. However, the uncertainty and hostility of the space environment demand that satellites be capable of more independent operation (Wen et al., Reference Wen, Zhang and Wu2012). This means that high-orbit satellites should still be able to function effectively if the satellite to ground tracking system connection is lost for long time periods. An autonomous navigation system with high reliability and stability is critical to achieving this.
For high-orbit satellites, Global Navigation Satellite Systems (GNSS) cannot be directly adopted, because of limited visibility and weak signals (Wen et al., Reference Wen, Zhang and Wu2012). Celestial navigation systems (CNS) using stellar angle measurement have been developed for autonomous navigation. They operate by measuring the angle between the line of sight from a given satellite to Earth and the line of sight from a reference star, and can be applied to all outer space. The sensors required by the method are of small size, and the corresponding measurement sampling period is short. However, the navigation performance of the method degrades as the orbit altitude of the satellite increases, and for high-orbit satellites, the best positioning accuracy achievable is of the order of around 1 km (Fang et al., Reference Fang and Ning2006).
X-ray pulsars are rapidly rotating neutron stars that are distant from Earth and generate periodic electromagnetic radiation in the X-ray band (Graven et al., Reference Graven, Collins and Sheikh2008a). The spin periods of X-ray pulsars are highly stable over long periods, and several millisecond pulsars can even match the quality of current atomic clocks for accuracy (Tian et al., Reference Tian, Tang and Yan2012). Adopting pulsars as accurate beacons to fix the positions of Earth-orbiting satellites has been investigated since the 1970s (Downs, Reference Downs1974; Hanson, Reference Hanson1996). In particular, the Advanced Concepts Team (ARIADNA) of the European Space Agency (ESA) studied the feasibility of spacecraft navigation relying on pulsar timing information in 2004, and the Defense Advanced Research Projects Agency (DARPA) of the United States initiated a program called X-ray Source-based Navigation for Autonomous Position Determination in 2005 (Sala et al., Reference Sala, Urruela and Villares2004; Graven et al., Reference Graven, Collins and Sheikh2008b).
The signals from X-ray pulsars can be detected over the whole of outer space. It has been found that the navigation performance of an X-ray pulsar-based navigation system (XNAV) would not be significantly affected by the increase of orbit altitude, and thus XNAV is a good candidate to provide reliable and stable position information of high-orbit satellites. However, the flux of pulsars is usually very low, and the pulsed signal from pulsars is weak and discrete. A pulse Time of Arrival (TOA) at a satellite can only be calculated after a period of observation, which usually lasts for several minutes. Furthermore, the navigation error accumulating during the period of observation cannot be ignored, and it would degrade the performance of XNAV. Moreover, an X-ray detector with a large area, which was usually set to be 1 m2, was assumed to be adopted to enhance the Signal Noise Ratio (SNR) of the signal detected (Xiong et al., Reference Xiong, Wei and Liu2009). However, it is impracticable for most satellites to carry three such X-ray detectors to fix position.
Using current techniques, many systematic errors exist in XNAV. Firstly, the direction of the pulsar cannot be determined with high accuracy (Liu et al., Reference Liu, Ma and Tian2010a). Moreover, for Earth-orbiting satellites, in order to fulfil the time transfer equation, which is the principle of XNAV, the position of Earth predicted by planetary ephemeris such as DE405 should be adopted, and the error existing in the planetary ephemeris would affect the performance of XNAV (Wang et al., Reference Wang, Zheng and Sun2013a). Furthermore, considering that the pulse TOA is recorded by a satellite-borne atomic clock, the error in this clock is also a factor affecting the navigation performance. These factors have been separately analysed previously.
As XNAV and CNS are complementary to each other, the concept of integrated navigation using pulsar and stellar angle measurements has been proposed (Liu et al., Reference Liu, Ma and Tian2010b; Wang et al., Reference Wang, Sun and Zheng2013b). In the work of Liu, XNAV and CNS work separately, and their results are fused by a federated unscented Kalman filter (FUKF). In addition, CNS operates during a period of pulsar observation, and only one X-ray detector with an area of 1 m2 is employed. However, in the case that only a pulsar is observed over the whole navigation process, the observability of XNAV is poor, and it would deliver an unsatisfactory navigation performance (Zheng et al., Reference Zheng, Liu and Qiao2008). Furthermore, the output of the federated filter is not globally optimal, because the two navigation subsystems use the same dynamic model and their results are correlative (Yang, Reference Yang2003; Yang, Reference Yang and Gao2004a). In our prior work, a joint unscented Kalman filter (JUKF) was used to fuse the original measurements. Although a JUKF can reach a globally optimal estimation, the computation burden is high, numerical conditioning difficulties may result and only the impact of clock error is considered.
In this paper, a new autonomous navigation method for high-orbit satellites is proposed that integrates XNAV and CNS measurements. Three X-ray detectors with small areas are adopted and an improved kinematic and static filter is derived to obtain globally optimal navigation results. In addition, the factors affecting the performance of XNAV will be handled by the technique of separate-bias.
The integrated navigation system is described in Section 2, the improved kinematic and static filter is presented in Section 3, simulations are described in Section 4 and followed by the conclusion in Section 5.
2. OVERALL DESIGN OF INTEGRATED NAVIGATION SYSTEM
2.1 Dynamic Model
An Earth-centred inertial coordinate system is selected, and the dynamic model of a high-orbit satellite can be described as
$$\left[ {\matrix{ {{\bi \dot r}} \cr {{\bi \dot v}} \cr}} \right] = \left[ {\matrix{ {\bi v} \cr {\bi a} \cr}} \right] + \left[ {\matrix{ {{\bi w}_r} \cr {{\bi w}_v} \cr}} \right]$$where r is the position vector of satellite, v is the velocity vector of satellite, w=[wrT, wvT]T is the process noise which can be assumed as a zero-mean Gaussian white noise process, and a=aTB+aNS+aT+aH.O.T is the acceleration of satellite about Earth including the following terms (Battin, Reference Battin1999).
aTB=−μ Er/||r||3 is the two-body gravitational acceleration of Earth, μ E is the gravitational coefficient of Earth, and ||•|| denotes the magnitude of vector.
aNS=∂U NSE/∂r is the non-spherical perturbation acceleration of Earth, and U NSE can be expressed as
$$U_{NSE} = - \displaystyle{{\mu _E} \over {\left\| {\bi r} \right\|}}\sum\limits_{n = 2}^\infty {\left( {\displaystyle{{R_e} \over {\left\| {\bi r} \right\|}}} \right)^n J_n P_n \sin \phi} - \displaystyle{{\mu _E} \over {\left\| {\bi r} \right\|}}\sum\limits_{n = 2}^\infty {\sum\limits_{m = 1}^n {\left( {\displaystyle{{R_e} \over {\left\| {\bi r} \right\|}}} \right)^n J_{n,m} P_{nm} \sin \phi \cos m(\lambda - \lambda _{n,m} )}} $$where R e is the radius of Earth, φ and λ are the latitude and longitude respectively, P n and P nm are the Legendre polynomials, J n is the zonal coefficient, λ n,m is the tesseral harmonic term, and J n,m is the tesseral harmonic coefficient.
${\bi a}_T = \sum\limits_{i = 1}^2 {\mu _i [{\bi r}_i /\left\| {{\bi r}_i} \right\|^3 - ({\bi r} - {\bi r}_i )/\left\| {{\bi r} - {\bi r}_i} \right\|^3 ]} $ is the third-body gravitational perturbation acceleration including the impacts of Sun and Moon, μ i is the gravitational coefficient of the ith celestial body, and ri is the position vector of the ith celestial body with respect to Earth.
aH.O.T represents all higher-order terms that may affect acceleration but are nominally considered negligible compared to the remaining effects.
2.2 XNAV Measurement Model
It is of note that the time coordinate system of the whole paper is assumed as the Barycentric Dynamical Time (TDB). Through processing the pulsed signal from the ith pulsar, the pulse TOA at the satellite, t Si, can be measured, and its corresponding arrival time at the Solar System Barycenter (SSB), t SSBi, can be predicted by the pulse timing model. Considering the geometric and relativistic effects, the time transfer equation can be expressed as (Sheikh et al., Reference Sheikh, Pines and Ray2006):
$$\eqalign{t_{SSB} ^i =\, & t_S ^i + \displaystyle{{{\bi n}^i \cdot {\bi r}_S} \over c} + \displaystyle{1 \over {2cD_0 ^i}} \left[ {({\bi n}^i \cdot {\bi r}_S )^2 - \left\| {{\bi r}_S} \right\|^2 + 2({\bi n}^i \cdot {\bi b})({\bi n}^i \cdot {\bi r}_S ) - 2({\bi b}\cdot {\bi r}_S )} \right] \cr & + \displaystyle{{2\mu _S} \over {c^3}} \ln \left| {\displaystyle{{{\bi n}^i \cdot {\bi r}_S + \left\| {{\bi r}_S} \right\|} \over {{\bi n}^i \cdot {\bi b} + \left\| {\bi b} \right\|}} + 1} \right|} $$where ni is the direction vector of the ith pulsar, rS is the position vector of the satellite relative to the SSB, μ S is the gravitational coefficient of Sun, b is the position of the SSB with respect to Sun, c is the speed of light, and D 0i is the distance between the ith pulsar and Sun.
For an Earth-orbiting satellite, rS can be obtained by
$${\bi r}_S = {\bi r} + {\bi r}_E $$where rE is the position vector of Earth relative to the SSB.
It can be seen from Equation (3) that the systematic errors including pulsar direction error, δni, Earth ephemeris position error, δrE, and clock error, δt Si, are the main factors that would significantly affect the performance of XNAV. Therefore, those systematic errors should be taken into consideration during the derivation of the XNAV measurement model.
Then, Equation (3) can be transformed into
$$\eqalign{\hskip-1.5pt t_{SSB} ^i = & t_S ^i + \displaystyle{{{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )} \over c} + \displaystyle{1 \over {2cD_0 ^i}} \{ [{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )]^2 - \left\| {{\bi r} + {\bi \tilde r}_E} \right\|^2 + 2({\bi \tilde n}^i \cdot {\bi b})[{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )] \cr & - 2[{\bi b}\cdot ({\bi r} + {\bi \tilde r}_E )]\} + \displaystyle{{2\mu _S} \over {c^3}} \ln \left| {\displaystyle{{{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \left\| {{\bi r} + {\bi \tilde r}_E} \right\|} \over {{\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|}} + 1} \right| \cr & + \displaystyle{1 \over c}B_1^i (\delta {\bi r}_E ) + \displaystyle{1 \over c}B_2^i (\delta {\bi n}^i ) + \delta t_S^i} $$where 
${\bi \tilde n}^i $ is the measured direction vector of ith pulsar, 
${\bi \tilde r}_E $ is the predicted position of Earth, which can be obtained by planetary ephemeris such as DE405, and B 1i(δrE) and B 2i(δni) are the system biases caused by Earth ephemeris error and the error in ith pulsar direction vector.
Ignoring the second and higher order terms, B 1i(δrE) and B 2i(δni) can be modeled by the following equations.
$$\eqalign{B_1 ^i (\delta {\bi r}_E ) = & {\bi \tilde n}^i \cdot \delta {\bi r}_E + \displaystyle{1 \over {D_0 ^i}} \{ - ({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi r}_E + [{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )]({\bi \tilde n}^i \cdot \delta {\bi r}_E ) - {\bi b}\cdot \delta {\bi r}_E \cr & + ({\bi \tilde n}^i \cdot {\bi b})({\bi \tilde n}^i \cdot \delta {\bi r}_E )\} + \displaystyle{{2\mu _S} \over {c^2}} \displaystyle{{{\bi \tilde n}^i \cdot \delta {\bi r}_E + [({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi r}_E ]/\left\| {{\bi r} + {\bi \tilde r}_E} \right\|} \over {{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \left\| {{\bi r} + {\bi \tilde r}_E} \right\| + {\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|}}} $$
$$\eqalign{B_2 ^i (\delta {\bi n}^i ) = & ({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi \tilde n}^i + \displaystyle{1 \over {D_0 ^i}} \{ [{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )][({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi \tilde n}^i ] \cr & + ({\bi b}\cdot \delta {\bi \tilde n}^i )[{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )] + ({\bi \tilde n}^i \cdot {\bi b})[({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi \tilde n}^i ]\} \cr & + \displaystyle{{2\mu _S} \over {c^2}} \displaystyle{{({\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|)[({\bi r} + {\bi \tilde r}_E )\cdot \delta {\bi \tilde n}^i ] - [{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \left\| {{\bi r} + {\bi \tilde r}_E} \right\|]({\bi b}\cdot \delta {\bi \tilde n}^i )} \over {[{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \left\| {{\bi r}_S} \right\| + {\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|][{\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|]}}} $$In Equation (7), δni is defined as
$$\delta {\bi n}^i = \left[ \hskip-2pt{\matrix{ { - \sin \theta \cos \varphi} \cr { - \sin \theta \sin \varphi} \cr {\cos \theta} \cr} \hskip-2pt} \right]\delta \theta + \left[ \hskip-2pt {\matrix{ { - \cos \theta \sin \varphi} \cr {\cos \theta \cos \varphi} \cr 0 \cr}} \hskip-2pt \right]\delta \varphi $$where θ and φ are the declination angle and right ascension angle of the ith pulsar respectively, and δθ and δφ are the errors in the declination angle and right ascension angle respectively.
Furthermore, for Earth-orbiting satellite, r≪r E, and then the term, 
${\bi r} + {\bi \tilde r}_E $, in Equations (6)–(7) can be approximated by 
${\bi \tilde r}_E $.
However, based on Equation (3), several systematic errors, which might affect the performance of XNAV and should be analysed further, are listed in Table 1.
Table 1. XNAV systematic errors contributing to navigation error.
In earlier works, the impacts of B 1(δrE) and B 2(δni) have been analysed carefully, and they were both modelled as slowly time-varying processes (Wang et al., Reference Wang, Zheng and Sun2013a; Liu et al., Reference Liu, Ma and Tian2010a). Moreover, for high-orbit satellites, the pulsar observation period is usually several minutes, and the drift of clock error can be ignored (Chen et al., Reference Chen, Huang and Lu2011). Therefore, the impacts of B 1i(δrE), B 2i(δni), and δt Si can be handled together.
Assume
$$B^i = B_1^i (\delta {\bi r}_E ) + B_2^i (\delta {\bi n}^i ) + c\delta t_S^i $$The state model of Bi can be described by
$$\dot B^i = 0 + w_B $$where w B is the process noise.
Assume the measurement Zp=[c(t SSB1−t S1), …, c(t SSBi−t Si), …, c(t SSBN−t SN)]T, the measurement noise Vp = [v 1, …, vi, …, vN]T, and the system bias Bp = [B 1, …, Bi, …, BN]T, where 1≤i≤N and N is the number of observed pulsars. The measurement model can be presented as:
$${\bi Z}_p = {\bi h}_P ({\bi x}) + {\bi B}_p + {\bi V}_p $$where hp(x) = [hp1(x), …, hpi(x), …, hpN(x)]T, and hpi(x) is the measurement equation observing ith pulsar. The expression of hpi(x) is:
$$\eqalign{{\bi h}_P ^i ({\bi x}) = \, & {\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \displaystyle{1 \over {2D_0 ^i}} \{ [{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )]^2 - \left\| {{\bi r} + {\bi \tilde r}_E} \right\|^2 + 2({\bi \tilde n}^i \cdot {\bi b})[{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E )] \cr & - 2[{\bi b}\cdot ({\bi r} + {\bi \tilde r}_E )]\} + \displaystyle{{2\mu _S} \over {c^2}} \ln \left| {\displaystyle{{{\bi \tilde n}^i \cdot ({\bi r} + {\bi \tilde r}_E ) + \left\| {{\bi r} + {\bi \tilde r}_E} \right\|} \over {{\bi \tilde n}^i \cdot {\bi b} + \left\| {\bi b} \right\|}} + 1} \right|} $$In Equation (11), Vp can be modeled as a zero-mean Gauss white noise whose standard deviation can be given by (Sheikh, Reference Sheikh2005):
$$\sigma _{TOA} = \displaystyle{{W\sqrt {[B_X + F_X (1 - p_f )]d + F_X p_f}} \over {2F_X p_f \sqrt {At_m}}} $$where W is the width of the pulse, B X is the X-ray background radiation flux, F X is the X-ray flux of the pulsar, p f is the pulsed fraction of the pulsar, d is the ratio of the pulse width to the pulse period P, A is the area of the X-ray detector, and t m is the period of XNAV measurement.
2.3 CNS Measurement Model
Figure 1 shows the stellar angle measurement, α. In Figure 1, s is the direction vector of the reference star with respect to the satellite. Then, the CNS measurement is (Qi et al, Reference Qi and Yang2006):
$$Z_{st} = {\bi h}_{st} ({\bi x}) + v_\alpha = \arccos \left( { - \displaystyle{{{\bi r}\cdot {\bi s}} \over {\left\| {\bi r} \right\|}}} \right) + v_\alpha $$where v α is a zero-mean Gaussian measurement noise with covariance Rst, which is determined by the accuracies of measurement sensors.
Figure 1. Stellar angle measurement.
3. IMPROVED KINEMATIC AND STATIC FILTER
3.1 Review of Kinematic and Static Filter
To avoid the problem that the output of a federated filter is suboptimal, Yang (Reference Yang2003) proposed an information fusion filter named “kinematic and static filter” and proved that the output of the proposed filter is globally optimal. The kinematic and static filter is composed of a kinematic filter and a static filter. The kinematic filter works based on the data from the dynamic model and the measurement from one measurement sensor, and the static filter works based on the output of the kinematic filter and the measurements from the other sensors. Therefore, the data from the dynamic model is just used once and the output of the static filter is globally optimal. For the case that two sensors are used, Figure 2 shows the structure of the kinematic and static filter.
Figure 2. Structure of kinematic and static filter.
As is shown in Figure 2, the predicted satellite state at time t k, xk−, and the corresponding error covariance matrix, Pk−, are provided by the dynamic model. The measurement from sensor 1 is y1,k, whose error covariance matrix is R1,k. The solution of the kinematic filter can be derived as:
$${\bi x}_{1,k} ^ + = [({\bi P}_k ^ - )^{ - 1} + {\bi H}_{1,k} ^T ({\bi R}_{1,k} )^{ - 1} {\bi H}_{1,k} ]^{ - 1} [({\bi P}_k ^ - )^{ - 1} {\bi x}_k ^ - + {\bi H}_{1,k} ^T ({\bi R}_{1,k} )^{ - 1} {\bi y}_{1,k} ]$$
$${\bi P}_{1,k} ^ + = ({\bi I} - {\bi K}_{1,k} {\bi H}_{1,k} ){\bi P}_k ^ - $$where
$${\bi K}_{1,k} = {\bi P}_k ^ - {\bi H}_{1,k} ^T ({\bi H}_{1,k} {\bi P}_k ^ - {\bi H}_{1,k} ^T + {\bi R}_{1,k} )^{ - 1} $$In Equations (15)–(17), H1,kT is the corresponding measurement matrix of y1,k.
At the static filter stage, the satellite state predicted by the dynamic model is no longer used. Instead the result of the kinematic filter, x1,k+ and P1,k+, are utilized to fuse with the measurement from sensor 2, y2,k, whose error covariance matrix is R2,k.
The state model in static filter is
$$\matrix{ {{\bi x}_{2,k} ^ - = {\bi x}_{1,k} ^ +} & {{\bi P}_{2,k} ^ - = {\bi P}_{1,k} ^ +} \cr} $$Based on Equation (18), the result of the static filter stage can be presented as:
$${\bi x}_{2,k} ^ + = [({\bi P}_{2,k} ^ - )^{ - 1} + {\bi H}_{2,k} ^T ({\bi R}_{2,k} )^{ - 1} {\bi H}_{2,k} ]^{ - 1} [({\bi P}_{2,k} ^ - )^{ - 1} {\bi x}_{2,k} ^ - + {\bi H}_{2,k} ^T ({\bi R}_{2,k} )^{ - 1} {\bi y}_{2,k} ]$$
$${\bi P}_{2,k} ^ + = ({\bi I} - {\bi K}_{2,k} {\bi H}_{2,k} ){\bi P}_{2,k} ^ - $$where
$${\bi K}_{2,k} = {\bi P}_{2,k} ^ - {\bi H}_{2,k} ^T ({\bi H}_{2,k} {\bi P}_{2,k} ^ - {\bi H}_{2,k} ^T + {\bi R}_{2,k} )^{ - 1} $$Considering only two sensors are adopted, x2,k+ and P2,k+ are the final estimated state and its corresponding error covariance matrix at time t k.
3.2 Improvement on Kinematic and Static Filter
As is shown in section 3.1, the kinematic and static filter can provide a more flexible structure to fuse the data from the dynamic model and that from different sensors. However, it requires that the sampling period of sensors should be the same; in practice, the period of pulsar observation is usually much longer than that of CNS measurement. Furthermore, the kinematic and static filter is designed for a linear or linearized system; the dynamic model, the XNAV measurement model and the CNS measurement model are all nonlinear equations. Moreover, to enhance the performance of an integrated navigation system, the impact of systematic error in XNAV must be reduced. The common method of reducing systematic error is the augmented-state method, but it increases the computation burden. Therefore, it is necessary to improve the kinematic and static filter.
3.2.1 Improvement with Different Sampling Periods
Suggested by the structure of the kinematic and static filter, the process velocity of the static filter can be slower than that of the kinematic filter. Considering that the sampling period of CNS measurement is much shorter than that of XNAV measurement, XNAV measurement should be handled by the static filter. Otherwise, the CNS measurements obtained during one pulsar observation period are not used in navigation but just wasted. Then, we improve the kinematic and static filter to be an integrated filter with two different sampling periods. Figure 3 shows the structure of the integrated filter.
Figure 3. Structure of integrated filter.
3.2.2 Improvement with Unscented Transformation
Unscented transformation (UT) is a method of approximating the way that the mean and covariance of a random variable changes when the random variable undergoes a nonlinear transformation (Julier et al., Reference Julier, Uhlmann and Hugh2000). The unscented Kalman filter (UKF), which utilizes a set of sigma points produced by an unscented transformation to capture the mean and covariance of the state, would not cause linearization error, and can be used to improve the kinematic and static filter.
Thus the kinematic and static filters can both be improved in the form of a UKF. Assume that the system model of integrated navigation system is:
$${\bi x}_k = {\bi f}({\bi x}_{k - 1} ) + {\bi w}_k $$
$$\left[ \hskip-2pt{\matrix{ {{\bi Z}_{st,k}} \cr {{\bi Z}_{p,k}} \cr}} \hskip-2pt \right] = \left[ \hskip-2pt {\matrix{ {{\bi h}_{st} ({\bi x}_k )} \cr {{\bi h}_p ({\bi x}_k ) + {\bi B}_p} \cr}} \hskip-2pt \right] + \left[ \hskip-2pt {\matrix{ {{\bi v}_{st,k}} \cr {{\bi v}_{p,k}} \cr}} \hskip-2pt \right]$$where
$$\matrix{ {E({\bi w}_k ) = 0} & {E({\bi w}_k {\bi w}_k ^T ) = {\bi Q}_k} \cr} $$
$$\matrix{ {E({\bi v}_{st,k} ) = 0} & {E({\bi v}_{st,k} {\bi v}_{st,k} ^T ) = {\bi R}_{st,k}} \cr} $$
$$\matrix{ {E({\bi v}_{\,p,k} ) = 0} & {E({\bi v}_{\,p,k} {\bi v}_{\,p,k} ^T ) = {\bi R}_{\,p,k}} \cr} $$In Equation (22), xk is the state of the satellite at time t k. In Equation (23), Zst,k and Zp,k are the measurements of CNS and XNAV at time t k, hst(xk) and hp(xk) are the measurement equations of CNS and XNAV at time t k, and vst,k, vp,k are the corresponding measurement noises.
Assume that the estimated satellite state and its corresponding error covariance matrix at time t k−1 are:
$$\matrix{ {{\bi \hat x}_{k - 1} = E({\bi x}_{k - 1} )} & {{\bi P}_{k - 1} = E[({\bi x}_{k - 1} - {\bi \hat x}_{k - 1} )({\bi x}_{k - 1} - {\bi \hat x}_{k - 1} )^T ]} \cr} $$Then, the improved filter is given by the following equations.
Kinematic filter:
Step 1. Structure of sigma points and weights
The set of sigma points, {χi,k−1| i=0,…,2n,k⩾1}, is
$$\left\{ \matrix{{\bi \chi} _{0,k - 1} = {\bi \hat x}_{k - 1} \cr \matrix{ {{\bi \chi} _{i,k - 1} = {\bi \hat x}_{k - 1} + \sqrt {n + \xi} \cdot \left( {\sqrt {{\bi P}_{k - 1}}} \right)_i} & {i = 1,2, \cdots, n} \cr} \cr \matrix{ {{\bi \chi} _{i + n,k - 1} = {\bi \hat x}_{k - 1} - \sqrt {n + \xi} \cdot \left( {\sqrt {{\bi P}_{k - 1}}} \right)_i} & {i = n + 1,n + 2, \cdots, 2n} \cr} } \right.$$where n is the number of components contained in state vector, ξ=α 2(n+κ)-n, in which α is used to control the distribution of sigma points and its value is between 0 and 1, besides κ equals 3-n, and 
$\sqrt {{\bi P}_{k - 1}} $ is the Cholesky factor of Pk−1.
The weights of mean values and covariance values are
$$\left\{ \matrix{\omega _0^m = \displaystyle{\xi \over {n + \xi}} \cr \omega _0^c = \displaystyle{\xi \over {n + \xi}} + \left( {1 - \alpha ^2 + \beta} \right) \cr \omega _i^m = \omega _i^c = \displaystyle{1 \over {2\left( {n + \xi} \right)}} } \right.$$where β is the parameter related with the prior distribution of the state and it is usually set to be 2 in the case of Gaussian distribution.
Step 2. Time update
$$\matrix{ {{\bi \chi} _{i,k|k - 1} = {\bi f}({\bi \chi} _{i,k - 1} )} & {{\bi x}_k^ - = \sum\limits_{i = 0}^{2n} {\omega _i^m {\bi \chi} _{i,k|k - 1}}} \cr} $$
$${\bi P}_k^ - = \sum\limits_{i = 0}^{2n} {\omega _i^c [{\bi \chi} _{i,k|k - 1} - {\bi x}_k^ - ][{\bi \chi} _{i,k|k - 1} - {\bi x}_k^ - ]^T + {\bi Q}_k} $$where xk− and Pk− are the predicted satellite state and its corresponding error covariance matrix at time t k.
Step 3. Measurement update
$$\matrix{ {{\bi Z}_{i,k|k - 1} ^{st} = {\bi h}_{st} ({\bi \chi} _{i,k|k - 1} )} & {{\bi Z}_{st,k}^ - = \sum\limits_{i = 0}^{2n} {\omega _i^m {\bi Z}_{i,k|k - 1} ^{st}}} \cr} $$
$${\bi P}_{st,{\bi \tilde Z}_{st,k} {\bi \tilde Z}_{st,k}} = \sum\limits_{i = 0}^{2n} {\omega _i^c [{\bi Z}_{i,k|k - 1} ^{st} - {\bi Z}_{st,k}^ - ][{\bi Z}_{i,k|k - 1} ^{st} - {\bi Z}_{st,k}^ - ]^T + {\bi R}_{st,k}} $$
$${\bi P}_{st,{\bi \tilde x}_k {\bi \tilde Z}_{st,k}} = \sum\limits_{i = 0}^{2n} {\omega _i^c [{\bi \chi} _{i,k|k - 1} - {\bi x}_k^ - ][{\bi Z}_{i,k|k - 1} ^{st} - {\bi Z}_{st,k}^ - ]^T} $$
$$\matrix{ {{\bi x}_{st,k} ^ + = {\bi x}_k^ - + {\bi K}_{st,k} ({\bi Z}_{st,k} - {\bi \hat Z}_{st,k}^ - )} & {{\bi P}_{st,k} ^ + = {\bi P}_k^ - - {\bi K}_{st,k} {\bi P}_{st,\tilde Z_{st,k} \tilde Z_{st,k}} {\bi K}_{st,k}^T} \cr} $$
$${\bi K}_{st,k} = {\bi P}_{st,{\bi \tilde x}_k {\bi \tilde Z}_{st,k}} {\bi P}_{st,{\bi \tilde Z}_{st,k} {\bi \tilde Z}_{st,k}} ^{ - 1} $$where xst,k+ and Pst,k+ are the results of the kinematic filter at time t k.
Static filter:
In the stage of static filter, Bp in Equation (23) is ignored and it will be handled in section 3.2.3.
Based on Equation (18), we have
$$\matrix{ {{\bi x}_{\,p,k} ^ - = {\bi x}_{st,k} ^ +} & {{\bi P}_{\,p,k} ^ - = {\bi P}_{st,k} ^ +} \cr} $$Another UT is applied by
$$\left\{ \matrix{{\bi \varepsilon} _{0,k - 1} = {\bi x}_{\,p,k} ^ - \cr \matrix{ {{\bi \varepsilon} _{i,k - 1} = {\bi x}_{\,p,k} ^ - + \sqrt {n + \xi} \cdot \left( {\sqrt {{\bi P}_{\,p,k} ^ -}} \right)_i} & {i = 1,2, \cdots, n} \cr} \cr \matrix{ {{\bi \varepsilon} _{i + n,k - 1} = {\bi x}_{\,p,k} ^ - - \sqrt {n + \xi} \cdot \left( {\sqrt {{\bi P}_{\,p,k} ^ -}} \right)_i} & {i = n + 1,n + 2, \cdots, 2n} \cr} } \right.$$where 
$\sqrt {{\bi P}_{\,p,k} ^ -} $ is the Cholesky factor of Pp,k−1−.
The static filter only contains the measurement update step, and it can be shown as
$$\matrix{ {{\bi Z}_{i,k|k - 1} ^p = {\bi h}_p ({\bi \varepsilon} _{i,k|k - 1} )} & {{\bi Z}_{\,p,k}^ - = \sum\limits_{i = 0}^{2n} {\omega _i^m {\bi Z}_{i,k|k - 1} ^p}} \cr} $$
$${\bi P}_{p,{\bi \tilde Z}_{p,k} {\bi \tilde Z}_{p,k}} = \sum\limits_{i = 0}^{2n} {\omega _i^c [{\bi Z}_{i,k|k - 1} ^p - {\bi Z}_{p,k}^ - ][{\bi Z}_{i,k|k - 1} ^p - {\bi Z}_{p,k}^ - ]^T + {\bi R}_{p.k}} $$
$${\bi P}_{\,p,{\bi \tilde x}_k {\bi \tilde Z}_{\,p,k}} = \sum\limits_{i = 0}^{2n} {\omega _i^c [{\bi \varepsilon} _{i,k|k - 1} - {\bi x}_{\,p,k} ^ - ][{\bi Z}_{i,k|k - 1} ^p - {\bi Z}_{\,p,k}^ - ]^T} $$
$$\matrix{ {{\bi x}_{p,k} ^ + = {\bi x}_{p,k}^ - + {\bi K}_{p,k} ({\bi Z}_{p,k} - {\bi Z}_{p,k}^ - )} & {{\bi P}_{p,k} ^ + = {\bi P}_{p,k}^ - - {\bi K}_{p,k} {\bi P}_{p,\tilde Z_{p,k} \tilde Z_{p,k}} {\bi K}_{p,k}^T} \cr} $$
$${\bi K}_{\,p,k} = {\bi P}_{\,p,{\bi \tilde x}_k {\bi \tilde Z}_{\,p,k}} {\bi P}_{\,p,{\bi \tilde Z}_{\,p,k} {\bi \tilde Z}_{\,p,k}} ^{ - 1} $$where xp,k+ and Pp,k+ are the outputs of the static filter at time t k.
Although the JUKF and the kinematic and static filter with UT can both achieve a globally optimal estimation, the computation burdens are different. Assume that the dimensions of satellite state, CNS measurements, and XNAV measurements are n, m st, and m p, respectively. We select the times of multiplication processed in a filter as an index to scale the computation cost by a filter. The indices of JUKF and the kinematic and static filter with UT are listed in Table 2.
Table 2. Computation burden indices of JUKF and the kinematic and static filter with UT.
To obtain a desired navigation performance, m st+m p should be greater than 3. In this case, the index of the kinematic and static filter with UT is less than that of the JUKF. Furthermore, compared with the JUKF, the kinematic and static filter with UT handles smaller matrices. Therefore, the kinematic and static filter with UT can reduce the computation burden and enhance the numerical stability.
3.2.3. Improvement with Separate-bias Estimation
The system bias Bp is handled in this subsection.
A common way to reduce the impact of systematic error is augmenting it into the state of the navigation system. However, the technique of augmented-state would cause the filter implementation to require computation with large matrices, which increases the likelihood of numerical conditioning difficulties (Friedland, Reference Friedland1969). To ensure the numerical stability of the implementation of the kinematic and static filter, we adopt the technique of separate-bias proposed by Friedland to reduce the impacts of systematic errors.
The essence of the technique of separate-bias is to decouple the augmented filter into two parallel filters. The first filter, called the bias-free filter, works based on the assumption that the system bias does not exist. The second filter, called the bias filter, works to estimate the bias. Finally, the outputs of the bias-free filter and bias filter can be used to reconstruct the original system state.
In the part of the static filter provided in section 3.2.2, the impact of Bp is ignored and the static filter can be used as the bias-free filter. The bias filter and the final estimation result are given as follows.
Assume that the estimated system bias and its corresponding error covariance matrix at time t k−1 are
$$\matrix{ {{\bi \hat B}_{\,p,k - 1} = E({\bi B}_{\,p,k - 1} )} & {{\bi P}_{k - 1} = E[({\bi B}_{\,p,k - 1} - {\bi \hat B}_{\,p,k - 1} )({\bi B}_{\,p,k - 1} - {\bi \hat B}_{\,p,k - 1} )^T ]} \cr} $$The part of system bias estimation is shown by the following equations.
$$\matrix{ {{\bi B}_{\,p,k} ^ - = {\bi \hat B}_{\,p,k - 1}} & {{\bi P}_{B,k} ^ - = {\bi \hat P}_{B,k - 1} + {\bi Q}_{B,k - 1}} \cr} $$
$$\matrix{ {{\bi P}_{B,{\bi \tilde x}_k {\bi \tilde Z}_{p,k}} = {\bi P}_{B,k} ^ - {\bi S}^T} & {{\bi P}_{B,{\bi \tilde Z}_{p,k} {\bi \tilde Z}_{p,k}} = {\bi SP}_{B,k} ^ - {\bi S}^T + {\bi P}_{p,{\bi \tilde Z}_{p,k} {\bi \tilde Z}_{p,k}}} \cr} $$
$$\matrix{ {{\bi \hat B}_k = {\bi B}_k ^ - + {\bi K}_{B,k} ({\bi Z}_{\,p,k} - {\bi \hat Z}_{\,p,k}^ - - {\bi SB}_k ^ - )} & {{\bi \hat P}_{B,k} = {\bi P}_{B,k}^ - - {\bi K}_{B,k} {\bi P}_{B,\tilde Z_{\,p,k} \tilde Z_{\,p,k}} {\bi K}_{B,k}^T} \cr} $$
$${\bi K}_{B,k} = {\bi P}_{B,{\bi \tilde x}_k {\bi \tilde Z}_{p,k}} {\bi P}_{B,{\bi \tilde Z}_{p,k} {\bi \tilde Z}_{p,k}} ^{ - 1} $$And then, the final estimations at time t k are:
$$\matrix{ {{\bi \hat x}_k = {\bi x}_{\,p,k} ^ + + {\bi V}_k {\bi \hat B}_{\,p,k}} & {{\bi \hat P}_k = {\bi P}_{\,p,k} ^ + + {\bi V}_k {\bi \hat P}_{B,k} {\bi V}_k ^T} \cr} $$
$$\matrix{ {{\bi S}_k = {\bi H}_k {\bi U}_k + {\bi I}_{m_p \times m_p}} & {{\bi V}_k = {\bi U}_k - {\bi K}_{\,p,k} {\bi S}_k} \cr} $$
$$\matrix{ {{\bi U}_k = {\bi V}_{k - 1}} & {{\bi H}_k = \left. {\displaystyle{{\partial {\bi h}_p ({\bi x})} \over {\partial {\bi x}}}} \right|_{{\bi x} = {\bi x}_{p,k} ^ -}} \cr} $$where I is a unit matrix.
Considering that the dimension of system bias is the same as that of XNAV measurement, the indices of the augmented-state method and the separate-bias method are shown in Table 3.
Table 3. Computation burden indices of augmented-state method and separate-bias method.
As shown in Table 3, based on the kinematic and static filter with UT, the separate-bias method costs less computation burden and handles smaller matrices compared to the augmented-state method.
4. SIMULATION AND RESULTS
To verify the performance of the proposed integrated navigation system, some simulations are described in this section.
Three existing satellite orbits of NTS 2, DOGE 1, and ATS 1 are investigated, and the initial orbital elements of the true orbits are shown in Table 4. The navigation pulsars adopted are selected from the X-ray pulsar database provided by Microcosm Incorporated, and the parameters are listed in Table 5. Fifty stars, distributed on the whole celestial sphere (visual magnitude⩽2m), are selected as reference stars. The data of star and orbit are generated by the Satellite Tool Kit (STK).
Table 4. Initial Orbital Elements.
Table 5. Parameters of pulsars selected by Microcosm Incorporated.
* 1kpc=3×1019m.
Assume that the navigation process lasts for 7 days, and the position errors of pulsars are 0·001 arcsec. The XNAV measurement noise standard deviation is determined by Equation (13), with the specified parameters A=0·3 m2, t m=0·5 h, and B X=0·005 ph/cm2/s. The accuracy of star sensor is 3″(1σ), and the accuracy of the optical camera is 0·05°(1σ).
For the three satellites listed in Table 4, the percentages of the time interval when the pulsars are visible for the whole navigation period are provided in Table 6. In the visibility analysis, the impacts of Sun, Moon, and Earth are considered.
Table 6. Percentages of the pulsar visible time interval to the whole navigation period [%].
It can be seen from Table 6 that most of the pulsars are visible during the whole navigation process and the visibility grows as the orbit altitude increases. The result is consistent with the analysis provided by Mao (Reference Mao2009). Among the pulsars that are visible during the whole navigation process, we selected navigation pulsars based on the following criteria: (1) navigation pulsars should not be binary pulsars whose measurement models are not considered in the paper; (2) the distributions of the navigation pulsars should not be close; (3) the flux of navigation pulsar should be as large as possible. Thus PSR B0531+21, PSR B1821-24, and PSR B1509-58 are selected as the navigation pulsars adopted in the following simulations.
Assume that the navigation errors existing in the initial state of the satellite are [1 km, 1 km, 1 km] and [1 m/s, 1 m/s, 1 m/s], and the clock error is 1 μs. Furthermore, DE405 ephemeris is used to predict the Earth position used in the navigation process, and the error contained in DE405 is approximated by the difference between DE405 and DE421, which was created after DE405. The covariance of the state process noise Q=diag(q 12, q 12, q 12, q 22, q 22, q 22), where q 1 is 2×10−5 m, q 2 is 6×10−4 m/s (Xiong et al., Reference Xiong, Wei and Liu2009). The covariance of the bias process noise Q B=diag(q 32, q 32, q 32), where q 3 is chosen as 0·01 m.
Figure 4 shows the performance comparison of the integrated navigation system, XNAV, and CNS. As is illustrated in Figure 4, the three navigation methods can all converge for the three satellites, but the integrated navigation system performs best. Table 7 gives the position accuracy comparison among the three navigation methods over 300 Monte Carlo trials. For the three satellites, we can see that the positioning RMS (Root Mean Square) of the integrated navigation systems are all less than 100 m but the RMS of the other two methods are greater than 200 m. Thus the integrated navigation system behaves better than XNAV and CNS. The reason for this is that the integrated navigation system works by fusing the XNAV measurement and CNS measurement. Moreover, it can be learnt from Figure 4 and Table 7 that the performance of the integrated navigation system is little affected by the increment of the orbit altitude. Thus the proposed integrated navigation system is suitable for high-orbit satellite autonomous navigation.
Figure 4. Positioning comparison of three navigation methods.
Table 7. Positioning RMSs of three navigation methods [m].
From the above results, we can also see that the NTS 2 has the lowest position accuracy among three satellites. To demonstrate the effectiveness of the proposed integrated navigation system, only the result of NTS 2 is considered in the remainder of this paper.
Figure 5 depicts the estimated position error and its corresponding 3σ positioning outlier for the integrated navigation system and XNAV. We can see that the errors of the estimated positions of the proposed integrated navigation system are smaller than 3σ during the whole navigation process, but some of the errors of estimated positions of XNAV are larger than that boundary. This means that the integrated navigation system can reduce the impacts of the systematic errors in XNAV.
Figure 5. Error of estimated positions and 3σ positioning outlier of integrated navigation system and XNAV.
Only the impact of satellite-borne clock error is considered. Figure 6 shows the navigation performance with different satellite-borne clock errors. We can see that the positioning accuracies of two navigation systems decline as the satellite-borne clock error increases but the integrated navigation system outperforms with the same clock error. This means that the performance of the proposed navigation system is little affected by the increment of systematic error. The impacts of pulsar direction error and Earth ephemeris error can also be analysed in the same way.
Figure 6. Error of estimated positions with different satellite-borne clock error.
The impact of linearization error that might occur in the kinematic and static filter is investigated. The navigation schemes designed are listed in Table 8, and the corresponding positioning accuracies over 300 Monte Carlo trials are provided in Table 9.
Table 8. Navigation schemes adopted to analyze the impact of linearization error.
Table 9. Positioning RMSs with different navigation schemes [m].
As shown in Table 9, linearization error would degrade the navigation performance. Compared with the linearization of measurement models, the linearization of the dynamic model would more easily degrade navigation performance. Therefore, to ensure the performance of the integrated navigation system, the dynamic model should not be linearized, but the XNAV measurement model and CNS measurement model can be linearized to reduce computation burden if the corresponding navigation performances are acceptable.
Finally, without consideration of systematic errors, Figure 7 presents the performance comparison between the integrated navigation system using three X-ray detectors. The area of each detector is 0·3 m2, and the integrated navigation system using one X-ray detector with area of 1 m2, which observes PSR B0531+21 during the whole navigation process. From Figure 7, it can be seen that the integrated navigation system using three small detectors outperforms the integrated navigation system using one detector with area of 1 m2. Although it can be seen from Equation (13) that the measurement noise would grow if the area of the detector decreases, three detectors can receive the signals from three pulsars at the same time, and can provide a better geometric structure compared to using only one detector. The impact of the increment of measurement noise can be offset by the performance improvement of a better geometric structure.
Figure 7. Error of estimated positions with different X-ray detector configurations.
5. CONCLUSION
In this paper, a high-orbit satellite autonomous navigation method by integrating X-ray pulsar-based navigation (XNAV) and celestial navigation system (CNS) is proposed. An improved kinematic and static filter is derived to fuse data. In the filter, the unscented transformation (UT) is used to reduce linearization error, and the separate-bias technique is adopted to reduce the impact of systematic errors contained in XNAV, considering that the sampling periods of the sensors are different. Compared with the XNAV and CNS, the proposed navigation system has an improved performance and can reduce the impact of systematic error effectively. The proposed navigation system is suitable for high-orbit satellite autonomous navigation.
ACKNOWLEDGEMENT
This work is supported by National Natural Science Foundation of China (10973048).
APPENDIX
In this appendix, we demonstrate that the result of the federated filter is suboptimal and that the result of the kinematic and static filter is globally optimal.
A1. THE RESULT OF THE FEDERATED FILTER IS SUBOPTIMAL
For simplification, the federated filter with the measurements from two sensors is analysed, and the result can be easily expanded to the case that the federated filter works based on the measurements form numbers of sensors.
The structure of the federated filter with two sensors is depicted in Figure A1. In Figure A1, y1,k and y2,k are the measurements from sensors 1 and 2 at time t k; x− is the predicted state obtained by the dynamic model; x1,k+ and x2,k+ are the estimations from local filters 1 and 2, and P1,k and P2,k are the corresponding covariance matrices. The outputs of two local filters would be fused in the master filter. And the global estimation is as follows (Yang, Reference Yang, Cui and Gao2004b; Yang, Reference Yang2006).
$${\bi x}_{g,k} = [{\bi P}_{1,k} ^{ - 1} + {\bi P}_{2,k} ^{ - 1} ]^{ - 1} [{\bi P}_{1,k} ^{ - 1} {\bi x}_{1,k} ^ + + {\bi P}_{2,k} ^{ - 1} {\bi x}_{2,k} ^ + ]$$
$${\bi P}_{g,k} = [{\bi P}_{1,k} ^{ - 1} + {\bi P}_{2,k} ^{ - 1} ]^{ - 1} $$Then, the global estimation is fed back to two local filters by the following equations.
$${\bi x}_{i,k} = {\bi x}_{g,k} (i = 1,2)$$
$${\bi P}_{i,k} = {\bi P}_{g,k} \beta _i $$
$$\beta _i = \displaystyle{{\left\| {{\bi P}_{i,k}} \right\|^{ - 1}} \over {\left\| {{\bi P}_{1,k}} \right\|^{ - 1} + \left\| {{\bi P}_{2,k}} \right\|^{ - 1}}} $$The expressions of x1+ and x2+ can be presented as
$${\bi x}_{i,k} ^ + = {\bi K}_{i,k} {\bi y}_{i,k} + ({\bi I} - {\bi K}_{i,k} {\bi H}_{i,k} ){\bi x}_k ^ - (i = 1,2)$$where Ki is the Kalman gain, Hi is the measurement matrix of ith local filter, I is a n×n unit matrix, and n is the dimension of state x.
Figure A1. Structure of the federated filter.
Therefore, the correlation matrix between x1,k+ and x2,k+ is
$${\mathop{\rm cov}} ({\bi x}_{1,k} ^ +, {\bi x}_{2,k} ^ + ) = ({\bi I} - {\bi K}_{1,k} {\bi H}_{1,k} ){\bi P}_k ^ - ({\bi I} - {\bi K}_{2,k} {\bi H}_{2,k} )$$As is shown in Equation (A7), the outputs of local filters correlate. However, the master filter fuses the outputs based on the least square criterion which requires the inputs to be uncorrelated. Therefore, the output of the master filter is not globally optimal.
A2. THE RESULT OF THE KINEMATIC AND STATIC FILTER IS GLOBALLY OPTIMAL
It is learnt that the output of the federated filter is suboptimal because the same dynamic model is used in the local filters. However, in terms of the kinematic and static filter, the dynamic model is just used in the kinematic filter, and the static filter works based on the output of the kinematic filter and the measurement of sensor 2. Therefore, the dynamic model is used once, and the result of the static filter is globally optimal.
