2.1. Inter-satellite Range Measurement Model

The inter-satellite range method cannot resolve the absolute position of constellation satellites. However, this method can resolve the relative position with high accuracy.

The dual-directional pseudoranges between two inter-visible satellites A and B are denoted by ρ^AB and ρ^BA as:

(1)$$\rho ^{AB} (t) = \rho _0^{AB} (t) + I^{AB} (t) + c(\delta t_A (t) - \delta t_B (t)) + \varepsilon (t)$$

(2)$$\rho ^{BA} (t) = \rho _0^{AB} (t) + I^{BA} (t) + c(\delta t_B (t) - \delta t_A (t)) + \varepsilon (t)$$

where:

• ρ ₀^AB(t) refers to the geometric distance between two satellites;

• c is the velocity of light;

• δt _A(t) and δt _B(t) denote the clock errors of satellites A and B, respectively;

• I^AB(t) and I^BA(t) refer to the ionosphere delay at reception time of satellites A and B, respectively;

• ε (t) is the one-way observable error or noise, which can be modelled as biased, first-order Gauss–Markov processes.

By summing up the two given equations, we obtain the pseudorange measurement:

(3)$$Z_\rho (t) = (\rho ^{AB} (t) + \rho ^{BA} (t))/2 = \rho _0^{AB} (t) + (I^{AB} (t) + I^{BA} (t))/2 + \varepsilon (t)$$

2.2. Pulsar-based Measurement Model

Pulsars can provide an external spatial reference to satellite constellations, potentially avoiding non-convergence resulting from the orbit determination on the basis of inter-satellite range measurements. As Figure 1 shows, two satellites in the constellation lock on the same X-ray pulsar source and simultaneously detect the same signal emitted from this source. Thus, the Relative Time of Arrival (RTOA) of the pulse to the two satellites can be acquired, the value of which is the difference between the two satellites' Time Of Arrival (TOA).

Figure 1. Relative TOA between two spacecraft that are observing the same pulsar.

The first order of the simplified RTOA equation is expressed as (Sheikh et al., Reference Sheikh and Pines2006):

(4)$$c\Delta t = {\bf n} \cdot ({\bf r}_{\bf A} - {\bf r}_{\bf B} ) = {\bf n} \cdot {\bf r}_{AB} = r_{AB} \cos \theta $$

Thus:

(5)$$\cos \theta = {\raise0.7ex\hbox{${(c\Delta t)}$} \!\mathord{\left/ {\vphantom {{(c\Delta t)} {r_{AB}}}} \right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${r_{AB}} $}}$$

where:

• c is the speed of light;

• Δt denotes the RTOA, which can be obtained from correlation algorithms (Sheikh et al., Reference Sheikh and Pines2006); r_AB is the relative position vector;

• r _AB is the range between the space vehicles A and B, which can be calculated from Equation (3);

• n is the normalized pulsar position vector, which can be acquired by astronomical observation (Sheikh et al., Reference Sheikh and Pines2006);

• θ is the angle contained by n and r_AB. As shown in Figure 2, it is between the direction of the pulsar source and direction of base line which is the line between two satellites detecting the same pulsar signal.

Figure 2. Angle between the pulsar direction and the baseline in the constellation.

To determine the direction of the baseline r_AB in the inertial space, that is, θ, measurements from two or more pulsars are needed.

2.3. Observability Analysis for the Pulsar-based Measurement Model

The observability of a system indicates whether the system states (here the system states are the Keplerian orbital elements) can be determined by the measurements. If the system is not observable, the system states cannot be determined, even if the noise level is negligible.

The three Keplerian orbital elements of a satellite in constellation, that is, semi-major axis (a), eccentricity (e), and mean anomaly (M), can be determined by inter-satellite range measurements. These are the Keplerian orbital elements related to the size and shape of the orbits, and the position of the spacecraft along the orbits. However, difficulties arise when attempting to estimate the other three Keplerian orbital elements that are related to orientation, namely, the right ascension of ascending node (Ω), the argument of perigee (ω), and the orbital inclination (i). (Keric, Reference Keric2007; Zhang, Reference Zhang2005; Liu, Reference Liu and Liu2001).

In order to illustrate whether the pulsar-based measurement model can determine the three orbital elements (Ω, ω, i), we need to analyse the partial differential equations of the pulsar-based measurement model with respect to the three orbital elements (Ω, ω, i). If the partial differential equations are different from zero, the three orbital elements can be determined.

These equations are analysed in detail below

2.3.1. Right Ascension of Ascending Node Ω

The partial differential equations of the pulsar-based measurement model with respect to the orbital element Ω can be expressed as:

(6)$$\eqalign{\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial \Omega}} = & \displaystyle{{\partial ({\bi n}_x \cdot \Delta {\bi x})} \over {\partial \Omega}} + \displaystyle{{\partial ({\bi n}_y \cdot \Delta {\bi y})} \over {\partial \Omega}} + \displaystyle{{\partial ({\bi n}_z \cdot \Delta {\bi z})} \over {\partial \Omega}} \cr & = {\bi n}_x \displaystyle{{\partial x_i} \over {\partial \Omega}} + {\bi n}_y \displaystyle{{\partial y_i} \over {\partial \Omega}} + {\bi n}_z \displaystyle{{\partial z_i} \over {\partial \Omega}} - \left( {{\bi n}_x \displaystyle{{\partial x_j} \over {\partial \Omega}} + {\bi n}_y \displaystyle{{\partial y_j} \over {\partial \Omega}} + {\bi n}_z \displaystyle{{\partial z_j} \over {\partial \Omega}}} \right)} $$

(7)$$\displaystyle{{\partial {\bi r}} \over {\partial \Omega}} = {\bi J}_z \times {\bi r} = (\matrix{ { - {\rm y}} & {\rm x} & 0 \cr} )\quad ({\rm Liuand}\,{\rm Liu},\,2001).$$

Substituting Equation (7) into Equation (6):

(8)$$\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial \Omega}} = - {\bi n}_x y_i + {\bi n}_y x_i - ( - {\bi n}_x y_j + {\bi n}_y x_j ) = {\bi n}_y (x_i - x_j ) - {\bi n}_x (y_i - y_j )$$

where:

• n is the normalized pulsar position vector, which can be obtained from astronomical observations; n_x,n_y and n_z are the components of the pulsar position vector n in the inertial coordinate system;

• r is the satellite position vector, Δr is the relative position vector between two satellites; Δx, Δy, Δz are the three components of the relative position vector in the inertial coordinate system; x _i, y_i, z_i, x_j, y_j, z_j are the three components of satellites i and j in the inertial coordinate system.

From Equation (8), when ${\textstyle{{(x_i - x_j )} \over {(y_i - y_j )}}} \ne {\textstyle{{{\bi n}_x} \over {{\bi n}_y}}} $, the partial differential equation of the pulsar-based measurement model with respect to the orbital element (Ω) is different from zero, and hence solvable with respect to Ω.

2.3.2. Orbital Inclination i

The partial differential equation of the pulsar-based measurement model with respect to the orbital element i is given by:

(9)$$\eqalign{\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial i}} = & \displaystyle{{\partial ({\bi n}_x \cdot \Delta \vec x)} \over {\partial i}} + \displaystyle{{\partial ({\bi n}_y \cdot \Delta \vec y)} \over {\partial i}} + \displaystyle{{\partial ({\bi n}_z \cdot \Delta \vec z)} \over {\partial i}} \cr & = {\bi n} \cdot \displaystyle{{\partial \vec r_i} \over {\partial i}} - {\bi n} \cdot \displaystyle{{\partial \vec r_j} \over {\partial i}}} $$

(10)$$\displaystyle{{\partial {\bi r}} \over {\partial i}} = \vec J_N \times {\bi r} = \left( \matrix{\cos \Omega \cr \sin \Omega \cr 0 } \right) \times \vec r\quad ({\rm Liuand}\;{\rm Liu},\;{\rm 2}00{\rm 1})$$

Substituting Equation (10) into Equation (9), we obtain:

(11)$$\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial i}} = {\bi n} \cdot ({\bi J}_{Ni} \times {\bi r}_i ) - {\bi n} \cdot ({\bi J}_{Nj} \times {\bi r}_j ) = ({\bi n} \times {\bi r}_i ) \cdot {\bi J}_{Ni} - ({\bi n} \times {\bi r}_j ) \cdot {\bi J}_{Nj} $$

Assuming Ω_i=Ω_j, then J_Ni=J_Nj=J_N, and:

(12)$$\eqalign{\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial i}} = & ({\bi n} \times {\bi r}_i ) \cdot {\bi J}_N - ({\bi n} \times {\bi r}_j ) \cdot {\bi J}_N = {\bi J}_N \cdot (({\bi n} \times {\bi r}_i ) - ({\bi n} \times {\bi r}_j )) = {\bi J}_N \cdot ({\bi n} \times ({\bi r}_i - {\bi r}_j )) \cr & = ({\bi n}_y \Delta {\bi r}_z - {\bi n}_z \Delta {\bi r}_y )\cos \Omega + ({\bi n}_z \Delta {\bi r}_x - {\bi n}_x \Delta {\bi r}_z )\sin \Omega} $$

From Equation (12), when${\textstyle{{({\bi n}_y \Delta {\bi r}_z - {\bi n}_z \Delta {\bi r}_y )} \over {({\bi n}_z \Delta {\bi r}_x - {\bi n}_x \Delta {\bi r}_z )}}} \ne {\textstyle{{\sin \Omega} \over {\cos \Omega}}} $, the partial differential equation of pulsar-based measurement model with respect to the orbital element i is different from zero, and hence solvable with respect to i.

2.3.3. Argument of Perigee ω

The partial differential equation of the pulsar-based measurement model with respect to the orbital element ω is given by:

(13)$$\displaystyle{{\partial ({\bi n} \cdot \Delta {\bi r})} \over {\partial \omega}} = {\bi n} \cdot ({\bi R}_i \times {\bi r}_i ) - {\bi n} \cdot ({\bi R}_j \times {\bi r}_j ) = ({\bi n} \times {\bi r}_i ) \cdot {\bi R}_i - ({\bi n} \times {\bi r}_j ) \cdot {\bi R}_j $$

where:

• R is the orbital normal vector given by: ${\bi R} = {\textstyle{1 \over {\sqrt {\mu p}}}} ({\bi r} \times {\bi \dot r})$, (Liu and Liu, Reference Liu and Liu2001).

When two satellites are in the same orbital plane, R_i=R_j, yielding:

(14)$$\displaystyle{{\partial ({\bi n} \cdot \Delta \vec r\hskip 2pt)} \over {\partial \omega}} = ({\bi n} \times {\bi r}_i ) \cdot {\bi R}_i - ({\bi n} \times {\bi r}_j ) \cdot {\bi R}_j = {\bi r}_i \cdot ({\bi n} \times {\bi R}_i ) - {\bi r}_j \cdot ({\bi n} \times {\bi R}_i ) = \Delta {\bi r}_{ij} \cdot ({\bi n} \times {\bi R}_i )$$

From Equation (14), if Δr_ij·(n×R_i)≠0, the partial differential equation of the pulsar-based measurement model with respect to the orbital element ω is different from zero, and hence solvable with respect to ω.

3.1. Dynamic Models

For a constellation comprising of n satellites, the state vector is:

$${\bi X}(t) = [{\bi X}_1 (t),{\bi X}_2 (t),{\bi X}_3 (t),{\bi X}_4 (t), \cdots {\bi X}_n (t)]$$

The common state model of the spacecraft navigation system is given as:

(15)$${\bi \dot X}_i (t) = \bar f\,({\bi X}_i (t) + {\bi \eta} (t) = [{\bi \dot r},{\bi \dot v}]^T = [v_x, v_y, v_z, a_x, a_y, a_z ]^T + {\bi \eta} (t)$$

where:

• X_i(t) is the state vector, ${\bi X}_i (t) = [{\bi \dot r},{\bi \dot v}]^T = [v_x, v_y, v_z, a_x, a_y, a_z ]^T $;

• r=[r _x, r_y, r_z]^T, and v=[v _x, v_y, v_z]^T are the position and the velocity vectors of the spacecraft with respect to the geocentric inertial frame (J2000);

• η(t) is the noise.

The total acceleration for the satellite is defined as:

(16)$${\bi a} = {\bi \ddot r} = {\bi a}_{\bi 1} {\bi + a}_{\bi 2} {\bi + a}_{\bi 3} {\bi + a}_{\bi H} $$

where:

• a ₁ is the gravitational acceleration of the Earth;

• a ₂ is the non-spherical perturbation acceleration of the Earth's gravitational acceleration;

• a ₃ is the gravitational acceleration of the third body;

• a _H is the tide perturbation acceleration and solar radiation pressure perturbation acceleration, corresponding to higher orders in the total acceleration;

• η(t) is the state process noise, assumed to be white noise with a zero-mean and covariance Q(t) (Sebastian, Reference Sebastian and Glenn2011).

3.2. Measurement Models

Two types of measurement models, inter-satellite range models and pulsar-based measurement models, provide information on the satellites' relative positions in the constellation and the direction of the constellation baseline in the inertial system.

For a constellation comprising of n satellites, the measurement equations for inter-satellite range models can be expressed as:

(17)$${\bi Y}_R (t) = H_R ({\bi X},t) + {\bi \varepsilon} (t)$$

In Equation (17) the range measurement Y_R(t) can be expressed as:

(18)$${\bi Y}_R (t) = \left\{ {\matrix{ {y_{i,1} = \rho _0^{i1} + (I^{i1} + I^{1i} )/2 + \varepsilon _1 (t)} \cr \vdots \cr {y_{i,n} = \rho _0^{in} + (I^{in} + I^{ni} )/2 + \varepsilon _n (t)} \cr}} \right.$$

where:

• ρ ₀ⁱⁿ is the geometric distance between the ith and nth satellites;

• Iⁱⁿ is the ionosphere-induced delay between the ith and nth satellites; and

• ε _n(t) is the random measurement noise.

In Equation (17) the measurement matrix H _R(X, t) is given by:

(19)$$\eqalign{& H_R^{i,j} = \left(\displaystyle{{\partial \rho} \over {\partial x_i}} \;\displaystyle{{\partial \rho} \over {\partial y_i}} \;\displaystyle{{\partial \rho} \over {\partial z_i}} \;0\;0\;0\; - \displaystyle{{\partial \rho} \over {\partial x_i}} \; - \displaystyle{{\partial \rho} \over {\partial y_i}} \; - \displaystyle{{\partial \rho} \over {\partial z_i}} \;0\;0\;0\right) \cr & \quad \quad \;\;(x_i \;y_i \;z_i \;v_{xi} \;v_{yi} \;v_{zi} \;x_j \;y_j \;z_j \;v_{xj} \;v_{xj} \;v_{xj} )^T} $$

The pulsar-based measurement is given as:

(20)$${\bi Y}_P (t) = H_P ({\bi X},t) + {\bi \varepsilon} (t)$$

In Equation (20) Y_P(t) can be determined from Equation (5) as:

(21)$${\bi Y}_P (t) = y_{n + k} = \cos \theta _k = \left( {c\Delta t_k} \right)/r_{ij} $$

From Equation (4), the measurement matrix H _P is computed as:

(22)$$H_P = ({\bi n}_{\bi x} \;{\bi n}_y \;{\bi n}_z \;0\;0\;0\; - {\bi n}_{\bi x} \; - {\bi n}_y \; - {\bi n}_z \;0\;0\;0)(x_i \;y_i \;z_i \;v_{xi} \;v_{yi} \;v_{zi} \;x_j \;y_j \;z_j \;v_{xj} \;v_{xj} \;v_{xj} )^T $$

and ε(t) is the zero-mean noise with covariance R(t).

3.3. UKF Algorithm

The UKF (Julier and Uhlmann, Reference Julier and Uhlmann1997) is an extension of the Kalman filter. The most important feature of the UKF is its capacity to compute the mean and the covariance of nonlinearly transformed variables from selected variables of so-called sigma points without the need for a linear approximation of the nonlinear function, as required by the extended Kalman filter. Thus, the UKF does not yield a linearization error. In this paper, the orbit dynamic model and the measurement models are nonlinear. Thus, the UKF is suitable for the pulsar positioning system. The UKF steps are shown in Figure 3.

Figure 3. UKF filter flowchart.

The first step is to initialize the state covariance matrix P ₀. In the second step, 2 l+1 sample sigma points are calculated based on the square-root of the augmented state covariance matrix, where χ is a matrix of state sigma points, and λ is the sigma point spread parameter. Weight vectors (ω^m, ω^c) for the mean and covariance are computed in the third step. In the following time step, the state prediction ${\bi \hat X}^ - $ and state covariance prediction ${\bi \hat P}^ - $ are computed using state dynamic equations and a weighted average of the transformed sigma points. Next, the output prediction ${\bi \hat Z}^ - $ is calculated using the current state sigma points in the observation equations and weighting vector. In the measurement update step, the Kalman gain matrix K _k is calculated using the output covariance matrix $P_{\tilde Z\tilde Z} $ and the covariance matrix $P_{\tilde X\tilde Z} $ between the state prediction and output estimate. Finally, the state and state covariance estimations are updated by the measurement vector Z and the Kalman gain matrix K _k.

4.1. Initial Condition

The simulation conditions are as follows. Assume that a constellation consists of four satellites, and the orbital altitude of all the satellites is 23 000 km. The initial orbit elements of the satellites are shown in the Table 1, and the true orbit ephemeris are provided by the 4×4 EGM96 model. The initial orbit elements of the four satellites in the constellation are presented in Table 1, and the parameters of the pulsars used are shown in Table 2.

Table 1. Initial orbit elements of the satellites.

Table 2. Parameters of pulsars.

The absolute position estimation errors versus time of satellite A are plotted in Figure 4, using pulsar measurements and inter-satellite range measurements. The position errors in the tangential T, and radial R directions converge to approximately 1·5 m and 0·7 m respectively. The position error in the normal direction N converges to a mean of approximately zero, with slow oscillations.

Figure 4. Position estimation errors when using pulsars and inter-satellite range measurement.

As a comparison to the aforementioned method, only the inter-satellite range measurement is utilized as the input to the UKF for the satellite state estimation, using the same initial condition. The position estimation errors for satellite A are shown in Figure 5. As the inter-satellite range measurement just provides radial direction measurements to the positioning system, only the position error in the radial direction is convergent.

Figure 5. Position estimation errors when using only inter-satellite range measurement.

From Figures 4 and 5, we can conclude that inter-satellite range measurement only restricts errors in the radial direction. The range measurement is unobservable to tangential and normal perturbation force. When the relative measurement based on X-ray pulsars is added to the measurement system, the relative measurement via pulsars is observable to the tangential and normal perturbation force, as discussed in Section 2.3. The abovementioned simulation results demonstrate that the combined measurement models accomplish the autonomous navigation of satellite constellation with high accuracy.

4.2. Effect of Pulsar Direction Error on the Measurement Model

Pulsar direction error is an important error source for pulsar positioning. The effect of the pulsar direction error on the relative pulsar-based measurement is analysed.

(23)$${\bi n} = \left[ {\matrix{ {\cos \delta \cos \alpha} \cr {\cos \delta \sin \alpha} \cr {\sin \delta} \cr}} \right]$$

where:

• α is the pulsar right ascension;

• δ is pulsar declination angle;

• n is the pulsar direction vector.

Assuming that the measured pulsar position with error is $(\tilde \alpha, \tilde \delta )$, and the error is (Δα, Δδ), the pulsar real position is given as:

(24)$$\left\{ \matrix{\alpha = \tilde \alpha + \Delta \alpha \cr \delta = \tilde \delta + \Delta \delta } \right.$$

Substituting Equation (24) into Equation (23), the direction vector of the pulsar is given as:

(25)$${\bi n} = \left[ {\matrix{ {\cos (\tilde \delta + \Delta \delta )\cos (\tilde \alpha + \Delta \alpha )} \cr {\cos (\tilde \delta + \Delta \delta )\sin (\tilde \alpha + \Delta \alpha )} \cr {\sin (\tilde \delta + \Delta \delta )} \cr}} \right]$$

Using a first order Taylor expansion (Jiang Reference Jiang and Li2007) in Equation (25), n can be simplified as:

(26)$${\bi n} = \left[ {\matrix{ {\cos \tilde \delta \cos \tilde \alpha} \cr {\cos \tilde \delta \sin \tilde \alpha} \cr {\sin \tilde \delta} \cr}} \right] + \left[ {\matrix{ { - \cos \tilde \delta \sin \tilde \alpha \cdot \Delta \alpha - \sin \tilde \delta \cos \tilde \alpha \cdot \Delta \delta} \cr {\cos \tilde \delta \cos \tilde \alpha \cdot \Delta \alpha - \sin \tilde \delta \sin \tilde \alpha \cdot \Delta \delta} \cr {\cos \tilde \delta \cdot \Delta \delta} \cr}} \right]$$

and defined as:

(27)$$\Delta {\bi n} = \left[ {\matrix{ { - \cos \tilde \delta \sin \tilde \alpha \cdot \Delta \alpha - \sin \tilde \delta \cos \tilde \alpha \cdot \Delta \delta} \cr {\cos \tilde \delta \cos \tilde \alpha \cdot \Delta \alpha - \sin \tilde \delta \sin \tilde \alpha \cdot \Delta \delta} \cr {\cos \tilde \delta \cdot \Delta \delta} \cr}} \right]$$

Therefore:

(28)$${\bi n} = {\bi \tilde n} + \Delta {\bi n}$$

Substituting Equation (27) into Equation (4), the measurement equation with direction error is given as:

(29)$$c\Delta t_{sc}^{ij} = ({\bi n} + \Delta {\bi n}) \cdot \Delta {\bi r}_{sc}^{ij} = {\bi n} \cdot \Delta {\bi r}_{sc}^{ij} + \Delta {\bi n} \cdot \Delta {\bi r}_{sc}^{ij} $$

where $\Delta {\bi n} \cdot \Delta {\bi r}_{sc}^{ij} $ represents the effect of pulsar direction error. Given that the baseline length is short, the measurement error is expected to be small. Assuming that the pulsar direction error is 1 milli-arcsecond (mas), corresponding to the current level of pulsar position observation accuracy (Hanson, Reference Hanson1996) and the baseline length is 60 000 km, the measurement error attributable to the pulsar direction error is only 0·3 m, compared to the measurement error of a hundred meters using absolute pulsar-based measurement (Sheikh, Reference Sheikh and Pines2006). The relative pulsar-based measurement largely weakens the pulsar direction error for the measurement error.