A BDGIM-Based Phase-Smoothed Pseudorange Algorithm for BDS-3 High-Precision Time Transfer

Abstract

Single point positioning (SPP) can meet the requirements of the majority of real-time time transfer applications. Meanwhile, a single-frequency (SF) receiver is cheaper than a dual-frequency receiver. However, SPP performance can be greatly affected by large pseudorange observation noise. Phase smoothing the pseudorange is an effective approach to reduce pseudorange noise. Since the classical phase-smoothed pseudorange algorithm does not account for the effect of ionosphere delay, we propose a BDGIM-based phase-smoothed pseudorange algorithm to eliminate the ionospheric delay and apply it to BeiDou Navigation Satellite System (BDS-3) SPP time transfer. In this paper, we first evaluate the performance of the BeiDou global ionospheric delay correction model (BDGIM) and compare it with that of the BeiDou Klobuchar model to determine if it is practical to incorporate the BDGIM into our suggested method. The performance of the BDGIM is better than that of the Klobuchar model. The mean RMS value of the BDGIM is 2.6 Total Electron Content Unit (TECU). The average ionospheric correction rate of the BDGIM is 75.5%. Then, we investigate the performance of the improved SF SPP time transfer. The performance of the improved SPP time transfer is much better than that of the traditional SPP time transfer. Compared with the traditional time transfer, the average Type A uncertainty of the improved time transfer is 2.08 ns, which is reduced by about 11.1% from the time transfer without it. Regarding frequency stability, the modified Allan deviations of the improved time transfer are 1.43E-12 and 1.68E-13 at 960 s and 61,440 s, with improvements of 51.2% and 59.9%, respectively.

1. Introduction

The law of motion and time scales are both characterized by the fundamental physical quantity known as time. High-precision time transfer techniques have been widely used in the fields of communication, power, financial services, autonomous driving, and aerospace, and they are a prerequisite for maintaining synchronization between time-keeping laboratories and time scales in different countries [1,2].

With the development of Global Navigation Satellite System (GNSS) technology, GNSS-based time transfer techniques have been enhanced. Currently, GNSS-based time transfer methods with pseudorange observations can be included, such as common-view (CV) [3], all-in-view (AV) [4,5], and two-way satellite time and frequency transfer (TWSTFT) [2] time transfer methods, which have greatly improved the accuracy of time transfer [3,6,7]. The aforementioned time transfer techniques have been the subject of extensive scholarly study. In 1980, Allan et al. [1] summarized the GPS time transfer methods and gave a detailed description of the error correction methods involved in time transfer. Defraigne et al. [8] presented error-processing strategies and the Common GPS GLONASS Time Transfer Standard related to CV time transfer. Wei et al. [9] developed a CV time transfer method using GEO satellites, and it had a time transfer performance comparable to that of TWSTFT. In addition, Zhang et al. [10] analyzed the effect of ionospheric delay on CV time transfer. Skakun et al. [11] considered the effect of phase ambiguity in CV time transfer, which greatly improved time transfer performance. AV time transfer experiments were carried out by Harmegnies et al. [12] utilizing GPS and GLONASS measurements, and the results demonstrated that AV time transfer performed better than CV time transfer. To overcome the drawbacks of the TWSTFT and PPP techniques, Wang et al. [13] combined both for time transfer research. Moreover, Zhang et al. [14] investigated the contribution of the Doppler shift to TWSTFT, and the experimental results showed that the effect of the Doppler shift on time transfer is so small that it can be ignored.

As is well-known, CV time transfer can determine the difference between the local atomic clock time and the GNSS satellite clock time of two remote stations using the GNSS clock time as a reference. The difference between the two atomic clock times is determined by comparing the two atomic clocks. Within a certain geographical distance, the satellite clock error, orbit error, and atmospheric delay can be eliminated, and the accuracy of remote time transfer improves, which has the benefits of high flexibility, all-weather service, and high accuracy. However, it has obvious disadvantages, as CV time transfer is limited by the geographical distance and the accuracy of pseudorange observations. As the distance between two stations increases, fewer identical satellites can be tracked, and the accuracy of time transfer decreases. With the release of IGS precise products, AV time transfer was formally employed for time comparison between international time laboratories in 2006 [5]. AV time transfer, in contrast to CV time transfer, is not limited by distance; the time difference between two locations may be estimated directly; and the accuracy of AV time transfer is comparable to that of CV time transfer. Unfortunately, however, AV time transfer is relatively dependent on ephemeris. AV time transfer can be negatively affected by satellite orbit errors, clock errors, and atmospheric delays near both stations. Similarly, the accuracy of AV time transfer is also limited by pseudorange observations. In addition to CV and AV time transfer, TWSTFT time transfer, which utilizes the characteristics of satellite signals transmitting and receiving in the same path and opposite direction to improve time transfer performance by eliminating the effects of station position error, satellite error, and atmospheric delay, is also an important means for the current international atomic time TAI calculation. Its time transfer accuracy is one order of magnitude higher than that of CV time transfer [15]. The quality of pseudorange observations has a significant impact on time transfer, as shown by the analysis above; hence, it is imperative to improve their quality. The accuracy of phase observations is about 100 times higher than that of pseudorange observations, and phase observations are less susceptible to the multipath effect. However, the carrier phase is affected by integer ambiguity. Pseudorange observations are the opposite of carrier phase observations. Therefore, smoothing the pseudorange using a high-precision carrier phase is an effective method to improve pseudorange accuracy. References [16,17,18] verified the feasibility of dual- and tri-frequency phase-smoothed pseudoranges. For the vast majority of civil users, pseudorange single point positioning (SPP) is sufficient to meet the needs of various application scenarios in daily life. For multi-frequency users, ionospheric delay errors may be removed by combining multi-frequency signals. However, for single-frequency (SF) users, ionospheric delay variation is the main factor contributing to the deterioration of phase-smoothed pseudoranges, so the key to improving the performance of SF phase-smoothed pseudoranges is to accurately calculate and deduct the ionospheric delay variation between adjacent epochs. Therefore, a phase-smoothed pseudorange method considering ionospheric delay variation is proposed in this paper, and it is applied to BDS-3 SF SPP time transfer, which belongs to AV time transfer. The ionospheric delay variation is calculated using the BeiDou global broadcast ionospheric delay correction model (BDGIM).

This paper is structured as follows: In Section 2, we provide a brief description of the BDGIM model, a detailed derivation of a phase-smoothed pseudorange model based on the BDGIM, and a presentation of the accompanying SF SPP time transfer method. In Section 3, the datasets used for SF SPP time transfer and its processing methods are given. In Section 4, we first evaluate the BDGIM’s performance compared with the Klobuchar model. Then, the performance of SF SPP time transfer is assessed. Finally, conclusions are summarized in Section 5.

2. BDGIM-Based Phase-Smoothed Pseudorange Algorithm

2.1. BDGIM

With a revised BDGIM that has 9 broadcast ionospheric parameters and 17 non-broadcast ionospheric parameters, BDS-3 was finally finished in 2020 after the project started in 2017 [19,20,21]. The ionospheric non-broadcast parameters are solidified to the receiver terminals, while the broadcast parameters are communicated via navigation messages. Twenty-four sets of coefficients corresponding to 24 h are provided at once by the BDGIM, which is updated once every day. The ionospheric delay can be obtained by using the following formula:

$$I=M_F·\frac{40.28\times {10}^{16}}{f^2}·[A_0+{\displaystyle \sum _{i=1}^9{\alpha }_i{A}_i}]$$

(1)

where $M_F$ is the mapping function, which can project the ionospheric delay from the station’s zenith direction to that of the satellites and stations in their line-of-sight direction; f denotes the frequency of the signal; $A_0$ is the predicted ionospheric delay determined from the ionospheric pierce point (IPP) and observation time; ${\alpha }_i(i=1~9)$ denotes the ionospheric broadcast parameter; and $A_i(i=1~9)$ denotes the Legendre function value associated with the geomagnetic latitude and longitude. The specific calculation of the symbols in Equation (1) is shown below.

The mapping function $M_F$ can be expressed as follows:

$$M_F=\frac{1}{\sqrt{1-{\left(\frac{R}{R+H}·\cos (Ele)\right)}^2}}$$

(2)

where H is the height of the ionosphere, equal to 400 Km; R is the radius of the Earth, equal to 6,378,137 m; and Ele is the satellite elevation angle.

$$A_0={\displaystyle \sum _{j=1}^{17}{\beta }_j}{B}_j$$

(3)

$$\{\begin{array}{l}{\beta }_j=a_{0,j}+{\displaystyle \sum _{k=1}^{12}(a_{k,j}·\cos {\omega }_k{t}_p+b_{k,j}·\sin {\omega }_k{t}_k)}\\ {\omega }_k=\frac{2\pi }{T_k}\end{array}$$

(4)

where $a_{0,j},a_{k,j}$, and $b_{k,j}$ represent the non-broadcast parameters of the BDGIM; $t_p$ stands for the odd hour of the day; and $T_k$ denotes the forecast period.

$$B_j=\{\begin{array}{c}{N}_{nj,mj}^{}·P_{nj,mj}^{}(\sin {\phi }^{\prime })·\cos (m_j·{\lambda }^{\prime })\hspace{1em}{m}_j\ge 0\\ N_{nj,mj}^{}·P_{nj,mj}^{}(\sin {\phi }^{\prime })·\sin (-m_j·{\lambda }^{\prime })\hspace{1em}{m}_j<0\end{array}$$

(5)

$$\{\begin{array}{rrrrrrrrrr}{n}_i=& 0& 1& 1& 1& 2& 2& 2& 2& 2\\ m_i=& 0& 0& 1& -1& 0& 1& -1& 2& -2\end{array}$$

(6)

where $N_{nj,mj}^{}$ represents the normalized function; $P_{nj,mj}^{}$ represents the standard Legendre function; and ${\phi }^{\prime }$ and ${\lambda }^{\prime }$ represent the geomagnetic latitude and longitude of the IPP in the sun-fixed coordinate system, respectively. $A_i$ is calculated in the same way as $B_j$.

Then, the vertical total electron content (VTEC) in the station’s zenith direction can be expressed as

$$VTEC=A_0+{\displaystyle \sum _{i=1}^9{\alpha }_i{A}_i}$$

(7)

In addition, the following restrictions are imposed in order to avoid negative VTECs calculated by the BDGIM:

$$VTE{C}_{\max }=\{\begin{array}{ll}\max ({\alpha }_0/10,VTEC)& {\alpha }_0\ge 35\\ \max ({\alpha }_0/8,VTEC)& 20\le {\alpha }_0<35\\ \max ({\alpha }_0/6,VTEC)& 12\le {\alpha }_0<20\\ \max ({\alpha }_0/4,VTEC)& else\end{array}$$

(8)

2.2. Classical Phase-Smoothed Pseudorange Principle

Hatch proposed a method of pseudorange smoothing using the phase difference between epochs in 1982 [22], and it essentially combines a high-precision phase and low-precision pseudorange observations. The method not only weakens the pseudorange noise and multipath effect but also avoids fixing the ambiguity, which allows the receiver observations to be utilized to their full potential. That is,

$${\tilde{P}}_k^{}={\omega }_k{P}_k^{}+(1-{\omega }_k)[{\tilde{P}}_{k-1}^{}+\lambda ({\varphi }_k^{}-{\varphi }_{k-1}^{})]$$

(9)

where k and k−1 denote the kth and k−1th epochs, respectively; ${\omega }_k$ is the smoothing weight factor; ${\tilde{P}}_k^{}$ and $P_k^{}$ represent the smoothed and the original pseudorange observations of the kth epoch; and ${\varphi }_k^{}$ and ${\varphi }_{k-1}^{}$ denote the phase observations of the kth and k−1th epochs. Since the phase observation noise is substantially less than the pseudorange observation noise, the accuracy of the smoothed pseudorange observation is greatly improved. However, it should be noted that Equation (9) ignores the effect of the ionospheric delay error variation in the kth and k−1th epochs, resulting in a less-than-optimal smoothing effect. The classical phase-smoothed pseudorange method usually treats the ionospheric delay variation as zero, in which case, the pseudorange accuracy is not effectively improved when the ionospheric delay variation is more drastic. Therefore, considering the variation in the ionospheric delay between epochs, we re-derived the principle of phase-smoothed pseudorange.

2.3. Improved Phase-Smoothed Pseudorange Principle

In general, the observation equations of the raw pseudorange P and raw phase L are expressed as follows [23,24]:

$$P=\rho +c(d{t}_r-d{t}^s)+T+I+c(d_r-d_{}^s)+{\epsilon }_P$$

(10)

$$L=\lambda \varphi =\rho +c(d{t}_r-d{t}^s)+T-I+\lambda (N+b_r-b_{}^s)+{\xi }_L$$

(11)

where the symbols r and s denote the receiver and satellite, respectively; $\rho$ denotes the distance from the satellite to the receiver (m); $d{t}_r$ and $d{t}^s$ represent the receiver and satellite clock errors (s), respectively; T denotes the tropospheric wet delay (m); I is the ionospheric delay; $d_r$ and $b_r^{}$ are the pseudorange and carrier phase hardware delay at the receiver end, respectively; $d_{}^s$ and $b_{}^s$ denote the pseudorange and carrier phase hardware delay at the satellite end, respectively; $\lambda$ is the wavelength on frequency $f$ (m); N is the integer ambiguity (cycle); and ${\epsilon }_P$ and ${\xi }_L$ indicate the unmodeled noises in meters.

We can see from Equation (11) that the carrier observations have integer ambiguity, which must be addressed before pseudorange smoothing can be applied. The effect of the ambiguity can be minimized by the carrier variation between the consecutive epochs, and high-precision carrier phase variation can also be obtained at the same time, supposing that cycle slip does not occur while the satellite is being tracked continuously. The following equations can be found by differentiating the pseudorange and carrier phase measurements of the consecutive epochs, respectively:

$$\Delta P_{k,k-1}^{}=P_k^{}-P_{k-1}^{}={\rho }_k^{}-{\rho }_{k-1}^{}+T_k^{}-T_{k-1}^{}+I_k^{}-I_{k-1}^{}+{\epsilon }_{P,k}^{}-{\epsilon }_{P,k-1}^{}$$

(12)

$$\Delta L_{k,k-1}^{}=L_k^{}-L_{k-1}^{}={\rho }_k^{}-{\rho }_{k-1}^{}+T_k^{}-T_{k-1}^{}-I_k^{}+I_{k-1}^{}+{\xi }_{L,k}^{}-{\xi }_{L,k-1}^{}$$

(13)

Note that Equations (12) and (13) have an ionospheric delay variation $I_k-I_{k-1}$ between the kth and k−1th epochs. This can be obtained using the BDGIM, and it can be expressed as $\Delta I_{BDGIM}=I_k-I_{k-1}$.

Finally, the improved phase-smoothed pseudorange algorithm can be written as follows:

$$\{\begin{array}{l}{\tilde{P}}_k^{}={\omega }_k{P}_k^{}+(1-{\omega }_k)[{\tilde{P}}_{k-1}^{}+\lambda ({\varphi }_k^{}-{\varphi }_{k-1}^{})+2\Delta I_{BDGIM}]\\ {\tilde{P}}_1^{}=P_1\end{array}$$

(14)

with ${\omega }_k=1/k$ [25].

Assume that the pseudorange observation error ${\epsilon }_P$ and the phase observation error ${\xi }_L$ are independent of each other and that their variances are ${\sigma }_P^2$ and ${\sigma }_L^2$, respectively. According to the error propagation law, we have ${\sigma }_{\tilde{P}}^2=\frac{1}{k}{\sigma }_P^2+\frac{k-1}{k}{\sigma }_L^2$. Since the carrier phase accuracy is much greater than the pseudorange accuracy, ${\sigma }_P^2\gg {\sigma }_L^2$, we have

$${\sigma }_{\tilde{P}}^2\approx \frac{{\sigma }_P^2}{k}$$

(15)

From the above equation, we can see that the accuracy of the smoothed pseudorange is approximately $\sqrt{k}$ times better than that of the original pseudorange. Therefore, we can infer that time transfer using the phase-smoothed pseudorange performs better. In this work, a weighted least squares algorithm is used to predict the receiver clock offset.

The principle of phase-smoothed pseudorange SF SPP time transfer based on the BDGIM is displayed in Figure 1. First, the BDS-3 observations and ephemeris products of the two stations are obtained. Note that the broadcast parameters of the BDGIM are included in the broadcast ephemeris. Then, we use the phase messages and ionospheric delay to smooth the pseudorange observations and to obtain the smoothed pseudorange observations. In addition, the receiver clock offset, i.e., the difference of local time ($T_1$ and $T_2$) relative to the GNSS reference time (Ref), is estimated in real time using SPP. Finally, the difference between the two receiver clock offsets is calculated.

3. Dataset and Processing Strategies for SF SPP Time Transfer

In this work, we collect observations from six MGEX stations, all of which can receive the signals of BDS-3. The dataset covers a five-day period from 27 March to 31 March 2022 (DOY 86–90, 2022). The broadcast ephemeris is provided by the China Satellite Navigation Office (CSNO) (ftp://ftp2.csno-tarc.cn/cnav/, accessed on 10 September 2022), and the GIM ionosphere products are downloaded from CODE. Figure 2 displays the distribution of the six stations. The information of the six stations is listed in Table 1. Note that the USUD is taken as the central node for time transfer, and five time links, namely, BRUX-USUD, LCK3-USUD, CUSV-USUD, STR1-USUD, and HOB2-USUD, are designed. Table 2 summarizes the processing strategies of SF SPP. The BDS-3 B1I signal is adopted to implement SF PPP time transfer with five consecutive days of ephemeris data.

4. Results and Discussion

For BDS single-frequency time transfer, the ionospheric delay is an error that cannot be ignored, and it directly determines the performance of time transfer. One of the main methods used to remove the ionospheric delay is the broadcast ionospheric model. Furthermore, the release of the BDS-3 BDGIM means that it will supersede the BeiDou Klobuchar, so we first evaluate the BDGIM before conducting time transfer experiments. The performance of the BDGIM is evaluated with reference to GIM products. By making a comparison with the Klobuchar model, the ionospheric delays over different sites are analyzed. Then, the performance of phase-smoothed pseudorange SF SPP time transfer based on the BDGIM is evaluated.

4.1. Evaluation of the BDGIM

The GIM product is used as a reference standard because it is updated every hour and has an inner precision of about 2~8 TECU Total Electron Content Unit (TECU), high accuracy, and continuity. GIM is the most representative global ionospheric VTEC grid product, with stable accuracy and reliability, currently released internationally, and it has high reliability for evaluating the broadcast ionosphere model [29,30]. The VTEC over the tracking station can be obtained by interpolation with the GIM product. The VTEC difference (dVTEC) between the broadcast ephemeris and the GIM is calculated. The root mean square (RMS) error and correction rate (CR) are employed to evaluate the ionospheric model, and they can be written as

$$RMS=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N{(VTE{C}_{model}^k-VTE{C}_{GIM}^k)}^2}}$$

(16)

$$CR=\frac{1}{N}{\displaystyle \sum _{k=1}^N(1-\frac{|VTE{C}_{model}^k-VTE{C}_{GIM}^k|}{VTE{C}_{GIM}^k})}$$

(17)

where N denotes the number of epochs involved in the calculation; $VTE{C}_{model}^k$ is the VTEC calculated with the broadcast ionosphere model; and $VTE{C}_{GIM}^k$ is the VTEC value calculated from the GIM product.

Figure 3 shows the time sequence of the dVTEC values of the BDGIM and Klobuchar models relative to the GIM product, which reflects the trend of the accuracy of the different models over time. We have three findings here. First, the BDGIM has the smallest error, and its dVTEC series has the smallest fluctuation range, fluctuating around the zero value, with a peak-to-peak dVTEC of about 10 TECU. The dVTEC series of the Klobuchar model has a relatively larger fluctuation range, with a peak-to-peak dVTEC of approximately 20 TECU. Second, the broadcast ionospheric error varies among the stations. For the BDGIM, the ionospheric errors are basically the same for each station. For the Klobuchar model, the ionospheric model errors vary greatly from station to station. Among all the stations, the ionospheric error of the BRUX station varies very drastically, which may be caused by insufficient BDS-3 observation data [19], while the ionospheric error of the STR1 station varies relatively steadily. This further indicates that the BDGIM outperforms the Klobuchar model. Third, the dVTEC values of both the BDGIM and Klobuchar models are characterized by a daily cyclic variation, which is because the strong sunlight intensity during the day increases the electron content, and the weak sunlight intensity at night results in lower electron content.

For better analysis and comparison, we calculated the RMS and CR values of the BDGIM and Klobuchar models with the GIM product as a benchmark, as displayed in Figure 4 and Figure 5, respectively. The statistical results indicate that the RMSs of the BDGIM are smaller than those of the Klobuchar model. The RMSs of the BDGIM are all within 4 TECU, with a mean value of 2.6 TECU, while the RMS values of the Klobuchar model are between 4 and 5 TECU, with a mean value of 4.3 TECU. For the ionospheric correction rate, the BDGIM has a better ionospheric correction rate. The CR values of the BDGIM are greater than 60% for all stations, with an average CR of 75.5%. The CRs of the Klobuchar model range from 55% to 75%, with a mean value of 68.2%. Interestingly, we found that the CR values for both the BDGIM and Klobuchar models are greater in the northern hemisphere than those in the southern hemisphere [31]. In summary, the BDGIM outperforms the Klobuchar model. Therefore, the BDGIM was used to conduct time transfer experiments in order to verify our proposed method.

4.2. Performance Evaluation of Phase-Smoothed Pseudorange SPP Time Transfer Based on BDGIM

In this subsection, we focus on verifying the feasibility of our proposed time transfer method. There are more than 120 monitoring stations distributed around the world, all of which contribute to the calculation of the GBM product. In addition, the GBM product is calculated using the B1I/B3I IF observations of a three-day arc, which have a high accuracy of about 2.0 cm for orbit and about 75.0 ps for clock offset. Moreover, the accuracy of precise point positioning (PPP) time transfer can reach a sub-nanosecond level [32]. Therefore, the time transfer results calculated from BDS-3 PPP using the GBM product are taken as a reference. Since dual-frequency PPP has better time transfer performance than SF SPP [33], the B1I/B3I ionosphere-free observations are used to carry out GBM PPP.

To better assess SF SPP time transfer performance, the following experiments are conducted with USUD and LCK3 stations as examples. Before performing time transfer, the number of BDS-3 visible satellites (NSAT) and the Position Dilution of Precision (PDOP) of the two stations are counted, as displayed in Figure 6. The average NSAT value of the USUD station is 6.9, and the mean PDOP is 4.7. For LCK3, the average NSAT value is 11.7, and the mean PDOP is 1.6.

Now turn to Figure 7 and Figure 8, which present the receiver clock errors of the USUD and LCK3 stations relative to the GBM PPP solutions. Note that scheme 1 indicates the traditional SF SPP, and scheme 2 indicates the improved SF SPP. From the two figures, three findings can be found. First, it is easy to see that the clock error sequence of scheme 1 has many discrete points, while the clock error sequence of scheme 2 is quite smooth, and the discrete points are significantly reduced, almost becoming a thin line, which can be observed in the local enlargement in Figure 7. This shows that the smoothing algorithm has considerable potential in time transfer, and the pseudorange observations after smoothing are closer to the real value, reducing the noise of the pseudorange observations. The receiver clock errors are more stable after smoothing, which is more beneficial to improving time transfer performance. Second, we note that there is a convergence period of the clock error during the initial period; this is because the PPP algorithm uses a Kalman filter to estimate the receiver clock offset, which takes some time to converge to the “true value”. However, there is no convergence in SPP, which uses the least squares method. When the receiver clock offsets solved by the SPP and PPP are differenced, there is a convergence problem in the clock errors, i.e., the left part of the red vertical line in Figure 7, which is removed to avoid the impact of the convergence problem on the performance evaluation of SPP time transfer. In this paper, the clock errors are considered to have converged when the clock error difference between 20 consecutive adjacent epochs is kept within 0.1 ns. Third, the clock error of the two stations does not fluctuate around the zero value. This is due to the various receiver hardware delays at various frequencies [34]. The receiver clock offset calculated by the GBM PPP absorbs the receiver hardware delay of the B1I/B3I signal, while that of the SF SPP absorbs the hardware delay of the B1I signal.

To represent the performance of the phase-smoothed pseudorange more comprehensively, we present the pseudorange residuals for both the USUD and LCK3 stations, as shown in Figure 9. Each color represents each satellite. The fluctuation range of the pseudorange residuals is significantly reduced after smoothing. The RMSs of the pseudorange residuals for the USUD station in scheme 1 and scheme 2 are 1.26 m and 0.95 m, respectively, reducing by 24.6%. The RMSs of the pseudorange residuals for the LCK3 station in scheme 1 and scheme 2 are 1.18 m and 0.97 m, respectively, reducing by 17.8%, which further indicates that the improved method enables better SPP performance.

As described above, the smaller clock errors indicate better time transfer results. We give the clock difference of SF SPP time transfer for the LCK3-USUD time link, as shown in Figure 10. It is not difficult to see that scheme 2 has a certain smoothing effect on scheme 1. Similarly, the convergence problem in the initial period is caused by the PPP.

Due to limited space, only the LCK3 and USUD stations are analyzed in detail above. Next, we count the STD of the receiver clock errors after convergence for each station and the improvement of scheme 2 compared with that of scheme 1 in STD, as shown in Figure 11. Two findings are concluded. First, for all stations, the STDs of scheme 2 are smaller than those of scheme 1. The STDs of scheme 1 range from 1.16 ns to 2.11 ns, with an average STD of 1.75 ns, while the STDs of scheme 2 range from 1.08 ns to 1.81 ns, with an average STD of 1.51 ns. The average improvement of the STD values of scheme 2 is 13.6% with respect to those of scheme 1. This further validates the feasibility of our proposed method. Second, the STDs of the clock errors of the mid- and low-latitude stations are smaller than those of the high-latitude ones. On the one hand, this may be due to the relatively poor quality of the BDS-3 observations. On the other hand, this may be due to the ionospheric correction rate being lower at high latitudes than at low latitudes [19], which can also be seen in Figure 5.

Taking the time transfer results with the GBM products as a reference, we count the clock difference box display of five time links, as shown in Figure 12. The clock differences of scheme 1 and scheme 2 are concentrated between ±5 ns, but the box length of scheme 2 is shorter than that of scheme 1, which suggests that scheme 2 has a smaller clock error fluctuation range and a higher time transfer accuracy than scheme 1. It is not difficult to see that the outliers in scheme 2 are significantly less than those in scheme 1, which indicates that, after smoothing, a large number of outliers are pulled back to normal values. Similarly, the peak-to-peak value of the clock difference is also reduced after smoothing.

To further evaluate the contribution of our proposed method to SF SPP time transfer, the STDs of the SF SPP solutions after convergence are calculated to evaluate the Type A uncertainty of the time links. Additionally, the improvement of time transfer in scheme 2 relative to that in scheme 1 in STD is also displayed, as shown in Figure 13. The STD values of scheme 1 range from 2.05 ns to 2.56 ns, with a mean STD of 2.34 ns, while the STD values of scheme 2 range from 1.92 ns to 2.40 ns, with a mean STD of 2.08 ns. The improvement of the STD values of scheme 2 ranges from 5.0% to 19.8% with respect to those of scheme 1, with an average value of 11.1%. In general, the BDGIM-based phase-smoothed pseudorange method provides a good improvement to the SF SPP time transfer performance in the time domain.

We analyze the performance of SF SPP time transfer in the time domain. Going one step further, we analyze the frequency stability of time transfer, which is another indicator for evaluating time transfer performance. We use the modified Allan deviation to indicate the frequency stability of the time link, which can be expressed as [35]

$$MDEV=\frac{1}{2n_{}^2{\tau }_0^2(N-3n+1)}\sum _{i=1}^{N-3n+1}[\sum _{j=i}^{i+n-1}(x_{j+2n}^{}-2x_{j+n}^{}+x_j^{}){]}^2$$

(18)

where N represents the number of sampling points; n denotes the smoothing factor, whose maximum value is less than half of N; ${\tau }_0^{}$ stands for the sampling interval; and $x_j$ denotes the clock difference.

Figure 14 shows the MDEV of the five time links. The frequency stability of all time links is greatly improved after smoothing. The average frequency stabilities of scheme 1 and scheme 2 are 2.93 × 10⁻¹² and 1.43 × 10⁻¹² at 960 s, respectively, with an improvement of 51.2%. The average frequency stabilities of scheme 1 and scheme 2 are 4.19 × 10⁻¹³ and 1.68 × 10⁻¹³ at 61,440 s, respectively, with an improvement of 59.9%. Therefore, the BDGIM-based phase-smoothed pseudorange method also provides a good improvement to the SF SPP time transfer performance in the frequency domain.

5. Conclusions

SPP has been widely used in the field of civil time and frequency. However, SPP relies solely on pseudorange-only observations, which is not conducive to the improvement of SPP accuracy because pseudorange observation noise is rather high. The accuracy of pseudoranges can be increased by smoothing them with a high-precision phase. However, the classical phase-smoothed pseudorange method ignores the effect of ionosphere delay. In this study, we develop a BDGIM-based phase-smoothed pseudorange algorithm, and we investigate the contribution of the method to BDS-3 SF time transfer. First, to verify the feasibility of the BDGIM, we evaluated the accuracy of the BDGIM by comparing it with the Klobuchar model using the GIM product as the benchmark. Then, we evaluated the improved SPP time transfer performance using the observation data from six stations over five consecutive days. The following conclusions can be drawn from the experimental results:

First, when comparing the broadcast ionosphere models, the VTEC error of the BDGIM fluctuates in a narrower range, with a peak-to-peak value within 10 TECU, while the peak-to-peak value of the VTEC error of the Klobuchar model is approximately 20 TECU. The RMS value of the BDGIM is 2.6 TECU, which is 1.7 TECU less than that of the Klobuchar model. The average CR value of the BDGIM is 75.5%, which is an improvement of 7.3%. Thus, the BDGIM outperforms the Klobuchar model.

Second, the receiver clock error sequence after smoothing is relatively more concentrated, and its STD value is decreased by 13.6%, from 1.75 ns to 1.51 ns. Regarding time transfer, Type A uncertainty with the improved SPP time transfer is reduced from 2.34 ns to 2.08 ns, with a reduction of 11.1%. Regarding frequency stability, the frequency stability of the improved SPP time transfer is greatly improved. The MDEVs at 960 s and 61,440 s are 4.19 × 10⁻¹³ and 1.68 × 10⁻¹³, with an increase of 51.2% and 59.9%, respectively.

Finally, we are aware that broadcast ephemeris and pseudorange observations are both used in SPP time transfer. Although the smoothed pseudorange helps to some extent with SPP time transfer accuracy, the quality of the satellite orbits and clock offsets provided by broadcast ephemeris is poor, which can potentially limit SPP time transfer performance. In order to further improve the capability of SPP time transfer, we will concentrate our efforts in the future on increasing the quality of broadcast ephemeris orbit and clock offsets.