Distributed Orbit Determination for Global Navigation Satellite System with Inter-Satellite Link

To keep the global navigation satellite system functional during extreme conditions, it is a trend to employ autonomous navigation technology with inter-satellite link. As in the newly built BeiDou system (BDS-3) equipped with Ka-band inter-satellite links, every individual satellite has the ability of communicating and measuring distances among each other. The system also has less dependence on the ground stations and improved navigation performance. Because of the huge amount of measurement data, the centralized data processing algorithm for orbit determination is suggested to be replaced by a distributed one in which each satellite in the constellation is required to finish a partial computation task. In the present paper, the balanced extended Kalman filter algorithm for distributed orbit determination is proposed and compared with the whole-constellation centralized extended Kalman filter, the iterative cascade extended Kalman filter, and the increasing measurement covariance extended Kalman filter. The proposed method demands a lower computation power; however, it yields results with a relatively good accuracy.


Introduction
For the global navigation satellite system (GNSS), the master control station (MCS) currently collects the satellite to monitor station measurement data, estimates the satellite ephemeris and clock offsets, and generates a time stream of navigation messages. The messages are then uploaded to satellites by ground antennas to broadcast to the user community [1]. However, the MCS as well as the other ground-based segments including monitor stations and ground antennas have the risk of destruction during a warfare or natural disaster. This is the case especially for the monitor stations which are distributed globally for increasing the accuracy of satellite orbit determination [2]. In order to enhance the viability of satellite navigation systems under the potentially fatal conditions, as early as in the 1980s, autonomous navigation techniques using inter-satellite link (ISL) measurements without support from the MCS were investigated for the global positioning system (GPS) [3]. If the ISL measurement was the only source for orbit determination and time synchronization, the datum mark would be insufficient [4]. This problem can be addressed by setting up a few ground anchorage stations (GASs) which provide a reference coordinate system and a time system [5]. Combining both the ISL and satellite-to-GAS measurements, the autonomous navigation system has several features: firstly, data processing will be completed by satellite onboard computers rather than the MCS; secondly,

Overview of Orbit Determination Algorithms
Ananda [7] proposed the framework of an autonomous navigation system without a MCS for the first time and validated the system by simulations. Rajan [3] introduced various autonomous navigation algorithms and presented preliminary on-orbit experiment results.
In the designing of the GPS Block IIR, the ISLs, programmable microprocessors, and redundant management were carried out. The following two major features are critical to ensure high precision autonomous navigation [3,7]: • the ISL communications and measurements in ultra high frequency (UHF) band; • a high-precision autonomous navigation algorithm which is adapted to the computing capacity of the satellite on-board computers.
A time division multiple access (TDMA) system, which has two frames, was employed for ranging and data transmission. In a ranging frame, each satellite is assigned a 1.5-s slot to make pseudo-range (PR) measurements with the visible satellites. Two frequencies were used in the measurement for ionospheric delay corrections. In a data frame, a 1.5-s slot was appointed to each satellite to transmit data that includes the PR measurements, the estimated satellite state vector, and the associated covariance matrix to all visible satellites.
The GPS Block IIF follows the design of the GPS Block IIR and improved the performance of ISL measurements and on-board data processing. Without contacting with the ground system, Block IIF can operate about 60 days in the autonomous navigation mode and provides navigation messages which are corrected by ISL measurements with a 3-m user range error (URE) (URE is a root mean square value and does not consider the impact of the polar motion and UT1). However, it is difficult to establish a precise prediction system because of the irregular polar motion and uncertain UT1. Therefore, in autonomous navigation mode, the URE is far greater than 3 m after a 60-day duration [8].
For GPS III, each satellite in the constellation will have the ability of ISL measurements and communications. It is designed that once there are enough satellites in the orbit, the GPS constellation will be able to operate autonomously in wartimes, but currently it is still under investigation [9].
Galileo navigation system is also planned to employ an autonomous navigation algorithm based on ISLs [10]. The spatial orientation problem is solved by the combination of the ISL and satellite-to-ground measurements. The simulation for the Galileo system showed that URE is on the level of decimeters [10], which is better than that of GPS.
For the distributed navigation algorithms, several techniques including the iterative cascade extended Kalman filter (ICEKF), increased measurement covariance EKF (IMCEKF), and Schmidt-Kalman Filter (SKF) are discussed by Schmidt, Park, and Ferguson [11][12][13]. The ICEKF is employed to processes a large number of space-borne GPS measurement data and a small quantity of ISL measurement data for low-Earth orbit formation flying satellites. In this method, the computation process will iterate for 3 to 4 times for convergence, and a good orbit accuracy is presented. For the method of IMCEKF, amendments are made based on ICEKF and it presents a better performance. On the other hand, the SKF yields results with less accuracy compared to IMCEKF while processing the GNSS ISL measurement data [4]. Recently, the International Association of Geodesy (IAG) initiated a GPS Dancer project which develops a distributed data processing algorithm to analyze the precision of GPS [14].
In China, distributed orbit determination and time synchronization algorithms based on ISL measurements were studied [15][16][17][18][19][20][21]. Distributed autonomous ephemeris updating was discussed for navigation satellites [22]. There are also many studies on developing higher precision orbit determination methods for new BeiDou (BDS-3) experimental satellites in the literature [23][24][25][26]. Currently, the developments for autonomous orbit determination, time synchronization, and autonomous operation and management are being widely investigated. More progress on designing an efficient distributed data processing algorithm, however, needs to be made.

Equations for Measurement
The position vectors, T , velocity vectors, T , and dynamic parameter vector, x i D , constitute the state vector which needs to be estimated: where the superscript T denotes the transpose of a matrix, and the superscript i denotes that it is for the ith satellite. The reference values for the state vector are stored in X i * , and the improving values for the state vector are written in: After correcting the hardware delay, ionospheric delay, relativistic effect, multi-path effect, and the antenna phase center offset [4], two one-way ISL PRs between the ith satellite and jth satellite need to be translated into the same measurement epoch (e.g., ranging frame epoch t). The PR equations are: where, c is the speed of light, δt i (t) and δt j (t) are clock errors for the ith satellite and the jth satellite, respectively, ε i→j (t) and ε j→i (t) are the measurement errors,  (3) and (4), the distance measurement equation that only contains orbit parameters is derived as: (3), the time measurement equation that only contains clock error parameters is deduced as: Following the steps above, the distance measurements and the clock bias measurements are decoupled. The orbit ephemeris and clock offsets can therefore be estimated independently.
Linearizing Equation (5) with Taylor expansion at the reference state vector X i * and X j * yields: After that, Equation (7) is converted into a linear measurement equation: where z ij (t) = ρ ij (t) − d X i * , X j * is the innovation. H i and H j are the measurement matrices.
Similarly, GAS can be considered as a pseudo-satellite. The PR needs to be corrected by an extra tropospheric delay compared to a normal satellite. The distance measurement equation that contains orbit parameters between the gth GAS and ith satellite is derived as: In the Equation (11), the array of state parameters of GAS X g (t) is known; only the state vector of the ith satellite is unknown. The reference ground coordinate is hence introduced into the satellite state by Equation (11), which overcomes the lack of the datum mark in data processing while only ISL measurements are utilized. The current linearized measurement equation, which is similar to Equation (8), becomes:

Equations for Motion
Satellites can be affected by a variety of factors when operated in orbit. For navigation satellites, only the gravitational forces, solar radiation pressure, and relativistic effects are considered [27]. The gravitational forces include the attractions from the earth, the moon, and other planets in the solar system. The dynamic equation for the ith satellite can be written as: where f c is a continuous function, w i is the system disturbances that have the following properties: E[·] denotes the expected value, and Q i (t) is a covariance matrix which is symmetric, non-negative, and definite. Here, the system disturbances are simulated by Gaussian white noise. Because the continuous function f c and the system disturbance w i are not coupled with each other, Equation (13) can be written as: In which G = 0 I 0 T is the coefficient matrix, I is the identity matrix, and 0 is the zero matrix.
Equation (15) is then linearized by Taylor expansion at the reference state vector X i * : From the equation above, the state increments are then derived as: , and F(t) is the dynamic partial derivative matrix [6]: Equation (17) is the state equation of the stochastic linear continuous systems and its general solution is: (19) in which Φ i (t, t 0 ) is the system state transition matrix and the solution of the following equations: where I is the identity matrix with the same dimensions as dynamics matrix, F(t).
The state transition matrix Φ i (t, t 0 ) has the following features: In the actual computation process, discretization needs to be implemented for Equation (19).
In a sampling interval from t k−1 to t k , the white noise w i k−1 (τ) can be considered as a constant. The integral coefficient denotes that: where Φ i k−1 denotes the state transition matrix from t k−1 to t k . According to Equation (24), the predicted state covariance matrix is:

Whole-Constellation Centralized Extended Kalman Filter
The whole-constellation centralized extended Kalman filter (WCCEKF) is one of the centralized data processing methods. According to this method, a main satellite and a back-up satellite are assigned to complete the task of data processing. The other satellites in the constellation need to send their measurement data, state vectors, and corresponding covariance matrices to the main satellite for orbit determinations.
For all the satellites in the constellation, the states and corresponding improving values from Equations (1) and (2) are collected and stored in a state vector X k and an improving values vector δx k [6,28], as in where n is the number of satellites in the system. In this way, the state equation for all the satellites can be obtained through Equation (24).
The state transition matrix Φ k and integral coefficient matrix G k are diagonal matrices and can be expressed as: The state noise vector, w k−1 , stores the noise for all the satellites in the constellation, and it has the statistical characteristics as follows: The measurement equation, which is a combination of Equations (8) and (12), is then derived as: where , and ε k = ε ij (t k ) for the ith satellite and jth satellite with ISL; , and ε k = ε ig (t k ) for GAS-to-ith satellite measurements. Next, the measurement covariance matrix R k = E(ε k ε T k ), the initial state vector , and the initial state vector covariance P 0 = E[X * 0 X * T 0 ] are defined. Finally, the method of WCCEKF that combines the satellite-to-satellite measurements and satellite-to-GAS measurements can be expressed as [29]: The dimension of state vector X k is where D i is the number of dynamic parameters for the ith satellite.
In the method of WCCEKF, each satellite is correlated with the other satellites through the state vector covariance matrix which has the dimension of N W × N W . Furthermore, matrix (H T k P k H k + R k ) with the dimension of m × m (m is the dimension of the measurement vector) needs to be inversed during the process, and a huge computation amount is expected. The computation amount for a process of WCCEKF is: If the WCCEKF algorithm is employed, the on-board computer of the main satellite would need to process all the ISL measurement and satellite-to-GAS measurement data to finish the task of orbit determination and navigation message generation for satellites in the constellation. Due to the huge computation amount and great complexity of communication, WCCEKF is difficult to be implemented in a satellite constellation with limited on-board computation ability.
In addition, the WCCEKF is also vulnerable. Once the main satellite and its backup satellite failed, the entire navigation constellation would stop working. To avoid the drawbacks in the WCCEKF method, many researches nowadays are focusing on developing a distributed data processing algorithm.

Distributed Orbit Determination
For distributed orbit determination algorithms based on ISL, data processing is assigned to each satellite. In this process, each satellite collects the ISL measurement data with respect to its visible satellite and estimates the self-related state vectors.

Reduced-Order Iterative Cascade EKF
For ith and jth satellites with ISL measurements, the iterative cascade EKF (ICEKF) [12] assumes that the state vector X j of the jth satellite is known. Thus, the measurement equation, which is similar to Equation (34), is derived as: For ISL measurement, the innovation is For satellite-to-GAS measurement, the innovation is z i k = z ig k (t k ), the measurement error is ε i k = ε ig (t k ), and measurement covariance matrix is GAS . An initial state Results from the method of ICEKF that combines the ISL measurement and satellite-to-GAS measurement can be obtained from: In this way, only the local state vector related to the ith satellite itself is included in the measurement equation. The dimension of state vector X i k is: The computation amount for the ICEKF algorithm is: As a result, the computational complexity is greatly reduced. However, the state vector of the ith satellite is only correlated with the measurement of itself and this method must be referred to as a reduced-order suboptimal filter. In order to improve the filtering accuracy, a common approach is to iterate the process above until convergence. In a data frame, after receiving the state vectors of the other visible satellites, the ith satellite updates its own state vectors and covariance matrix by Equations (42)-(47). Other satellites do the same process in turn and iterate until the state vector of each satellite is converged.
However, the method of ICEKF assumes that the state vectors of the other satellites have no errors, but this is not the case. Therefore, the method of ICEKF needs an uncertain number of iterations to approach convergence. In a constellation with a large number of satellites, reaching convergence could be time-consuming [12,13].

Reduced-Order Increased Measurement Covariance EKF
To accelerate the data-processing in the ICEKF method, the reduced-order increased measurement covariance EKF (IMCEFK) [12] is carried out. This method includes the error of the state vectors of the jth satellite into the measurement covariance matrix R i k between the ith satellite and the jth satellite. If we regenerate the ISL measurement error in Equation (42) , then the corresponding measurement covariance matrix is: where P j k is the state vector covariance matrix of the jth satellite. Equation (50) implies that ISL measurements contain not only the measurement errors, but also the jth satellite state vector error; thus, the measurement covariance matrix is assembled as: where the subscript assembled indicates that it is an assembled measurement covariance matrix. Next, the R i k in Equation (44) is replaced by R i k amp , and the gain matrix becomes: After repeating the steps in Equations (42), (43), (51), (45)-(47), the orbits for ith satellites are determined. In this way, a reduction of the number of iterations is expected. The computation amount of IMCEKF is: In some situations, iteration may not be required [18]. To summarize, compared to ICEKF, the IMCEKF is a reduced-order approach and needs to transmit not only the local state vector but also its covariance matrix to the other satellites.

Balanced Extended Kalman Filter
In one computation cycle, the ICEKF and IMCEKF algorithm only improve the state vector on one end of the ISL. To increase the efficiency and accuracy, the balanced extended Kalman filter (BEKF) is proposed. For the ith and jth satellite, the satellites' state vectors on the both ends of the ISL can be improved simultaneously. To keep the balance of the accuracy increments on both satellites, the improving state vectors should be adjusted by: With the constraint of Equation (54), the BEKF can be derived from Equations (8), (12), and (24), and it can be completed by the following steps: The dimension of state vectors X i k and X j k of the ith and jth satellites, respectively, is: The computation amount of BEKF is: The method has the following features: 1.
The method of BEKF collects the data of ISL measurements, satellite-to-GAS measurements (if they are available), and the satellite state vectors as well as their covariance matrices on both ends of the ISL. After the calculation of the BEKF algorithm, the improved state vectors and their covariance matrices are sent to the other visible satellites. The BEKF method modifies the denominator of the gain matrix K which is similar to the method of IMCEKF. Furthermore, it modifies the gain matrix by a factor of (I − M C ). Therefore, the BEKF algorithm is expected to yield results with higher precision.

2.
It seems that BEKF requires more ISL processes than the other EKFs. In fact, the state vectors and their covariance matrices on both ends are improved at the same time. It is unnecessary to repeat the ISL process for the same two satellites. The computation load of BEKF is similar to that of IMCEKF.

3.
The iteration process that is implemented in the ICEFK algorithm is not required in the BEKF method to achieve high accuracy. 4.
The improving state vectors are balanced in such a way that the satellite with lower-state precision will undergo more increments in accuracy while the satellite with higher-state precision will have fewer adjustments.

5.
Compared to the other EKFs, in Equation (64), M C P i k 0 0 P j k is subtracted from the state vectors' covariance matrices of the two satellites. Therefore, the values in the matrices are reduced and the accuracy of the state vectors is improved.

6.
The two satellites are correlated byP ij k andP ji k in Equation (65); however, these two matrices have to be ignored in this distributed filter. As a result, the current method should be categorized as a reduced-order sub-optimal orbit determination method.

Simulations and Analyses
In order to compare the performance of the abovementioned methods, navigation constellation simulations were carried out with the parameters of Walker 24/3/2:55 • , 22,116 km [30]. The dynamic model applied to satellite orbit was: (1) the earth's gravitational effects of 70 × 70, (2) the lunar, solar, and other planetary gravitational perturbations, (3) the solar radiation pressure, and (4) the other general relativistic forces.
The eighth-order Runge-Kutta method was employed for orbit integration. The IERS96 model was adopted for the Earth orientation parameters [31]. Besides, the TDMA mode was adopted in ISL with measurement frames and data transmitting frames. The total error of Ka-band ISL PR was 0.5 m (1σ). To avoid ground atmospheric disturbance, the ISLs with a vertical distance of less than 1000 km to the Earth surface were not considered. Eight GASs that are located in Xiamen, Kashi, Beijing, Lhasa, Sanya, Urumqi, Jiamusi, and Xi'an in China were set up in the simulation. With a minimum elevation of 10 • , the Hopfield -Marini model [1] were employed in tropospheric delay correction for the satellite-to-GAS PR which had a total error of 3 m (1σ).
The impact of complex factors was not considered in the simulations. In addition, the ISL PR measurement noise was assumed to be normally distributed without pollution, to have better comparisons for the different algorithms.
The orbit determination simulations were carried out in two steps: (1) An analytical orbit was generated and the corresponding ISL and satellite-to-GAS PRs were calculated. (2) Using the abovementioned PRs, the satellite orbits were calculated by the different methods and compared with the analytical orbit to find out orbit determination precisions.
The position error, radial error, along track error, and cross track error versus time normalized by day for satellite SV-01 computed from different methods are presented in Figures 1-4. It should be noted that plots of the errors from the other satellites in the constellation are excluded since the errors are similar to those of satellite SV-01.
(4) the other general relativistic forces.
The eighth-order Runge-Kutta method was employed for orbit integration. The IERS96 model was adopted for the Earth orientation parameters [31]. Besides, the TDMA mode was adopted in ISL with measurement frames and data transmitting frames. The total error of Ka-band ISL PR was 0.5 m (1 σ ).
To avoid ground atmospheric disturbance, the ISLs with a vertical distance of less than 1000 km to the Earth surface were not considered. Eight GASs that are located in Xiamen, Kashi, Beijing, Lhasa, Sanya, Urumqi, Jiamusi, and Xi'an in China were set up in the simulation. With a minimum elevation of 10°, the Hopfield -Marini model [1] were employed in tropospheric delay correction for the satellite-to-GAS PR which had a total error of 3 m (1 σ ).
The impact of complex factors was not considered in the simulations. In addition, the ISL PR measurement noise was assumed to be normally distributed without pollution, to have better comparisons for the different algorithms.
The orbit determination simulations were carried out in two steps: (1) An analytical orbit was generated and the corresponding ISL and satellite-to-GAS PRs were calculated. (2) Using the abovementioned PRs, the satellite orbits were calculated by the different methods and compared with the analytical orbit to find out orbit determination precisions.
The position error, radial error, along track error, and cross track error versus time normalized by day for satellite SV-01 computed from different methods are presented in Figures 1-4. It should be noted that plots of the errors from the other satellites in the constellation are excluded since the errors are similar to those of satellite SV-01.      For each method, the orbit determination errors tended to oscillate steadily after poor initial results. However, differences can be observed among the four algorithms. To have quantitative comparisons, the root mean square (RMS) of different errors for all the satellites in the constellation was calculated for the stable section. The average RMSs were then obtained by averaging the data from all the satellites and are plotted in Figure 5. It is worth pointing out that the average RMS of position error for the WCCEKF, ICEKF, IMCEKF, and BEKF algorithms was around 1.6 m, 4.5 m, 2.9 m, and 1.9 m, respectively.    For each method, the orbit determination errors tended to oscillate steadily after poor initial results. However, differences can be observed among the four algorithms. To have quantitative comparisons, the root mean square (RMS) of different errors for all the satellites in the constellation was calculated for the stable section. The average RMSs were then obtained by averaging the data from all the satellites and are plotted in Figure 5. It is worth pointing out that the average RMS of position error for the WCCEKF, ICEKF, IMCEKF, and BEKF algorithms was around 1.6 m, 4.5 m, 2.9 m, and 1.9 m, respectively.  For each method, the orbit determination errors tended to oscillate steadily after poor initial results. However, differences can be observed among the four algorithms. To have quantitative comparisons, the root mean square (RMS) of different errors for all the satellites in the constellation was calculated for the stable section. The average RMSs were then obtained by averaging the data from all the satellites and are plotted in Figure 5. It is worth pointing out that the average RMS of position error for the WCCEKF, ICEKF, IMCEKF, and BEKF algorithms was around 1.6 m, 4.5 m, 2.9 m, and 1.9 m, respectively. Next, the computation amounts for the different methods are summarized in Table 1 and visualized in Figure 6.

Methods
Computation Amount (FLOPS: Floating-Point Operations Per Second) Next, the computation amounts for the different methods are summarized in Table 1 and visualized in Figure 6.

Methods
Computation Amount (FLOPS: Floating-Point Operations Per Second)   Figure 6. Computation amounts.
The WCCEKF algorithm yielded optimal results with the highest precision. Among the three distributed orbit determination algorithms (ICEKF, IMCEKF, and BEKF), the highest orbit estimation precision was observed in the method of BEKF. When it comes to computation amounts for the four methods in Table 1 and Figure 6, the largest calculation amount was required in the WNCEKF algorithm, which achieved the best orbit accuracy. In the methods of ICEKF, IMCEKF, and BEKF with sub-optimal distributed orbit determination, significant reductions of the computation amount were observed when compared to the method of WNCEKF.
Finally, the performances of the WCCEKF, ICEKF, IMCEKF, and BEKF algorithms are summarized and provided in Table 2, considering different aspects.  The WCCEKF algorithm yielded optimal results with the highest precision. Among the three distributed orbit determination algorithms (ICEKF, IMCEKF, and BEKF), the highest orbit estimation precision was observed in the method of BEKF. When it comes to computation amounts for the four methods in Table 1 and Figure 6, the largest calculation amount was required in the WNCEKF algorithm, which achieved the best orbit accuracy. In the methods of ICEKF, IMCEKF, and BEKF with sub-optimal distributed orbit determination, significant reductions of the computation amount were observed when compared to the method of WNCEKF.
Finally, the performances of the WCCEKF, ICEKF, IMCEKF, and BEKF algorithms are summarized and provided in Table 2, considering different aspects.

Conclusions
The fundamental theory of satellite orbit determination for autonomous navigation is introduced. Four algorithms for autonomous navigation with onboard data processing, i.e., whole-constellation centralized extended Kalman filter (WCCEKF), iterative cascade extended Kalman filter (ICEKF), increased measurement covariance extended Kalman filter (IMCEKF), and balanced extended Kalman filter (BEKF), are illustrated. The WCCEFK technique processes the measurement data for orbit determination by a main satellite while the other three algorithms distribute the computation on every satellite in the constellation. The simulation results show that the method WCCEKF has the optimal orbit determination accuracy among the four methods but demands the largest computation loads. Similar computation amounts are observed in the three distributed algorithms. Compared to the ICEKF technique, covariance matrices of the other satellites are absorbed into the measurement covariance matrix in the method of IMCEKF. As a result, a smaller orbit error is observed in the IMCEKF algorithm than in the ICEKF. Since the BEKF method estimates the state vectors of the satellites on both ends of the ISL in a balanced mean and increases their accuracy simultaneously, this method yields the best results among the three distributed estimation algorithms. The BEKF can be considered as an appropriate distributed data processing algorithm for a GNSS.