Evaluation of BDS-3 Orbit Determination Strategies Using Ground-Tracking and Inter-Satellite Link Observation

Abstract

Dual one-way inter-satellite link (ISL) pseudoranges of BDS-3 satellites can be introduced as an additional measurement along with L-band pseudoranges and phases to improve the accuracy of precise orbit determination (POD). In the existing research, although the clock-free or geometry-free ISL observables are derived from the raw two one-way pseudoranges, only the clock-free observables are adopted for the ISL joint POD (Joint 1 POD) without considering the geometric-free observables. An improved joint (Joint 2 POD) strategy making full use of the clock-free and geometry-free observables is applied in this contribution. The orbits of ground-only POD, ISL-only POD, Joint 1 POD, and Joint 2 POD are comprehensively compared by the orbit overlapping differences, the Satellite Laser Ranging (SLR) residuals, and the characteristics of the satellite clock offsets estimated simultaneously. The comparisons prove that the performance of the Joint 2 POD strategy is better than that of the other three POD strategies regardless of the types of satellites. Moreover, this paper discusses ISL’s contribution to the station selection strategy in terms of the number and distribution. The experimental results show that, when there are more than 20 stations, each additional 10 stations contributes to a maximum of 7.5%, 3.9%, and 2.8% improvement on MEO, IGSO, and GEO satellites 3D accuracy, respectively. When the number of stations reaches 50, the precise orbits achieve similar accuracy to the results using 80 stations. In addition, after adding ISL data, the orbits estimated using 10 regional stations and 10 global stations are greatly improved, and the accuracy between them is only 0.9 cm in 3D errors.

1. Introduction

The third-generation BeiDou Navigation Satellite System (BDS-3), which is currently being deployed, is a global satellite navigation system. It is the third step of China’s satellite navigation construction. The complete BDS-3 constellation will consist of 24 satellites in three orbital planes in medium Earth orbit (MEO) with an orbit inclination angle of 55°; 3 satellites in geostationary Earth orbit (GEO) located at longitudes 80°, 110.5°, and 140°; and 3 additional satellites in inclined geosynchronous orbits (IGSO) with an inclination of 55° over the Asia-Pacific region [1,2]. By the end of 2019, 28 BDS-3 satellites have been launched into space.

BDS-3 adds three new frequencies to the regional BeiDou satellite navigation system (BDS-2) [3], and three types of BDS-3 satellites are given different functions. Due to the geostationary characteristics of GEO satellites and the resulting lower orbit accuracy, BDS-3 GEO satellites are mainly used to provide some featured services, such as satellite-based augmentation (SBAS), short message communication (SMC), and precise point positioning (PPP) [1]. The rest of the 24 MEO satellites and 3 IGSO satellites play a key role in providing the radio navigation satellite service (RNSS). However, the other side of the coin is that this hybrid constellation challenges precision orbit determination (POD) because MEO satellites require tracking stations distributed globally and IGSO/GEO satellites are only tracked by stations in the Asia-Pacific region.

The normal operation of navigation satellites is accomplished by means of continuous tracking and control by the ground segment. To get rid of the dependence on the ground-based support, GPS has utilized an ultrahigh-frequency (UHF) inter-satellite link (ISL) technology to achieve autonomous navigation since the GPS IIR satellites [4,5,6,7,8]. Research shows that autonomous navigation can maintain navigation accuracy for six months without ground contact [8]. For BDS, the monitoring stations are distributed in mainland China for various reasons. Once the satellite loses the ground contact, the satellite will not be controlled, and thus orbit and clock accuracy will drop sharply. Therefore, the new BDS-3 satellites are equipped with ISL payload following GPS. This new payload enables intel-satellite communication and ranging by the Ka-band signal, which is suffered from less multipath and noise issues. ISL can be used to route the data from satellites out of the coverage region to satellites that are within the coverage region. Benefit from these new payloads, the satellite running outside mainland China can also be tracked and controlled indirectly. The operational control system (OCS) can make full use of both L band signals and Ka-band signal observations to generate the broadcast messages and to update messages [9].

Moreover, the ISL payloads provide BDS-3 with the potential for autonomous navigation. Many results of the BDS-3 experimental or BDS-3 satellites orbit determination with the real ISL measurements are published. Ren et al. discussed the orbit determination with prior orbit constraints [10]. Yang et al. presented a detailed evaluation of the combined orbit determination method, which indirectly achieves a 37–95% improvement in the orbit-only signal-in-space range error accuracy [11]. Tang et al. applied the ISL data from BDS experiment satellites on the centralized autonomous orbit determination [12]. Ren et al. analyzed the application of batch processing and extended Kalman filtering algorithm on BDS autonomous navigation [13]. Guo et al. proposed an onboard extended Kalman filter constrained with high predicted orbit inclination and longitude of ascending node to resolve the constellation rotation error issue [14].

Precise orbit is the key to high-precision positioning applications and related scientific research. The other advantage of ISL is that the accuracy of the BDS-3 satellite orbit can be improved by using both ISL and L-band measurements. Xie et al. used 16 global distributed stations and 6 China regional stations to verify the improvement of orbit determination accuracy when adding ISL measurement and finally achieved a great improvement of 42% and 83%, respectively [15]. Yang et al. explored that with the help of ISL, the 3D root mean square (RMS) of eight satellites orbit position computed using 10 regional stations was decreased from 1 to 0.5 m [16]. In general, orbit accuracy at the few decimeters level could be achieved by the different studies.

Nevertheless, there are still many aspects of BDS-3 ISL that have not been studied. For example, there are two ISL observation models applied for joint POD using ISL and ground-tracking data. One is to treat a pair of raw one-way ISL pseudoranges as two independent observables [17], and the other is to derive the ISL clock-free observables from a dual one-way pseudoranges [2,11,15,18]. These derived clock-free observables contain only geometric information, so the observation model excludes the satellite clock offset parameters. In the case of joint POD, it is unreasonable to use only the clock-free observables. Because the satellite clock offsets and the orbit in the radial component are highly correlated in the L-band POD results [19], the addition of ISL data can reduce the orbit impact on the estimation of satellites clock offsets and adjust the clock parameters. Besides, as more and more receivers are upgraded, more than 100 ground tracking stations have access to receive the signal broadcast by BDS-3 satellites. Considering the addition of a large amount of ISL data, there will be a lot of L-band redundant data for joint orbit determination, which cannot contribute to the accuracy, but reduce the calculation efficiency. On the other hand, since BDS-3 is a hybrid constellation, the distribution of stations required for the POD of satellites operating in different orbital planes is not the same. The addition of ISL measurements can optimize the geometry of the observation model. Therefore, further research is needed to explore whether the introduction of ISL to joint POD has requirements on the distribution of stations.

This article continues to analyze the issue of ISL’s contribution to the joint POD. An improved joint POD strategy using ground-tracking and ISL observations of the BDS-3 satellites is applied. This strategy makes full use of both clock-free and geometry-free ISL observables derived from the raw dual one-way ISL observations, which has not been used by previous studies. The satellite orbit parameters, as well as satellite clock offsets estimated at the same time, are evaluated and compared. Through this improved strategy, the ISL’s contribution to joint POD in the case of the selection of stations, including the distribution of stations and the number of stations, is analyzed separately for three types of BDS-3 satellites.

2. Methodology

The topology of the BDS-3 satellite links is a dynamic network. Based on a pre-designed timeslot scheduling, a satellite can establish a connection with another satellite in 3 s following a time division multiple access (TDMA) scheme. Since the round-trip measurement period is short, the dual one-way pseudoranges are usually normalized to the same reference time to derive geometry-free or clock-free observables. Once the derived observables are obtained, joint orbit determination can be performed with L-band code and phase measurements.

2.1. Dynamic Inter-Satellite Network

ISL measurement is a dual one-way observation generated by a pair of satellites. The two satellites are each other’s transmitter and receiver. The phased array antenna equipped on BDS-3 satellites can realize the beam scanning in a wide space range of −60° to 60° [15]. The pointing agility and fast-tracking characteristics of phased array antenna make it possible to complete dual unidirectional ranging in a short time of 3 s, which can establish a dynamic inter-satellite network [20]. Phased array antennas work following predesigned timeslot schedule by ground control personnel. The timeslot schedule records the chain building information for a week, including when and how one satellite connects with the other satellites.

The timeslot schedule of the BDS-3 ISL takes 1 h as a pre-designed period. Each pair of satellites is assigned a short-time segment of 3 s to finish a dual one-way ranging. A forward link is established in the first 1.5 s, and a backward link is established in the second 1.5 s [20]. Therefore, a satellite can establish contact with up to 20 satellites or anchor stations in one minute without regard to visibility. In the following 59 min, the order and rules for establishing the links remain the same as in the first minute.

2.2. L-Band Observation Model

Both carrier phase measurements and the code measurements are adopted in the POD processing. Ionospheric-free (IF) combination is used to remove the first-order ionospheric effect. The L-band observation equations about IF phase measurements $P_{k,IF}^j$ and code measurements ${\Phi }_{k,IF}^j$ are expressed as

$$P_{k,IF}^j={\rho }_k^j-\mathrm{C}\cdot \left(d{t}^j-d{t}_k\right)+m_k^j\cdot T{r}_k+{\epsilon }_{k,p}^j,$$

(1)

$${\Phi }_{k,IF}^j={\rho }_k^j-\mathrm{C}\cdot \left(d{t}^j-d{t}_k\right)+m_k^j\cdot T{r}_k-{\lambda }_{IF}{N}_{k,IF}^j+{\epsilon }_{k,L}^j,$$

(2)

where ${\rho }_k^j$ is the distance from the satellite j to the L-band antenna phase centers of the receiver k; $d{t}^j$ and $d{t}_k$ are satellite and receiver clock offset, respectively; $T{r}_k$ is the non-dispersive tropospheric delay in zenith direction, which can be mapped to satellite-to-receiver direction by multiply the mapping function $m_k^j$; C is the light speed; ${\lambda }_{IF}$ is the wavelength of the IF observation and $N_{k,IF}^j$ is the corresponding ambiguity; and ${\epsilon }_{k,p}^j$ and ${\epsilon }_{k,L}^j$ represent the receiver noise of code and phase measurement, respectively.

2.3. Ka-Band Observation Model

As a pair of one-way pseudoranges, ISL measurement can be used alone or along with L-band data to determine satellite orbits and clock offsets. Suppose that satellite A receives the Ka-band signal transmitted by satellite B at time t₁ and generates the observation ${\rho }_{BA}\left(t_1\right)$ and satellite B receives the Ka-band signal transmitted by satellite A at time $t_2$ and generates the observation ${\rho }_{AB}\left(t_2\right)$. Then, the Ka-band observation model can be formulated as

$$\begin{array}{c}{\rho }_{BA}\left(t_1\right)=\left|{\overrightarrow{r}}_A\left(t_1\right)-{\overrightarrow{r}}_B(t_1-\frac{{\rho }_{BA}}{\mathrm{C}})\right|+\mathrm{C}\cdot \left(d{t}_A\left(t_1\right)-d{t}_B(t_1-\frac{{\rho }_{BA}}{\mathrm{C}})\right)+\mathrm{C}\\ \hspace{1em}\hspace{1em}\cdot \left({\tau }_A^{Rcv}+{\tau }_B^{Send}\right)+\mathsf{\Delta }{\delta }_{BA}+{\epsilon }_{BA}\end{array}$$

(3)

$$\begin{array}{c}{\rho }_{AB}\left(t_2\right)=\left|{\overrightarrow{r}}_B\left(t_2\right)-{\overrightarrow{r}}_A(t_2-\frac{{\rho }_{AB}}{\mathrm{C}})\right|+\mathrm{C}\cdot \left(d{t}_B\left(t_2\right)-d{t}_A(t_2-\frac{{\rho }_{AB}}{\mathrm{C}})\right)+\mathrm{C}\\ \hspace{1em}\cdot \left({\tau }_B^{Rcv}+{\tau }_A^{Send}\right)+\mathsf{\Delta }{\delta }_{AB}+{\epsilon }_{AB}\end{array}$$

(4)

where ${\overrightarrow{r}}_A$ and ${\overrightarrow{r}}_B$ is the Ka-band satellite antenna phase center vectors in Cartesian coordinates at the certain observation epoch; $d{t}_A$ and $d{t}_B$ are satellite clock offset; ${\tau }_A^{Send}$ and ${\tau }_B^{Send}$ are the transmitting hardware delay; ${\tau }_A^{Rcv}$ and ${\tau }_B^{Rcv}$ are receiving hardware delay; C represents the speed of light; $\mathsf{\Delta }{\delta }_{AB}$ and $\mathsf{\Delta }{\delta }_{BA}$ refer to corrections that can be easily modeled for the measurement, such as relativistic effects and Ka-band phase center offsets; and ${\epsilon }_{BA}$ and ${\epsilon }_{AB}$ are the measurement noises.

We can treat them as two independent one-way observations and model them directly similar to the L-band pseudoranges [17], or we can transform a pair of one-way observations to two observations at the same time for processing, as

$${\rho }_{BA}\left(t_0\right)={\rho }_{BA}\left(t_1\right)+\mathsf{\Delta }{\rho }_{BA},$$

(5)

$${\rho }_{AB}\left(t_0\right)={\rho }_{AB}\left(t_2\right)+\mathsf{\Delta }{\rho }_{AB},$$

(6)

where

$${\rho }_{BA}\left(t_0\right)=|{\overrightarrow{r}}_A\left(t_0\right)-{\overrightarrow{r}}_B\left(t_0\right)|+\mathrm{C}\cdot \left(d{t}_A\left(t_0\right)-d{t}_B\left(t_0\right)\right)+\mathrm{C}\cdot \left({\tau }_A^{Rcv}+{\tau }_B^{Send}\right)+\mathsf{\Delta }{\delta }_{BA}+{\epsilon }_{BA},$$

(7)

$${\rho }_{AB}\left(t_0\right)=|{\overrightarrow{r}}_B\left(t_0\right)-{\overrightarrow{r}}_A\left(t_0\right)|+\mathrm{C}\cdot \left(d{t}_B\left(t_0\right)-d{t}_A\left(t_0\right)\right)+\mathrm{C}\cdot \left({\tau }_B^{Rcv}+{\tau }_A^{Send}\right)+\mathsf{\Delta }{\delta }_{AB}+{\epsilon }_{AB},$$

(8)

and the correction $\mathsf{\Delta }{\rho }_{BA}$ and $\mathsf{\Delta }{\rho }_{AB}$ can be computed with broadcast messages following:

$$\begin{array}{l}\mathsf{\Delta }{\rho }_{BA}={\rho }_{BA}\left(t_0\right)-{\rho }_{BA}\left(t_1\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}=|{\overrightarrow{r}}_A\left(t_0\right)-{\overrightarrow{r}}_B\left(t_0\right)|-\left|{\overrightarrow{r}}_A\left(t_1\right)-{\overrightarrow{r}}_B\left(t_1-\frac{{\rho }_{BA}}{\mathrm{C}}\right)\right|+\mathrm{C}\\ \hspace{1em}\hspace{1em}\hspace{1em}\cdot \left(d{t}_A\left(t_0\right)-d{t}_B\left(t_0\right)\right)-\mathrm{C}\cdot \left(d{t}_A\left(t_1\right)-d{t}_B(t_1-\frac{{\rho }_{BA}}{\mathrm{C}})\right)\end{array}$$

(9)

$$\begin{array}{l}\mathsf{\Delta }{\rho }_{AB}={\rho }_{AB}\left(t_0\right)-{\rho }_{AB}\left(t_2\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}=|{\overrightarrow{r}}_B\left(t_0\right)-{\overrightarrow{r}}_A\left(t_0\right)|-\left|{\overrightarrow{r}}_B\left(t_2\right)-{\overrightarrow{r}}_A\left(t_2-\frac{{\rho }_{AB}}{\mathrm{C}}\right)\right|+\mathrm{C}\\ \hspace{1em}\hspace{1em}\hspace{1em}\cdot \left(d{t}_B\left(t_0\right)-d{t}_A\left(t_0\right)\right)-\mathrm{C}\cdot \left(d{t}_B\left(t_2\right)-d{t}_A(t_2-\frac{{\rho }_{AB}}{\mathrm{C}})\right)\end{array}$$

(10)

The accuracy of $\mathsf{\Delta }{\rho }_{BA}$ and $\mathsf{\Delta }{\rho }_{AB}$ depends on the changing rate of satellite orbit and clock offsets computed by broadcast messages, and the accuracy loss is much smaller than the ISL residuals, which can be ignored in a short time [11]. Obviously, the sum and difference of Equations (7) and (8) are equal to

$$\frac{{\rho }_{BA}\left(t_0\right)+{\rho }_{AB}\left(t_0\right)}{2}=\left|{\overrightarrow{r}}_A\left(t_0\right)-{\overrightarrow{r}}_B\left(t_0\right)\right|+ \mathrm{C}\cdot \frac{{\tau }_A^{Rcv}+{\tau }_B^{Send}+{\tau }_B^{Rcv}+{\tau }_A^{Send}}{2}+\frac{{\epsilon }_{BA}+{\epsilon }_{AB}}{2}$$

(11)

$$\frac{{\rho }_{BA}\left(t_0\right)-{\rho }_{AB}\left(t_0\right)}{2}=\mathrm{C}\cdot \left(d{t}_A\left(t_0\right)-d{t}_B\left(t_0\right)\right)+\mathrm{C}\cdot \frac{{\tau }_A^{Rcv}+{\tau }_B^{Send}-{\tau }_B^{Rcv}-{\tau }_A^{Send}}{2}+\frac{{\epsilon }_{BA}-{\epsilon }_{AB}}{2}$$

(12)

Compared to the L-band observation model expressed by Equations (1) and (2), the clock offsets or geometry distance are removed in Equation (11) or Equation (12). This feature allows us to perform POD or clock offsets estimation separately. It is very meaningful for clock characteristic evaluation. In the past, the satellite clock offsets computed by the L-band data contained unmodeled orbit errors, which affected the evaluation results. Through the geometry-free Equation (12), the influence of the orbit error on the clock offsets will be greatly eliminated. In the previous studies about BDS-3 POD, only Equation (11) is adopted in Ka-band POD. Therefore, the information and functional relationship implied in Equation (12) are rudely discarded. The improved joint POD method in this paper makes full use of both geometry-free Equation (11) and clock-offsets-free Equation (12) to achieve joint orbit determination.

2.4. Processing Strategy

Four POD modes, namely ground-only POD, ISL-only POD, Joint 1 POD, and Joint 2 POD, are utilized to analyze the contribution of BDS-3 ISLs to orbit determination in this research. For ground-only POD, the code and carrier measurements in ionosphere-free combination expressed as Equations (1) and (2) are adopted. For ISL-only POD, only the clock-free pseudoranges described as Equation (11) are used. In particular, only in this mode, two ground anchor stations are added to reduce the constellation rotation resulting from the ISL-only POD [8,12,21]. For the Joint 1 POD, only ISL clock-free observations (Equation (11)) along with the IF code and carrier measurements (Equations (1) and (2)) are adopted. The improved joint POD applied in this contribution is expressed as Joint 2 POD, and the IF code and carrier measurements (Equations (1) and (2)) along with ISL geometry-free and clock-free observations (Equations (11) and (12)) are used simultaneously to determine the precise orbit and clock parameters as well as other relevant parameters. Some necessary correction models and settings of L-band and Ka-band observation are summarized in

Table 1. Processing strategy for the POD.

Correction	L-Band Model	Ka-Band Model
Observation type	Code and carrier observations of B1 and B3	Code observations
Measurements combination	Ionosphere-free combination	Geometry-free combination and clock-free combination
Elevation cutoff	10°	-
Sampling rate	300 s	300 s
Weight	Elevation-dependent weight factor with a standard deviation of 0.002 m and 0.2 m for raw phase and code, respectively	0.5 m for the ISL geometry-free and clock-free measurements and independent of satellite nadir angle
Tropospheric delay	Saastamoinen model [22], global mapping function [23], and two-hourly Zenith Tropospheric Delays (ZTD) without gradients.	Same as L-band model when considering the anchor stations in ISL-only POD, and no correction is needed for satellite-to-satellite link
Phase wind up	Corrected [24]	-
Satellite antenna phase center	Published on BDS official website (http://en.beidou.gov.cn/SYSTEMS/Officialdocument)	Provided by manufacturer
Tide displacement	Solid Earth tide, pole tide, and ocean tide loading; according to International Earth Rotation Service (IERS) Conventions 2010 model [25]	Same as L-band model when considering the anchor stations in ISL-only POD
Relativity effect	IERS Conventions 2010 model [25]
Geopotential	EIGEN_GL04C up to 12 × 12 degrees and orders [26]
N-body gravitation	Jet Propulsion Laboratory (JPL) DE405 ephemeris model
Solar radiation	ECOM five-parameter model [27]

.

It is important to emphasize the differences between the observables of L-band and Ka-band. First, the different types of observations have different weight settings. A prior standard deviation for the raw L-band phase, code, and Ka-band measurements are 0.002, 0.2, and 0.5 m, respectively. After that, the final weight of L-band observations needs to be adjusted according to the satellite elevation, while the weight of the Ka-band is independent of the satellite nadir angle. Second, although a pair of one-way ranging is completed at a 3-s interval, the observables are still transformed to the reference time with a sampling interval of 300 s, which is to maintain consistency with the L-band data to facilitate the estimation of the satellite clock offsets. However, due to the accuracy of broadcast ephemeris, the precision of the transformation (Equations (9) and (10)) is limited by the length of the extrapolation time of reference epoch in the broadcast ephemeris, thus only the ISL data in a 120-s time window centered on the reference time are utilized [15]. Moreover, for ISL data, the ionosphere delay is negligible because of the high frequency (23 GHz) of the Ka-band and the low electron density on the signal propagation path. Similarly, there is no need to consider tropospheric correction, if there are no ISL measurements with regard to anchor stations. Besides, the L-band antenna phase center offset (PCO) values are provided by the BDS official website while the Ka-band PCO values are provided by manufacturers.

Four groups of precise BDS-3 orbits are processed by revised PANDA software [28], which has an additional function of outlier detection and adjustment for ISL data. The procedure of the Joint 1/2 POD can be divided into three steps. First, the POD using ground code and phase observation is carried out. This step is mainly to eliminate the L-band abnormal data and detect carrier cycle slips. Then, the broadcast orbit is used as a reference for detecting outliers in the ISL data through Equation (11). Once the above two steps are completed, the derived ISL measurements will be added to the least squares to estimate the orbital and clock parameters together with L-band measurements. After a few iterations of least squares adjustment, the orbital parameters along with satellite clock and other estimated parameters can be obtained.

The general method to evaluate the internal consistency of satellite orbits is the direct comparison of positions from different orbit solutions, e.g., the last two days of a three-day arc with the first two days of the next three-day arc. The residuals are usually given in a satellite-centered orbit system in along-track, radial, and cross-track directions. Besides, the Satellite Laser Ranging (SLR) data can be used as an independent tool to evaluate orbit accuracy.

2.5. Data Collection

The experimental data of BDS-3 satellites come from a total of 80 MGEX and iGMAS ground stations and 190 pairs of satellite-to-satellite links during the period of day-of-year (DOY) 198–280 in 2019. The distribution of these stations is marked in Figure 1. Note that during this period only 21 BDS-3 satellites were in orbit, including 18 MEO satellites, 2 IGSO satellites (C38 and C39), and 1 GEO satellite (C59). However, the IGSO and GEO satellites are still in a test stage, during which the signal can only be decoded by the iGMAS receivers. The observation of satellite C39 started from DOY 227. A brief comparison of the amount of two types of measurements in BDS Week 715 (DOY 251–257) is plotted in Figure 2 as an example. Figure 2a shows the number of available one-way ISL measurements in pairs, where the outliers and single one-way observations that cannot form a dual one-way link are discarded. Figure 2b shows the number of available code/phase measurements after screening outliers and deleting short arcs less than 1200 s. It is obvious that there are many more ISL data than L-band measurements despite the different types of orbits. The number of L-band measurements of IGSO and GEO satellites is less than that of MEO satellites.

3. POD Results

In this section, the orbits estimated by the Joint 2 POD method are presented and compared with other POD methods, in order to analyze the ISL’s contribution to the orbit determination. Then, we discuss the enhancement of introducing ISL in the selection of L-band tracking stations, including the distribution of stations and the number of stations.

3.1. Comparison among Different POD Strategies

According to the four POD modes mentioned above, the precise orbits derived from three-day arcs are estimated. The two-day overlapping differences between two adjacent three-day solutions in along-track (A), radial (R), and cross-track (C) direction are used to evaluate the internal consistency of the orbits. The RMS results of MEO and IGSO satellites are presented in Figure 3, and the RMS value along with the errors of GEO satellites C59 are listed separately in

Table 2. RMS of two-day arc overlapping difference of ground-only orbit, ISL-only orbit, Joint 1 orbit, and Joint 2 orbit based on three-day arcs solutions according to the satellite type.

Type	Strategy	Along-Track (cm)	Cross-Track(cm)	Radial (cm)	3D (cm)
MEO	Ground-only orbit	5.8	3.1	1.4	6.7
	ISL-only orbit	13.7	14.0	2.1	19.8
	Joint 1 orbit	5.8	4.2	1.7	7.4
	Joint 2 orbit	5.3	3.9	1.5	6.7
IGSO	Ground-only orbit	17.6	10.3	4.7	20.9
	ISL-only orbit	21.0	22.2	4.3	30.9
	Joint 1 orbit	8.4	7.6	3.2	11.8
	Joint 2 orbit	8.6	9.0	3.1	12.8
GEO	Ground-only orbit	87.3	6.7	6.9	87.9
	ISL-only orbit	20.6	18.6	6.0	28.4
	Joint 1 orbit	14.2	9.2	5.7	17.9
	Joint 2 orbit	14.4	9.2	5.4	17.9

because they are much larger than those of MEO and IGSO satellites.

According to Figure 3 and Table 2, it can be found that the addition of ISL data has a positive effect on the orbit determination. For the ground-only POD, the overlapping accuracy of MEO satellites is similar to that of Joint 2 POD and higher than Joint 1 POD because there are currently enough globally distributed ground stations to observe MEO satellites. On the other hand, the overlapping errors as an indicator of internal coincidence precision reflect the internal consistency of the orbit estimates, not the true accuracy. Besides, the improper weight of ISL data may also affect the contribution of ISL data on joint orbit determination, which needs further research in the future. Different from MEO satellites, the IGSO satellites can only be observed by iGMAS stations located in the Asia-Pacific region during the evaluation period. Therefore, the accuracy of Joint 2 has been significantly improved by 38.8% and 43.5% after adding ISL data for satellites C38 and C39, respectively. As a GEO satellite, satellite C59 has the characteristic of geostationary and high orbit height, and it can only be tracked by ground stations near the Asia-Pacific region similar to IGSO satellites. With the help of ISL, the weak geometry of the GEO satellite is strengthened, and the number of measurements is increased. Compared with ground-only results, the Joint 1/2 orbit accuracy of satellite C59 is greatly increased by 79.6%. As for ISL-only POD results, the orbit accuracy for the GEO satellite is greatly improved compared to the orbits obtained by ground-only POD strategy. It is noted that the ISL-only POD errors of the whole constellation in the along-track and cross-track direction are larger than those of the ground-only POD, except for the along-track direction of GEO satellites, which may be caused by the weaker orbital orientation constrain from only two anchor stations [12].

For the two different joint POD results, the Joint 2 orbits show a radial improvement of 10%, 2.5%, and 1.4% for MEO, IGSO, and GEO satellites, respectively, compared to the Joint 1 orbit. As a result, the 3D accuracy of MEO has increased by 9%. However, in the case of IGSO and GEO satellites, the addition of ISL clock information increases the errors in the cross-track or along-track direction and even leads to a decrease of 8.5% in the 3D accuracy of the IGSO satellites. Considering the fact that the IGSO and GEO satellites are still in the test phase and there is no measurement observed from stations other than iGMAS stations, a further study on the orbit differences between the two joint POD method needs to be carried out in the future.

The SLR data can be used as a reliable independent tool to evaluate orbit accuracy. At present, there are only four MEO satellites that are regularly tracked by International Laser Ranging Service (ILRS) [29,30]. The residuals (observed minus computed) are presented in

Table 3. SLR residuals to BDS-3 MEO satellites derived from ground-only orbit, ISL-only orbit, Joint 1 orbit, and Joint 2 orbit based on three-day arcs.

Strategy	PRN	No. of Normal Points	Mean (cm)	STD (cm)	RMS (cm)
Ground-only orbit	C20	797	1.0	5.5	5.6
	C21	766	3.6	5.6	6.6
	C29	568	−3.3	3.8	5.0
	C30	572	−3.5	3.9	5.2
ISL-only orbit	C20	797	−0.2	7.8	7.8
	C21	766	2.5	7.6	8.0
	C29	568	−1.2	10.9	10.9
	C30	572	0.7	10.6	10.6
Joint 1 orbit	C20	797	−1.1	3.3	3.5
	C21	766	1.9	3.5	4.0
	C29	568	−3.5	3.5	5.0
	C30	572	−2.2	3.4	4.1
Joint 2 orbit	C20	797	−0.9	3.0	3.1
	C21	766	2.1	3.1	3.7
	C29	568	−3.0	2.8	4.1
	C30	572	−3.0	3.1	4.3

. It follows that the ISL-only orbits have the worst standard deviations (STDs) while the ground-only orbits have the largest biases in the SLR residuals. When using both L-band and Ka-band data, the SLR residuals are significantly reduced. The RMS of the Joint 2 results is 8% smaller than that of the Joint 1 results.

Since the orbit radial component is highly correlated with the clock offset parameters, the performance of the satellite clock offset estimated with the orbit parameters at the same time can also be used as an indicator to evaluate the accuracy of the orbits. Note that the ISL-only strategy does not solve the clock offset parameter. The clock overlapping differences of the two-day arc are computed based on two adjacent three-day arcs, and the RMS value is plotted in Figure 4. It can be found that the performance of the MEO satellite clock offsets became worse when applying the Joint 1 strategy, but the accuracy of clock offsets computed by the Joint 2 strategy shows a better performance than that of the ground-only strategy with an improvement of 27.0%, 49.5% and 122.3% for MEO, IGSO, and GEO satellites, respectively.

Existing studies have shown that the orbit error contained in the clock offsets estimated by L-band observations can be observed in the Allan variance [31], while the application of ISL-only measurements can reduce this error in clock offsets [32]. Further analysis of the stability characteristic of the satellite clock is shown in Figure 5. The overlapping Allan deviation values are computed by the satellite clock offsets that have been detrended by removing a daily quadratic trend. In the case of MEO satellites, the obtained clock stability of the two joint POD with ISL is improved compared to the ground-only clock. The bumps at around 10,000–20,000 s that are usually caused by the orbit errors lying in the ground-only clock offsets become much smoother in the results of the Joint 2 solutions. For the IGSO and GEO (C59) satellites, the stability derived from the Joint 2 clock is greatly improved compared to the other POD strategies. It is noted that the clock stability of satellite C59 obtained by the Joint 1 POD without geometry-free observables shows unexpected worse performance than that of the ground-only POD. Therefore, the performance of the satellite clock offsets estimated by the Joint 2 strategy indicates the rationality and necessity of adding both clock-free and geometry-free ISL observables.

Although the accuracy of IGSO and GEO satellites using Joint 2 strategy is slightly worse than that of the Joint 1 strategy on the basis of the orbit overlapping differences, which only represents an internal consistency accuracy, the radial component in orbit overlapping differences, the SLR residuals and the characteristic of satellite clock offsets demonstrate that the orbit determination accuracy of the Joint 2 POD strategy is higher than that of Joint 1 POD. In the following sections, the Joint 2 POD strategy is adopted for three types of BDS-3 satellites to analyze whether the number of stations required for orbit determination can be reduced and whether the influence of stations distribution can be reduced or even eliminated after the introduction of ISL.

3.2. Impact of the Number of the Ground Stations

The number of ground tracking stations is crucial for the quality of POD results. On the one hand, sufficient tracking measurements are the guarantee of the accuracy of POD. On the other hand, redundant data will not improve calculation accuracy but will reduce calculation efficiency. As a rule of thumb, with 50 stations, an orbit with satisfactory accuracy can be determined [33]. Under the premise of high POD accuracy, we select different station networks to analyze whether the addition of ISL data can reduce the dependence of POD on the number of ground stations. These tracking station networks consist of different numbers of stations that are evenly distributed globally. With 10 stations as a unit interval, the number of stations is increased from 10 to 50 for the orbit determination. Different networks of tracking stations marked with different colors are plotted in Figure 6.

The overlapping orbit errors of BDS-3 satellites are plotted in Figure 7. The 3D errors of three types of satellites show an expected trend, that is, the more stations used, the higher the obtained accuracy. The accuracy of orbits obtained using 50 stations is similar to the accuracy of the results from 80 stations (see Figure 3). By using more stations, visible improvement is found in the along-track and cross-track directions, while the radial accuracy of these three satellites does not show much change. More importantly, according to the slope of the broken line in Figure 7, it can be inferred that the orbit accuracy changes most when the number of stations decreases from 20 to 10.

Table 4. Percentage improvement when increasing the number of tracking stations.

Types	Directions	Percentage Improvement
		10 to 20	20 to 30	30 to 40	40 to 50
MEO	A	18.1%	−0.8%	3.0%	10.4%
	C	24.0%	10.5%	6.4%	3.3%
	R	4.5%	0.0%	1.0%	4.3%
	3D	19.7%	3.3%	4.4%	7.5%
IGSO C38	A	25.4%	3.6%	−0.2%	2.7%
	C	15.0%	4.6%	−0.3%	1.5%
	R	6.9%	0.0%	4.0%	3.6%
	3D	20.0%	3.9%	0.0%	2.1%
IGSO C39	A	9.7%	3.2%	2.3%	−1.2%
	C	11.7%	2.3%	1.8%	−0.7%
	R	3.3%	2.5%	1.2%	−1.2%
	3D	10.1%	2.7%	2.0%	−1.0%
GEO C59	A	14.5%	3.3%	1.5%	3.4%
	C	6.4%	2.0%	−0.3%	1.3%
	R	3.1%	1.8%	0.8%	2.8%
	3D	11.5%	2.8%	0.9%	2.7%

lists the improvement percentage of orbit accuracy obtained by increasing the number of tracking stations. For MEO satellites, the subsatellite points fall between 55° N and 55° S on the ground, thus increasing the number of global tracking stations can continue to improve the orbit accuracy, especially when the number of stations increases from 10 to 20 (19.7%) and from 40 to 50 (7.5%). With regard to GEO and IGSO satellites, since the subsatellite points are only over the Asia-Pacific region, the addition of globally distributed stations contributes to a slight improvement in orbit accuracy. Starting with 20 stations, each additional 10 stations bring a maximum of 3.9% improvement in the 3D RMS. Therefore, with the supplement of sufficient ISL data, an appropriate ground tracking network should be selected for the three types of satellites according to the required accuracy, which can improve the calculation efficiency of POD and ensure accuracy.

3.3. Improvement on POD Using Regional Stations

At present, the ground facilities of the BDS are all located in China. To verify the ISL improvement on POD when using only regional stations, which is crucial to the generation of broadcast ephemeris, a regional tracking network consisting of seven stations in China and three stations in Australia is selected. The results with or without ISL data are compared to precise orbits obtained from 10 globally distributed stations. The distribution of regional and global tracking stations, as well as the subsatellite points of three types of satellites, are plotted in Figure 8 in which the red dots represent regionally distributed stations and the green dots represent globally distributed stations.

Figure 9 shows the comparison between the orbits estimated with different tracking networks, based on the two-day overlapping differences between two consecutive three-day orbit arcs. Obviously, after adding ISL data, the joint POD accuracy of all three types of BDS-3 satellites is significantly improved regardless of the distribution of stations. With the help of ISL, the orbits of the whole constellation obtained using 10 regionally distributed stations are comparable to the orbits obtained using 10 globally distributed stations with a difference of only 0.9 cm. For MEO satellites, due to the reasonable station distribution, the orbit accuracy based on global ground stations alone is higher than the results of using only regional ground stations. However, this difference resulting from the different distributions of stations is greatly reduced after adding ISL data. The 3D error using regional stations is decreased by 69.7% to 10.7 cm.

In terms of IGSO satellites, subject to the insufficient number of stations located in the Asia-Pacific region in our selected global tracking network, the orbit accuracy using the 10 regional stations (the green bars) in Figure 9b,c shows better performance than that using the 10 global stations (the yellow bars). Especially in the along-track direction, the RMS of IGSO satellites using the regional stations is 41.8% smaller than that using the global stations. After the addition of ISL data to POD, the accuracy obtained by the two strategies is comparable and the 3D error using regional stations is reduced to 15.2 cm.

The orbit of GEO satellite C59, located at longitude 140°, shows similar performance as IGSO satellites due to the global network selected. For the ground-only POD results, the along-track error generated from 10 regional stations is worse than that generated from 10 global stations. This may be attributed to the fact that the selected regional stations are mainly concentrated in two narrow regions: China in the northern hemisphere and Australia in the southern hemisphere. The GEO satellite fixed above the equator is stationary relative to the ground. In this case, this tracking station distribution is not sensitive to the change of the orbits in the along-track direction. The addition of ISL effectively reduces the along-track error of the regional station POD orbit from 87.4 to 19.3 cm, a significant decrease of 78%. The 3D accuracy is therefore improved by about 73.9%. It follows that ISL can significantly improve the accuracy of the orbits computed with observations from regional stations.

4. Conclusions

The evaluation of the BDS-3 joint orbit determination using the ground-tracking and ISL observations is presented by using an improved joint POD strategy, which is referred to as the Joint 2 POD strategy. This strategy makes full use of the derived ISL geometry-free observables and clock-free observables, thus the performance of computed satellite orbits and clock offsets is improved. The experimental results show that the accuracy of the Joint 2 orbits is higher than ground-only orbits, ISL-only orbits, and Joint 1 POD orbits. For the two joint POD strategies with ISL data, the 3D accuracy of MEO satellites estimated using the Joint 2 strategy is improved by about 9%. In terms of IGSO and GEO satellites, the addition of ISL geometric-free observables increases the errors in the cross-track or along-track direction and results in a slight decrease of 8.5% in the 3D accuracy of the IGSO satellites. However, the results of the independent SLR check and the performance of satellite clock offset estimated with orbit parameters demonstrate that the orbit determination accuracy of the Joint 2 POD strategy is higher than that of Joint 1 POD regardless of the types of orbits.

We also prove that the addition of ISL data can reduce the dependence on the ground stations. The results show that the accuracy of the orbit improves as the number of stations increases. When the number of stations reaches 50, the corresponding orbits achieve similar accuracy as the results using 80 stations. However, as the number of stations is larger than 20, each additional 10 stations contributes to a maximum of 7.5%, 3.9%, and 2.8% improvement for MEO, IGSO, and GEO satellites 3D accuracy, respectively.

Further, the Joint 2 orbits estimated from only 10 regional stations or 10 global stations are analyzed to verify ISL’s enhancement to POD. With the use of ISL measurements, the orbit accuracy of all three types of BDS satellites is greatly improved regardless of the distribution of stations, and the orbits using the global stations are slightly better than those of regional stations with a difference of only 0.9 cm in 3D errors. For the orbits estimated using 10 regional stations, the 3D error of MEO satellites is decreased by 69.7%. In the case of IGSO satellites, the 3D error using regional stations is reduced to 15.2 cm. For GEO satellites, the addition of ISL effectively reduces the along-track error of regional station POD orbits by 78%.

Based on the above analysis, we conclude that the impact of ISL for MEO satellites is rather minor when there are enough ground tracking stations. As for GEO satellites, it is obvious that any improvement of geometry resulting from ISL ranging will have a positive impact on POD. Therefore, an appropriate ground tracking network should be selected for the three types of BDS-3 satellites according to the required accuracy and calculation efficiency.