Autonomous Positioning of Satellite Constellations via X-ray Pulsar Measurements

Abstract

A novel autonomous positioning approach based on X-ray pulsars is proposed in this paper. First, the principles of the pulsar–based measurement model and the inter-satellite range model in the autonomous positioning of satellite constellations are presented. The observability of the pulsar-based measurement model is then shown. Second, the autonomous positioning algorithms, including measurement models and orbital dynamics models, are formulated using an unscented Kalman filter to estimate the position vectors of satellites. Finally, the feasibility of the proposed measurement scheme compared with an inter-satellite range scheme is illustrated by numerical simulation. The results show that the proposed approach can keep the satellite state convergent, and achieve position accuracies of 2 m. The proposed scheme provides a promising approach for autonomous absolute positioning of constellation systems by using X-ray pulsars.

1. INTRODUCTION

Traditional methods of orbit determination in the three Global Navigation Satellite Systems (GNSS) rely extensively on ground-based operations (Kintner and Ledvina, 2005; Beutler and Moore, 2009). This working mode not only increases ground stations' burden and maintenance expenses, but also lacks long-term autonomous positioning capacity. Thus, autonomous positioning is a direction for the development of satellite constellations. Autonomous positioning mode implies that satellites confirm their own positions and maintain their time reference by using their own measurements without any support from ground stations for a long period.

To develop autonomous positioning technology for satellite constellations, a number of methods have been proposed. The inter-satellite range (Ananda et al., 1990; Zhao et al., 2011) method, which is used in the Global Positioning System (GPS) autonomous positioning mode, is a common method for the automatic positioning of satellite constellations. However given the lack of an absolute reference, this method can only guarantee the accuracy of the relative position but cannot ensure the correctness of constellation rotation (Menn and Bernstein, 1994) in the inertial frame. Inter-satellite range is unobservable for two error types, namely, orbital initial error and normal forces. Based on inter-satellite range measurement, the constellation's orientation measurement (Cai et al., 2006) is considered by using star sensors to resolve the problem of constellation rotation. However, this method requires star sensors with a high accuracy, far beyond current hardware capabilities. An anchor station (Wen and Zhu, 2009) on the Earth is also applied to estimate the position of satellites. However, a system that obtains support from the Earth cannot truly be considered an autonomous positioning system.

To achieve high-accuracy autonomous positioning, we propose a novel positioning method based on X-ray pulsar and inter-satellite range measurements for constellations. Pulsars are rapidly rotating, highly magnetized neutron stars that produce stable and predictable signatures (Sheikh and Pines, 2006; Bernhardt et al., 2011). Pulsars can provide an external spatial reference that differs from the inter-satellite range measurement. Moreover, based on X-ray pulsars, the position estimation accuracy of relative navigation is higher than estimation accuracy of absolute navigation (Emadzadeh et al., 2007). Given these special characteristics, X-ray pulsars are potential candidates for use in autonomous positioning systems. Therefore, we adopt a relative measurement based on X-ray pulsars to acquire the orientation information of satellite constellations. This measurement is then combined with inter-satellite range measurement that ascertains the relative position, to estimate the absolute position of the satellites in the constellation.

The outline of this paper is as follows: Section 2 presents the autonomous positioning measurement models as well as their mathematical formulations. Section 3 discusses the total autonomous positioning algorithms; the Unscented Kalman Filter (UKF) is utilized for its efficient performance in nonlinear estimation. Section 4 outlines the numerical simulations of these algorithms and presents an analysis of error factor.

2. AUTONOMOUS POSITIONING MEASUREMENT MODELS

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., 2006):

(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., 2006); 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., 2006);

• θ 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, 2007; Zhang, 2005; Liu, 2001).

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, 2001).

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. AUTONOMOUS POSITIONING ALGORITHMS

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, 2011).

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, 1997) 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. SIMULATION RESULTS AND ANALYSIS

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 2007) 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, 1996) 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, 2006). The relative pulsar-based measurement largely weakens the pulsar direction error for the measurement error.

4. CONCLUSIONS

An autonomous positioning approach for satellite constellations based on X-ray pulsars is proposed. This approach, using baselines between two satellites in the constellation, builds the system measurement model by combining relative pulsar-based and inter-satellite range measurements. It is shown that the satellites' three orbit elements, Ω, i, ω, which cannot be determined from inter-satellite range measurement, can be observed by using relative pulsar-based measurement. The absolute position of all satellites in the constellation can be determined efficiently by the above measurement models. Compared to the conventional method using absolute pulsar-based measurements, the proposed method significantly reduces the effect of pulsar direction errors on the satellite positioning performance.