Integrated orbit determination and time synchronization using inter-satellite observations for BDS-3 satellites with satellite clock drifts estimated simultaneously

Abstract

By simultaneously estimating satellite clock drifts (SCDs) as either constant parameters or piece-wise parameters, we present an improved integrated orbit determination and time synchronization approach for BDS-3 satellites with raw inter-satellite link (ISL) observations. Experiments with L-band data from seven monitoring stations in China and ISL data from eight satellites of the third-generation Beidou Navigation Satellite System (BDS-3) were carried out and the two SCD estimation strategies are validated. It is demonstrated that, with SCDs estimated, the quality of orbits and clock offsets is comparable to those with SCDs corrected using predicted values. The accuracy of the estimated orbits and clocks are up to 0.019 m (radial) and 0.095 ns, respectively, with improvements of 95% and 90%, when compared with the results using the L-band data alone. It is also demonstrated that estimating SCDs time slice by time slice is slightly worse in accuracy but superior in coping with possible frequency jump of satellite clocks.

1. Introduction

On 31 July 2020, the third generation of the Beidou Navigation Satellite System (BDS-3) was declared to have achieved full operational capability. Unlike GPS and Galileo, whose monitoring stations are globally distributed, the monitoring stations of BDS-3 are only located in a small region. inter-satellite link (ISL) is a trump card to help BDS-3 to provide competitive global positioning and timing service (Yang et al., 2018, 2019).

The concept of ISL for navigation satellites was initially proposed in the early 1980s to enhance survivability with auto navigation (AutoNav) capability (Codik, 1985; Ananda et al., 1990) if the operational control segment (OCS) fails. However, the ISL payloads onboard BDS-3 satellites are state of the art (Yang et al., 2017), as they work on Ka-band and enable high-precision ranging and fast beam switching, and the phased array antenna employed can scan a large space range by controlling the phase relation of each antenna's elements. As introduced in detail by Yang et al. (2017) and Tang et al. (2018), each BDS-3 satellite can perform bidirectional ranging and communication with others in an asynchronous or exactly concurrent spatial time division (CSTD) mode according to the preset connecting schedule. The ISL not only enables BDS-3 AutoNav capability but also enhances the orbit determination (OD) and time synchronization (TS) in the OCS.

Yang et al. (2017) first demonstrated OD for BDS-3 satellites with Ka-band ISL data and L-band monitoring data, and similar work was done by Wang et al. (2019), Xie et al. (2019) and Yang et al. (2020). It is demonstrated that the orbit accuracy of BDS-3 Medium Earth-Orbiting (MEO) satellites achieved with data from regional stations and ISL is better than 5 cm in the radial direction (Xie et al., 2020). Orbit determination using ISL data alone has also been verified, e.g. by Tang et al. (2017, 2018) and Ren et al. (2017, 2019). Pan et al. (2018) presented a method of TS for BDS-3 satellites with ISL data and L-band two-way satellite time-frequency transfer (TWSTFT) data, which shows the accuracy of predicted clocks is about 0.5 ns. TS with ISL data alone was demonstrated to be superior to that with L-band data, as they are independent of orbit errors (Xie et al., 2020).

Despite different details of the studies quoted, they share a common characteristic that raw ISL observations are not directly used. Instead, the non-simultaneous one-way ISL measurements are converted into two kinds of independent dual-one-way observations: the derived ephemeris observation (DEO) and the derived clock observation (DCO) which are respectively used for OD and TS separately. Shortcomings of this method include: the resulting satellite orbits and clock offsets lacking self-consistency, being quite tedious and time-consuming to match raw one-way observations for each satellite-pair to compute DEOs and DCOs, and so on (Ruan et al., 2020). Therefore, Ruan (2018) proposed an approach of integrated OD and TS directly using raw one-way ISL data. Orbit determination and TS can be achieved using raw ISL data combined with data from monitoring stations or anchor stations or both (Ruan, 2018). The key idea is to model satellite clocks with polynomials in a series of short time slices, correcting clock variations during each time slice using predicted satellite clock drifts (SCDs).

The SCDs are essential to correct satellite clock variations during each time slice (Ruan et al., 2020), as well as for computing the DEOs and DCOs (Tang et al., 2018). A SCD with an accuracy about 1 × 10⁻¹³ s/s is usually obtainable, i.e., from the broadcast navigation messages. However, in practical, predicted SCDs may not always be available or may suffer from unbearable errors due to unpredictable frequency jump or artificial frequency adjustment. In this study, we present a slightly improved approach of integrated OD and TS using raw ISL data with simultaneous estimation of SCDs.

The remaining part of this paper is organised as follows: the second section introduces joint observation equations for integrated OD and TS; in the third section, experiments are carried out to validate the proposed SCD estimation strategies; the fourth section analyses the adaptability of the two SCD estimation strategies to different time slice lengths and possible frequency jump; and conclusions are provided in the final section.

2. Observation equations

This section presents the observation equations for the L-band ionosphere-free (IF) observations and the original Ka-band ISL pseudoranges, as well as the joint observation equations for integrated OD and TS.

2.1 L-band pseudorange and carrier phase

Dual-frequency IF combinations of pseudoranges and carrier phases from monitoring stations are conventionally used for OD and TS of navigation satellites. The IF pseudorange, $P_r^i$, and phase, $L_r^i$, for the satellite i observed at the station r at the $k$-th epoch ${t_k}(k=1,2, \cdots )$ are written as (Kleusberg and Teunissen, 1996; Collins et al., 2010):

(1)\begin{equation}\begin{aligned} & P_r^i({t_k})=\rho _r^i+ {\delta _{r,k}} - \delta _k^i+ {b_r}+ {b^i}+ T_r^i+ \xi \\ & L_r^i({t_k})=\rho _r^i+ {\delta _{r,k}} - \delta _k^i+ {\eta _r}+ {\eta ^i}+ T_r^i+ a_r^i+ \varepsilon \end{aligned}\end{equation}

where all symbols are in distances; $\rho _r^i=|{{\boldsymbol{R}_r}({t_k}) - {\boldsymbol{R}^i}({t_k} - \tau_r^i)} |$ is the geometric distance, with $\tau _r^i$ the propagation time, and ${\boldsymbol{R}_r}$ and ${\boldsymbol{R}^i}$ the positions of the station and the satellite, respectively; ${\delta _{r,k}}$ and $\delta _k^i$ are the clock offsets of the station and satellite, respectively; ${b_r}$ and ${b^i}$ are, respectively, the group hardware delays (HDs) at the receiver and satellite, while ${\eta _r}$ and ${\eta ^i}$ are the phase HDs considered stable over time; $a_r^i$ is the IF phase bias, superposition of the ambiguities and the non-zero initial phases (Kleusberg and Teunissen, 1996); $T_r^i$ is the troposphere delay; $\xi$ and $\varepsilon$ are the measurement noises for the two kind of observations, including higher-order ionosphere delays. In Equation (1), terms such as phase wind-up (Wu et al., 1992), relativistic effect on propagation (Kleusberg and Teunissen, 1996) and the phase centre offsets (PCOs) and phase centre variations (PCVs) of the antennas (Schmid et al., 2016) are omitted but should be corrected with accurate models or data.

2.2 Ka-band intersatellite pseudoranges

The ISL payloads on BDS-3 satellites can only either receive or transmit signals at any time. Each satellite performs bidirectional ranging and communication with other visible satellites successively according to the preset connect schedule. A 3-s time slot is assigned to the planed link of an individual satellite pair, and each of the two satellites costs half the time slot, i.e. 1.5 s, to receive and transmit signals alternately (Tang et al., 2018). The total time required for every linkable satellite pair to complete a bidirectional link is referred to as a polling period whose length depends on the link topology and the scale of the constellation.

Taking the link between satellites j and i, for example, the forward pseudorange measurement, ${K^{i,j}}$_, observed at time ${t^{i,j}}$ at satellite i can be expressed as (Tang et al., 2018):

(2)\begin{equation} {K^{i,j}}({t^{i,j}})= {\rho ^{i,j}}+ {\delta ^i}({t^{i,j}}) - {\delta ^j}({t^{i,j}} - {\tau ^{i,j}}) + {\Delta ^j}+ {\beta ^i}+ \omega \end{equation}

where ${\rho ^{i,j}}=|{{\boldsymbol{R}^i}({t^{i,j}}) - {\boldsymbol{R}^j}({t^{i,j}} - {\tau^{i,j}})} |$ is the geometric distance, while ${\tau ^{i,j}}$ is the signal propagation time from the transmission to the reception; ${\boldsymbol{R}^i}$ and ${\boldsymbol{R}^j}$ are the positions of the two satellites, respectively; ${\delta ^i}$ and ${\delta ^j}$ are the clock offsets of the two satellites; ${\Delta ^j}$ and ${\beta ^i}$ are the group HDs for the ISL transmitter and receiver, respectively, which are considered stable over time; $\omega$ is the measurement noise. In Equation (2), the omitted terms, such as the relativistic effect and PCOs of phased array antennas, should be properly corrected. Similarly, the backward pseudorange measurement is:

(3)\begin{equation}{K^{j,i}}({t^{j,i}}) = {\rho ^{j,i}}+ {\delta ^j}({t^{j,i}}) - {\delta ^i}({t^{j,i}} - {\tau ^{j,i}}) + {\Delta ^i}+ {\beta ^j}+ \omega\end{equation}

Since the ISL of BDS-3 works on the Ka-band, the pseudoranges between satellites are referred to as Ka-pseudoranges. According to the published connect schedule of the BDS-3 ISL (Yang et al., 2017; Tang et al., 2018), the time tags of Ka-pseudoranges follow the relationships:

(4)\begin{equation}\left\{ \begin{array}{@{}l} \min |{{t^{i,j}} - {t^{j,i}}} |=1.5 s\\ {t^{i,j}} \ne {t^{i,k}} \ne {t^{j,k}} \cdots ({\forall i \ne j \ne k} )\end{array} \right.\end{equation}

Assuming that one satellite performs inter-satellite link with other 10 satellites in a polling period, the observation time difference between the first and the last satellites may be larger than 30 s, so the Ka-pseudoranges corresponding to these satellites should no longer be regarded as simultaneous observations, as is the L-band monitoring data. As shown in Figure 2, the clock drifts for some BDS-3 satellites may exceed ${10^{ - 11}}\,s/s$ in magnitude, and variations of satellite clocks in this time scale may exceed 0.3 ns, which cannot be neglected for OD and TS with the highest accuracy requirements.

The traditional method to process the nonsimultaneous ISL observations for OD or TS is to convert the two (forward and backward) one-way ranging observations for any link in the same polling period to two independent derived observations, i.e. the DEO and the DCO (see Appendix A). And the DEOs and DCOs are used for OD and TS separately. The shortcomings of the traditional method have been discussed in detail by Ruan et al. (2020).

2.3 Joint observation equations for integrated OD and TS

The integrated OD and TS method proposed by Ruan (2018) was used in this paper. To form joint equations for the simultaneous L-band observation data and the non-simultaneous ISL data for a short time slice with length of w centred on ${t_k}$(the sampling epoch of the L-band data), satellite clock offsets are modelled with linear polynomials, and the joint observation equations with L-band and Ka-band observations are written as:

(5)\begin{equation}\begin{aligned} & P_r^i({t_k})=\rho _r^i+ {\delta _{r,k}} - \delta _k^i+ {b_r}+ {b^i}+ T_r^i+ \xi \\ & L_r^i({t_k})=\rho _r^i+ {\delta _{r,k}} - \delta _k^i+ {\eta _r}+ {\eta ^i}+ T_r^i+ a_r^i+ \varepsilon \\ & {K^{i,j}}({t^{i,j}})={\rho ^{i,j}}+ [{\delta_k^i+ \dot{\delta }_k^i \ \cdot \ \Delta t_k^{i,j}} ]- [{\delta_k^j+ \dot{\delta }_k^j \ \cdot \ (\Delta t_k^{i,j} - {\tau^{i,j}})} ]+ {\Delta ^j}+ {\beta ^i}+ \bar{\omega } \end{aligned}\end{equation}

where $\Delta t_k^{i,j}={t^{i,j}} - {t_k}$ is the time difference between the observation time of the Ka-pseudoranges and the sampling epoch ${t_k}$, and $|{\Delta t_k^{i,j}} |\le \frac{w}{2}$; $\dot{\delta }_k^i$ and $\dot{\delta }_k^j$ are SCDs. The term $\bar{\omega }$ includes the errors caused by linear modelling for the satellite clock offsets.

In Equation (5), the clock offsets and the group HDs in the equations of IF- and Ka-pseudoranges are inseparable due to their linear correlativity, as are the phase bias, clock offsets and phase HDs in the IF phase equation, which means these terms are not unbiased estimable (De Jonge, 1998). Consequently, Equation (5) should be rewritten as (Ruan et al., 2020):

(6)\begin{equation}\begin{aligned} & P_r^i({t_k})=\rho _r^i+ {\theta _{r,k}} - \theta _k^i+ T_r^i+ \xi \\ & L_r^i({t_k})=\rho _r^i+ {\theta _{r,k}} - \theta _k^i+ T_r^i+ A_r^i + \varepsilon \\ & {K^{i,j}}({t^{i,j}})={\rho ^{i,j}}+ [{\theta_k^i+ \dot{\theta }_k^i \ \cdot \ \Delta t_k^{i,j}} ]- [{\theta_k^j+ \dot{\theta }_k^j \ \cdot \ (\Delta t_k^{i,j} - {\tau^{i,j}})} ]+ \Delta _b^j+ \beta _b^i+ \bar{\omega } \end{aligned}\end{equation}

where $\theta _k^\ast =\delta _k^\ast{-} {b^\ast }$ is the clock offset parameter for satellite *, combination of the pure clock offset and the L-band group HD at the satellite; ${\theta _{{\ast} ,k}}={\delta _{{\ast} ,k}}+ {b_\ast }$ is the clock offset parameter for the station *, combination of the pure clock offset and L-band group HD at the receiver; $\dot{\theta }_k^\ast =\dot{\delta }_k^\ast$ is the SCD parameter for the satellite *; $A_r^i=a_r^i+ {\eta _r}+ {\eta ^i} - {b_r} - {b^i}$ is the estimable phase bias, including L-band IF group and phase HDs of both the satellite and the receiver; $\Delta _b^\ast ={\Delta ^\ast } - {b^\ast }$ is the transmitting Ka-bias for the satellite *, the difference between the ISL transmitter group HD and the L-band group HD; $\beta _b^\ast ={\beta ^\ast }+ {b^\ast }$ is the receiving Ka-bias for the satellite *, the summation of the ISL receiver group HD and the L-band group HD.

Although the satellite clock variations during each time slice can be corrected with predicted SCDs (Ruan et al., 2020), in this study, both the clock drift and the clock offset parameters will be estimated as unknown parameters. The satellite and station clock offsets will be solved time slice by time slice. The SCDs will be solved time slice by time slice as piece-wise parameters or as a constant parameter during each OD solution arc. With estimation of SCDs, no external SCD information is need for OD and TS processing.

3. Validation

To validate the proposed method and evaluate the performance of OD and TS enhanced with ISL under the condition of regional monitoring stations, experiments are carried out and the results are analysed in this section.

3.1 Data

Pseudorange and carrier phase data on B1I and B3I of 8 BDS-3 satellites (B19-22 and B27-30) from 16 monitoring stations operated by the International GNSS Monitoring and Assessment System (iGMAS) (Jiao, 2014) during 1–12 June 2018 (DOY152-163) were collected. Figure 1 shows the distribution of the iGMAS stations and the seven stations in red are located in the territory of the P. R. China.

Figure 1. Distribution of the iGMAS stations used in this study

3.2 Configure of experiments

The Satellite Positioning and Orbit Determination System (SPODS) software (Ruan et al., 2014) was used to perform OD and TS processing on a daily basis with three-day data started on each day. After pre-processing, the original 30-s sampled L-band data are re-sampled at 5-min intervals and formed into IF-combinations. Time slices with 1-min length are assigned centred on every sampling epoch. The IF-pseudoranges and phases are weighted according to the function, ${\sin ^2}e$, of the elevation e, while the Ka-pseudoranges are equally weighted. The a priori precision for the IF pseudoranges, phases and Ka-pseudoranges are set as 2.0 m, 0.02 m and 0.1 m, respectively.

The coordinates of the iGMAS stations are fixed to known values. The Earth orientation parameters (EOPs) data from the IERS website are used. The ECOM1 (Springer, 1999) model is used to model the solar radiation pressure with the five parameters (${D_0},{Y_0},{B_0},{B_c}$ and ${B_s}$) estimated together with the initial state vector. The Saastamoinen model (McCarthy and Petit, 2004) is used to calculate the zenith troposphere delay for stations, with zenith delay parameter estimated every 2 h and the Global Mapping Function (GMF) (Beohm et al., 2006) mapping the troposphere delay from the zenith to the line of sight. One phase-bias parameter in the IF-phase observations is estimated as a float parameter for each data pass free of cycle slips. Two Ka-biases for each satellite are also estimated, with the receiving Ka-bias of satellite B21 as reference. The clock offsets of each station and satellite in each time slice are estimated with one station clock selected as reference.

For GNSS navigation users, the predicted orbits and clock offsets are more concerned, so the accuracy of predicted orbits and clocks for BDS-3 satellites based on the improved integrated OD and TS method is of interest. Once the solution of orbits and clocks are obtained, predictions are conducted for the future 24 h. For this purpose, the 72-h clock offset series are fitted with quadratic polynomials.

As shown in Figure 2, the magnitude of SCDs (predicted values) of BDS-3 satellites do not exceed $5 \times {10^{ - 11}}\,s/s$, the variation of the satellite clock would be smaller than 3.0 ns during a 1-min time slice. As the time slices are centred on the sampling epochs, the uncertainty of Ka-pseudoranges caused by satellite clock variation will be smaller than 0.45 m. To validate the proposed estimation strategies for SCDs, several experiments of integrated OD and TS will be carried out, with different processing strategies for the SCD parameters, including: (1) the SCDs are ignored, i.e. regarded as zero (denoted by EXP1); (2) the SCDs are corrected with the predicted values from the broadcast navigation data (denoted by EXP2); (3) the SCDs are estimated as constant parameters over the three-day solution arc (denoted by EXP3); (4) the SCDs are estimated time slice by time slice (denoted by EXP4). In these 4 experiments, only the seven stations in red were used in order to simulate the distribution of monitoring stations in the BDS-3 OCS as much as possible. For comparisons, experiments without ISL data using the seven stations in China (denoted by EXP5) and all the 16 stations (denoted by EXP6) were also performed to demonstrate the contribution of the ISL.

Figure 2. Predicted SCDs of BDS-3 satellites from the broadcast navigation data

3.3 Results analysis

This section analyses posterior residuals of the Ka-pseudoranges obtained with different strategies and then discusses the precision of orbits, clocks and Ka-biases successively.

3.3.1 Posterior residuals of the Ka-pseudoranges

Figure 3 shows the posterior Ka-pseudorange residuals of the typical link between satellite B27 and B29 as a function of the time difference, $\Delta t_k^{i,j}={t^{i,j}} - {t_k}$. Only the residuals of the three-day solution arc during DOY 152–154 are presented. It can be seen that, in panel (a), the range of residuals in EXP1 is about from −0.4 to 0.4 m with an obvious trend of linear variation. Moreover, the residuals of the individual directional link are centrally symmetrical in general, and there seems to be negative correlativity between residuals of the forward and backward data. These phenomena are caused by ignoring the SCDs and thus disappeared in other panels. As shown, the residuals in EXP2, EXP3 and EXP4 are almost no more than 0.2 m in magnitude.

Figure 3. Ka-pseudorange residuals between satellites B27 and B19 as a function of the time difference between the observation time and the sampling epoch of each 1-min time slice

Figure 4 shows the root-mean-square (RMS) of the Ka-pseudorange residuals with different strategies for each reception satellite over the whole data arc. It can be seen that the RMSs in EXP1 range from 0.11 to 0.25 m, much larger than those in other experiments, while the remains are close to each other, not exceeding 0.08 m without exception. It is noted that the RMSs of EXP4 are slightly smaller than those of EXP2 and EXP3, as the SCD parameters are estimated on a time slice basis leading to decreased redundancy. The overall RMSs are 0.181, 0.054, 0.054 and 0.053 m for EXP1, EXP2, EXP3 and EXP4, respectively. Compared to EXP1, the RMSs in the other three experiments are significantly decreased by more than 70%, which strongly verifies both the estimation strategies for SCDs in EXP3 and EXP4, as the RMSs are as small as those in EXP2 with SCDs corrected.

Figure 4. RMS of Ka-pseudorange residuals for each reception satellite in different experiments over the whole data arc

3.3.2 Orbits

The overlap orbit differences (OODs) in the radial (R), tracking (T) and normal (N) directions are calculated for each 48-h overlapped arcs. The predicted orbits errors (POEs) are computed in a similar way for each 24-h prediction arcs.

Table 1 shows the overall RMS of OODs and POEs for each experiment. It can be seen that in EXP2, EXP3 and EXP4, the statistics in each component for both OODs and POEs are nearly the same with differences no more than 2 mm. Compared to the results obtained in EXP1, the RMS OODs in R, and 3D position are decreased by at least 9.5% and 19.7%, respectively, while the RMS POEs decreased by at least 5.4% and 15.0%, respectively. This confirms that both the two SCD estimation strategies are feasible.

Table 1. Overall RMSs of OODs and POEs (unit: m)

Compared to the results obtained in EXP5, the RMS OODs in R and 3D position are significantly decreased by at least 95.7% and 93.2%, respectively; the RMS POEs are decreased by at least 95.5% and 96.5%, respectively; even for the EXP1, the RMSs of OODs are decreased by 95.2% and 91.6%, respectively, for R and 3D positions, while the RMSs of POEs decreased by 95.3% and 95.9%, respectively. The accuracies of orbits with ISL data, no matter which strategies is adopted to deal with the SCDs, are also significantly better than those achieved in EXP6.

3.3.3 Clocks

When calculating overlap clock differences (OCDs), the mean bias over all satellites in each time slice is removed from the original differences to avoid possible deviation in the reference clocks (Ge et al., 2012). The predicted clock errors (PCEs) are computed in the same way. Besides, clock-fitting residuals (CFRs) are analysed as another indicator to assess the quality of the estimated clock offsets, although they are also largely related to the behaviour of satellite clocks.

Table 2 lists the overall RMSs of OCDs, CFRs and PCEs for the satellite clock solutions obtained in different experiments. It is noted that differences between the RMSs of EXP4 and EXP3 are not larger than 0.002 ns for the OCDs and CFRs, and not larger than 0.005 ns for the PCEs, while the RMSs in EXP3 and EXP2 are identical, which confirms the feasibility of the proposed approaches. Compared with EXP1, the RMSs of OCDs, CFRs and PCEs in EXP3 are decreased by 10.5%, 16.2% and 7.2%, respectively, while in EXP4, the RMSs are decreased by 8.7%, 15.2% and 6.2%, respectively.

Table 2. Overall RMSs of the OCDs, CFRs and PCEs (unit: ns)

Compared with EXP5, the RMSs of OCDs, CFRs and PCEs in EXP1 are decreased by 89.4%, 75.2% and 69.8%, respectively, while in EXP2, EXP3 and EXP4, the RMSs are decreased by at least 90.3%, 78.9% and 71.1%, respectively. Compared with the results of EXP6, the clock accuracy in EXP3 and EXP4 are improved by 84.8% and 84.5%, respectively, while the accuracy of predicted clocks improved by 41.7% and 41.1%, respectively.

3.3.4 Ka-biases

Figure 5 shows the variations of the estimated Ka-biases obtained in EXP3. It is demonstrated that the Ka-biases are stable over time, with variation smaller than 1 ns for individual satellite. There seems to be systematic variation in the time series of the estimated Ka-biases, which may partly due to the fact that the Ka-biases are not unbiased estimable. The variations of the estimated Ka-biases superimpose the variation of the selected reference Ka-biases, the receiving Ka-bias of B21. Moreover, the Ka-biases are the combinations of the hardware delays of Ka-band ISL signal and L-band navigation signal, and the variations of the estimated Ka-biases are superposition of both of the two hardware delays. We can find the variations are negative correlative for the receiving and transmitting Ka-biases of the same satellite, since the L-band hardware delay appears in them with the opposite sign.

Figure 5. Variations of the satellite-specific Ka-biases obtained in EXP3

The estimated Ka-biases obtained with different strategies are quite consistent with each other. For example, the differences between the mean estimated satellite-specific Ka-biases obtained in EXP3 and those in EXP1, EXP2 and EXP4 are not larger than 0.1 ns, 0.003 ns and 0.03 ns, respectively. An effective indicator to assess the quality of the estimated Ka-biases is repeatability. The standard deviations (SDevs) of the estimated Ka-biases obtained in different experiments are presented in Figure 6. It can be seen that the SDevs for each satellite obtained with SCDs estimated are generally smaller than 0.5 ns, which are comparable to those with SCDs corrected and generally smaller than the results in EXP1. The repeatability of the estimated Ka-biases obtained in EXP1, EXP2, EXP3 and EXP4 is 0.398, 0.367, 0.367 and 0.362 ns, respectively.

Figure 6. Standard deviation (SDev) of the estimated Ka-biases of each satellite obtained in different experiments

4. Adaptability to different time slice lengths and possible frequency jump

In this section, experiments are carried out to compare the adaptability of the two SCD estimation strategies to different time slice lengths and possible frequency jump of satellite clocks.

4.1 Influence of different time-slice lengths

In the experiments in the previous section, the time slice length is 1 min, while the sampling interval is 5 min, which leads to using only about 20% of the ISL data. For a specific sampling interval, increasing the time-slice length means more ISL data would be used; it is interesting to note how time slice lengths affect the accuracy of orbit and clock solutions. The experiments with the two SCD estimation strategies, i.e. EXP3 and EXP4, were carried out again with different time-slice length, e.g. 30 s, 50 s, 60 s, 90 s, 120 s, 180 s, 240 s and 300 s. As there are a total of eight satellites, 30 s is selected as the minimum time-slice length to cover at least an entire polling period. Figure 7 gives the overall RMSs of 3D OODs with different time-slice lengths and SCD estimation strategies, while Figure 8 shows the overall RMSs of OCDs. It can be found that the RMSs of the OODs and OCDs are slightly decreased as the time-slice length increases. However, in EXP3, the differences for 3D OODs with different time-slice length are no more than 0.02 m and no more than 0.025 ns for the OCDs, while in EXP4, the differences are smaller than 0.019 m and 0.028 ns, respectively. For different time-slice lengths, the results of orbits and clocks in EXP3 are generally slightly better than those in EXP4. This is probably due to the fact that the number of SCD parameters in EXP4 is 864 times more than that in EXP3, which greatly decreases the redundancy. However, the differences of the RMSs for OOD and OCD are no more than 0.01 m and 0.008 ns, respectively.

Figure 7. Overall RMSs of 3D OODs with different time slice lengths

Figure 8. Overall RMSs of OCDs with different time slice lengths

4.2 Influence of frequency jump

In practice, during the operation of the system, satellite clocks may experience frequency (or phase) adjustment to keep the clock offsets in sync with the system time within a certain range. A more extreme case may also occur is switch between the primary clock and the backup clock. In these cases, predicted SCDs is usually not available, and estimating the SCDs time slice by time slice is expected to have advantages. Unfortunately, in the short period covered by the experimental data, such evens did not really happen. To compare the ability of the two strategies to cope with such situations, a frequency jump for each satellite was simulated at the middle of each three-day solution arc, and EXP3 and EXP4 were reconducted, with their results denoted by EXP3′ and EXP4′, respectively. The magnitudes of the simulated frequency jumps are twice that of the predicted clock drifts, as shown in Figure 2, ranging from $1.5 \times {10^{ - 12}}$ to $4.5 \times {10^{ - 11}}$ $s/s$.

Table 3 shows the overall RMSs of ISL residuals, OODs and OCDs achieved in EXP3′ and EXP4′. Compared to the results in EXP3′, the RMSs of OODs in EXP4′ are decreased by 9.5%, 12.5%, 8.3% and 9.7% for R, T and N components and 3D positions, while the RMS of OCDs decreased by 29.2%; the RMS of the ISL residuals decreased by 70.5%. It can be found that the results achieved in EXP4′ are nearly identical to those in EXP4, while the results of EXP3′ are significantly worse than those of EXP3. Therefore, in the case of possible frequency jump of satellite clocks, estimating the SCDs on a time-slice basis is preferred.

Table 3. Overall RMSs of OODs and OCDs with different SCD estimation strategies in case of frequency jump

5. Discussion

Satellite clock frequency stability is an important factor for the presented approach. The baseline of the presented method has proven to be highly flexible and robust because adapts to different strategies for handling the SCDs. As demonstrated, the results in EXP2, EXP3, and EXP4 are nearly the same; even when the SCDs are roughly ignored, orbits and clocks with considerable accuracy are also achievable. This should be attributed to the fact that the actual SCDs are fairly small and the idea to divide the solution arc into a series of short time slices. The latter makes the errors caused by the (estimated) SCDs keep in a tolerable range. Furthermore, the midpoint of each time slice is selected as the reference time of the SCDs. Therefore, in each time slice, errors cause by the (estimated) SCDs for the data before and after the midpoint are opposite in sign and their influence on the results is cancelled out to a large extent.

For the most precise OD and TS processing, the broadcast SCDs from the navigation data are usually accurate enough. However, estimating the clock drifts would be preferred for it does not rely on any external SCD information. Normally, it is appropriate to estimate one SCD parameter for each satellite for the whole solution arc, but estimating SCDs on a time-slice basis may have the advantage of dealing with unstable satellite clock behaviour and even possible satellite clock frequency jumps (or adjustments).

6. Conclusions

This paper proposes an improved integrated OD and TS approach using the asynchronous ISL data via simultaneously estimation of the satellite clock drifts (SCDs) as either constant parameters or piece-wise parameters. With the improved method, no additional SCD information is need. Using the L-band monitoring data from seven iGMAS stations in China and the ISL data from eight BDS-3 satellites, OD and TS experiments with different strategies to deal with the SCDs were carried out. The precision of the orbits and clocks is evaluated with overlapping comparison. It is demonstrated that the quality of orbits, clocks and Ka-biases achieved with the two SCD estimation strategies are comparable to those with SCDs corrected. The precisions of the estimated orbits and clocks are about 0.019 m (radial) and 0.095 ns, respectively. The improvements are 95.7% and 90.3% when compared with the results using the L-band data alone. The accuracies of predicted orbits and clocks are 0.035 m (radial) and 0.498 ns, respectively, while the improvements are more than 95% and 71%. We also demonstrated that the accuracies of orbits and clocks with regional stations and ISL data are even significantly greater than those with global stations. Both the two SCD estimation strategies show good adaptability to different time-slice lengths and estimating SCDs as piece-wise parameters is superior in coping with possible clock frequency jumps.