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, Reference Kleusberg and Teunissen1996; Collins et al., Reference Collins, Bisnath, Lahaye and Heroux2010):

(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, Reference Kleusberg and Teunissen1996); $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., Reference Wu, Wu, Hajj, Bertiger and Lichten1992), relativistic effect on propagation (Kleusberg and Teunissen, Reference Kleusberg and Teunissen1996) and the phase centre offsets (PCOs) and phase centre variations (PCVs) of the antennas (Schmid et al., Reference Schmid, Dach, Collilieux, Jäggi, Schmitz and Dilssner2016) 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., Reference Tang, Hu, Zhou and Liu2018). 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., Reference Tang, Hu, Zhou and Liu2018):

(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., Reference Yang, Yang, Li, Zhou and Tang2017; Tang et al., Reference Tang, Hu, Zhou and Liu2018), 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. (Reference Ruan, Jia, Feng, Zhu, Huyan, Li and We2020).

2.3 Joint observation equations for integrated OD and TS

The integrated OD and TS method proposed by Ruan (Reference Ruan2018) 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, Reference De Jonge1998). Consequently, Equation (5) should be rewritten as (Ruan et al., Reference Ruan, Jia, Feng, Zhu, Huyan, Li and We2020):

(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., Reference Ruan, Jia, Feng, Zhu, Huyan, Li and We2020), 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.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, Reference Jiao2014) 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., Reference Ruan, Jia, Wu, Feng and Zhu2014) 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, Reference Springer1999) 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, Reference McCarthy and Petit2004) 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., Reference Beohm, Neill, Tregoning and Schuh2006) 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

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