**1. Introduction**

The BeiDou Navigation Satellite System (BDS) provides positioning, navigation, and timing (PNT) information to global users. The system was developed in three phases. The first is the demonstration system (BDS-1), developed in 2003 and mainly providing services through two first-generation experimental geostationary Earth orbit (GEO) satellites. The second phase is a regional system (BDS-2) for the Asia–Pacific region, operational since 25 October 2012. This phase comprises five GEO, five inclined geosynchronous orbit (IGSO), and four medium-altitude Earth orbit (MEO) satellites. The third phase is the global system (BDS-3), operational since July 2020, which comprises 30 satellites, including 3 GEO satellites, 3 IGSO satellites, and 24 MEO satellites [1,2]. The BDS-2 is expected to remain in service for at least another decade [3,4], although BDS-3 is already fully operational. In combination, the system provides exceptional potential for PNT users as it employs more BDS satellite resources.

Multi-GNSS constellations have known benefits of time and frequency transfer with respect to precision, integrity, and availability, because of the increased number of available satellites [5–10]. In particular, the multi-GNSS carrier phase technique (CP) has been

**Citation:** Zhang, P.; Tu, R.; Tao, L.; Wang, B.; Gao, Y.; Lu, X. Preliminary Analysis of Intersystem Biases in BDS-2/BDS-3 Precise Time and Frequency Transfer. *Remote Sens.* **2022**, *14*, 4594. https://doi.org/ 10.3390/rs14184594

Academic Editors: Shengfeng Gu, Xiaopeng Gong, Yidong Lou and Chuang Shi

Received: 19 August 2022 Accepted: 9 September 2022 Published: 14 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

proven to outperform a single GNSS with respect to remote time and frequency transfer with precision at the nanosecond level [11,12]. Accordingly, this technique is applied widely by international time laboratories for the campaign of time and frequency transfer. However, multi-GNSS precise time and frequency transfer is not free of challenges. In this regard, the unification of coordinate frames and time benchmarks and the differences in receiver hardware delays related to using signals from different systems should be considered [13–15]. Such biases are called intersystem biases (ISBs) and they affect data processing when combined multi-GNSS are employed for time and frequency transfer [16].

Abundant satellite sources are provided by BDS-3 and BDS-2, with common coordinate frames and time benchmarks, namely Geodetic Coordinate System 2000 (CGCS2000) and BDS Time (BDT). The combination of BDS-3 and BDS-2 signals has aroused the interest of numerous scholars. For instance, Li et al. analyzed the ISB between the BDS-2 and BDS-3 experimental systems using the International GNSS Monitoring and Assessment System (iGMAS) and the Multi-GNSS Experiment (MGEX) observations [17]. These authors point out that there a systematic bias between BDS3 and BDS2, and the systematic biases are different for stations in different networks. Pan et al. evaluated the multi-GNSS positioning performance with a priori ISB constraint. According to these authors, the ISBs in BDS-2/BDS-3 could not be ignored [18]. Zhao et al. showed that the ISB values of stations with the same type of receiver were similar, while a substantial difference existed for different receiver types [19]. Fu et al. introduced the ISB parameter between BDS-2 and BDS-3 and improved the standard deviation (STD) of all satellite clocks, even using overlapping B1I/B3I measurements [20]. Although the performance of BDS-2/BDS-3 time and frequency transfer has been investigated thoroughly [21–25], few studies have focused on the characteristics and mechanism of ISBs in BDS-2/BDS-3.

It is well-known that the GNSS CP technique relies heavily on precise satellite products for time and frequency transfer. However, it is not clear whether ISB characteristics also depend on satellite clock products. Accordingly, estimating the ISBs in BDS-2/BDS-3 and determining their influence on the performance would facilitate optimal utilization of the BDS system.

In this study, we discuss the temporal characteristics, parameter composition, generation mechanism, and the effect of ISBs in BDS-2/BDS-3 on time and frequency transfer. This study starts with a brief description of mathematical models of BDS-2/BDS-3 time and frequency transfer with the CP technique, after which ISB estimation methods are discussed. Subsequently, the data sources with several types of receivers, antennas, external time and frequency references, and data-processing strategies are introduced. In addition, the temporal characteristics of ISB based on satellite products from different analysis centers and the reasons for changes are analyzed. Furthermore, the influence of different ISB estimation methods on the time and frequency transfer results is evaluated. Finally, we present our summary and conclusions.

#### **2. Materials and Methods of Multi-GNSS CP Technique**

The classical CP time transfer comprises two components for the pseudorange with time information and carrier phase measurements with millimeter-scale noise levels. The ionosphere-free (IF) linear combination with dual-frequency observation is used to eliminate the effect of the first-order ionosphere. The observation model can be expressed as follows, taking the single BDS as an example [12]:

$$\begin{cases} \begin{aligned} P^i\_{IF} &= \rho^i\_r + c \cdot \left( dt^i\_r - dt^{s,i} \right) + T\_{trop} + b^i\_{r,IF} - b^{s,i}\_{IF} + c^i\_{IF} \\\ L^i\_{IF} &= \rho^i\_r + c \cdot \left( dt^i\_r - dt^{s,i} \right) + T\_{trop} + \lambda\_{IF} \cdot \left( N^i\_{IF} + B^i\_{r,IF} - B^{s,i}\_{IF} \right) \\\ &\quad + \varepsilon^i\_{IF} \end{aligned} \end{cases} \tag{1}$$

where the superscript *i* denotes the BDS system (i.e., BDS-2, BDS-3). *P<sup>i</sup> IF* and *<sup>L</sup><sup>i</sup> IF* are the IF combination of pseudorange and carrier phase observation, respectively; *ρ<sup>i</sup> <sup>r</sup>* is the geometric distance between the phase center of the satellite and receiver antenna; *c* is the speed of light in vacuum; *dt<sup>i</sup> <sup>r</sup>* and *dts*,*<sup>i</sup>* represent the receiver and satellite clock offsets, respectively. *Ttrop* is the tropospheric delay; *b<sup>i</sup> <sup>r</sup>*,*IF* and *<sup>B</sup><sup>i</sup> <sup>r</sup>*,*IF* are the IF combination of reciever pseudorange and phase hardware delay, respectively; *bs*,*<sup>i</sup> IF* and *<sup>B</sup>s*,*<sup>i</sup> IF* are the IF combination of satellite pseudorange and phase hardware delay. *λIF* is the wavelength of the IF combination; *Ni IF* is the phase ambiguity of the IF combination; and *<sup>e</sup><sup>i</sup> IF* and *<sup>ε</sup><sup>i</sup> IF* are measurement noise for the pseudorange and carrier phase observation, respectively. The GNSS satellite and receiver phase center offset and variation, phase wind-up, solid tide, ocean load, pole tide, and relativistic delay should also be considered, although these terms are not listed in Equation (1).

In precise time and frequency transfer employing the CP technique, the satellite orbit and clock products provided by the International GNSS Service (IGS) are used to reduce the orbit and clock errors. The IF combination is used in the data processing of the IGS to determine the satellite orbit and clock parameter. The satellite hardware delay *bs*,*<sup>i</sup> IF* is absorbed in the satellite clock offset, providing a reference for the receiver clock offset. Therefore, the receiver IF combination pseudorange hardware delay *b<sup>i</sup> <sup>r</sup>*,*IF* is assimilated into the receiver clock offset. The CP delays *B<sup>i</sup> <sup>r</sup>*,*IF* and *<sup>B</sup>s*,*<sup>i</sup> IF* at the receiver and satellite are related closely to the ambiguity parameter *N<sup>i</sup> IF* and lumped with the estimated ambiguity parameter. Therefore, Equation (1) can be written further as:

$$\begin{cases} \overline{P}\_{IF}^i = \rho\_r^i + c \cdot \overline{dt}\_r^i + T\_{trop} + e\_{IF}^i\\ \overline{L}\_{IF}^i = \rho\_r^i + c \cdot \overline{dt}\_r^i + T\_{trop} + \lambda\_{IF} \cdot \overline{N}\_{IF}^i + \varepsilon\_{IF}^i \end{cases} \tag{2}$$

where *P<sup>i</sup> IF* and *L i IF* are the actual pseudorange and CP observations when using the IGS precise satellite orbit and clock products. Therefore,

$$
\overline{P}\_{IF}^i = P\_{IF}^i + c \cdot dt^{s\_r i} + b\_{IF}^{s,i} \tag{3}
$$

$$
\overline{L}\_{IF}^i = L\_{IF}^i + c \cdot dt^{s\_r i} + b\_{IF}^{s\_i i} \tag{4}
$$

$$\overline{dt}\_r^i = dt\_r^i + b\_{r,IF}^i \tag{5}$$

$$\overline{N}\_{IF}^i = N\_{IF}^i + B\_{r,IF}^i - \frac{b\_{r,IF}^i}{\lambda\_{IF}} - B\_{IF}^{s,i} + \frac{b\_{IF}^{s,i}}{\lambda\_{IF}} \tag{6}$$

where *dt<sup>i</sup> <sup>r</sup>* and *<sup>N</sup><sup>i</sup> IF* are the new receiver clock offset and ambiguity parameter, lumped with the corresponding hardware delays.

When BDS-2 and BDS-3 are combined for precise time and frequency transfer, an additional ISB parameter (*ISBC*3,*C*2) between the two systems is introduced to obtain a common receiver clock offset that references a unique system time scale [26–28]. Therefore, the BDS-2/BDS-3 time and frequency transfer model can be written as:

$$\begin{cases} \overline{P}\_{IF}^{\complement2} = \rho\_r^{\complement2} + c \cdot \overline{dt}\_r^{\complement2} + T\_{\text{trop}} + \varepsilon\_{IF}^{\complement2} \\ \overline{L}\_{IF}^{\complement2} = \rho\_r^{\complement2} + c \cdot \overline{dt}\_r^{\complement2} + T\_{\text{trop}} + \lambda\_{IF} \cdot \overline{N}\_{IF}^{\complement2} + \varepsilon\_{IF}^{\complement2} \\ \overline{P}\_{IF}^{\complement3} = \rho\_r^{\complement3} + c \cdot \overline{dt}\_r^{\complement2} + ISB^{\complement3,\complement2} + T\_{\text{trop}} + \epsilon\_{IF}^{\complement3} \\ \overline{L}\_{IF}^{\complement3} = \rho\_r^{\complement3} + c \cdot \overline{dt}\_r^{\complement2} + ISB^{\complement3,\complement2} + T\_{\text{trop}} + \lambda\_{IF} \cdot \overline{N}\_{IF}^{\complement3} + \epsilon\_{IF}^{\complement3} \end{cases} (7)$$

where superscript *C*2 and *C*3 denote the BDS-2 and BDS-3 system. Among the parameters, the unique receiver clock offset *dtC*<sup>2</sup> *<sup>r</sup>* is the most interesting parameter for precise time transfer, which is determined jointly by the BDS-2 and BDS-3 observations, although it is denoted simply as the BDS-2 system. Two stations, A and B, located at different places on

Earth, are equipped with their corresponding time and frequency references. The operation of time transfer between the two references can be obtained using the following expression:

$$\begin{aligned} \Delta T\_{A,B} &= \overline{dl}\_r^{\mathbb{C}2}(A) - \overline{dl}\_r^{\mathbb{C}2}(B) \\ &= \begin{pmatrix} t\_A^{\mathbb{C}2} - BDT + b\_{r,IF}^{\mathbb{C}2}(A) \end{pmatrix} - \begin{pmatrix} t\_B^{\mathbb{C}2} - BDT + b\_{r,IF}^{\mathbb{C}2}(B) \end{pmatrix} \\ &= t\_A^{\mathbb{C}2} - t\_B^{\mathbb{C}2} + b\_{r,IF}^{\mathbb{C}2}(A) - b\_{r,IF}^{\mathbb{C}2}(B) = \Delta t\_{A,B}^{\mathbb{C}2} + \Delta b\_{r,IF}^{\mathbb{C}2}(AB) \end{aligned} \tag{8}$$

where *t <sup>C</sup>*<sup>2</sup> is the external time and frequency reference when BDS-2 observation is used. The term BDT is the BeiDou time scale, which uses the international system of units (SI) second without leap seconds, connects with universal time coordinated (UTC) through UTC (NTSC, national time service center), and the deviation of BDT to UTC is maintained within 50 nanoseconds. The initial epoch of BDT is 00:00:00 on January 1, 2006, of UTC (BeiDou Navigation Satellite System Open Service Performance Standard, Version 3.0, May 2021). The Δ*t C*2 *<sup>A</sup>*,*<sup>B</sup>* is the clock difference between two external time and frequency references at stations A and B; Δ*bC*<sup>2</sup> *<sup>r</sup>*,*IF*(*AB*) is the delay difference of receiver pseudorange between stations A and B, usually calibrated using the time-transfer link calibration or receiver calibration approaches [29,30]. After the combined observation equation (Equation (7)) has been transformed and linearized, the unknown parameter vector X can be summarized as:

$$\mathbf{X} = \left[ \mathbf{x}, y, z, \overline{d\mathbf{t}}\_r^{\mathbf{C2}}, \mathbf{I} \mathbf{S} \mathbf{B}^{\mathbf{C3}, \mathbf{C2}}, \mathbf{T}\_{\text{trop}\prime} \overline{\mathbf{N}}\_{\text{IF}\prime}^{\mathbf{C2}} \overline{\mathbf{N}}\_{\text{IF}}^{\mathbf{C3}} \right] \tag{9}$$

where (x, y, z) is the station coordinate parameter.

In order to further clarify the origin of ISB, referring to Equations (5), (7), and (8), the ISB can be written further as:

$$\begin{split} ISB^{\rm C3,C2} &= c \cdot \overline{dt}\_r^{\rm C3} - c \cdot \overline{dt}\_r^{\rm C2} = c \cdot dt\_r^{\rm C3} - c \cdot dt\_r^{\rm C2} + b\_{r,IF}^{\rm C3} - b\_{r,IF}^{\rm C2} \\ &= c \cdot \Delta dt\_r^{\rm C3,C2} + \Delta b\_{r,IF}^{\rm C3,C2} \end{split} \tag{10}$$

The *ISB* theoretically comprises two components, namely the time difference of two receiver clock offsets with different GNSS observations, and the difference in the receiver hardware delays for two GNSS systems [31]. For the former, Δ*dtC*3,*C*<sup>2</sup> *<sup>r</sup>* is a function of the receiver clock offset, which is the difference between the external time and frequency reference and GNSST (GNSS time, GNSST), as discussed previously. Unlike the combination of different GNSSs, such as BDS, GPS, Galileo, and GLONASS, with their individual system time scales, BDT (BeiDou system time, BDT), GPST (GPS system scale, GPST), GST (Galileo time scale, GST), and UTC (SU) include the component of time deviation between different GNSSTs for the term. If Δ*dtC*3,*C*<sup>2</sup> *<sup>r</sup>* does not contain this term, the formula can be expressed as follows:

$$
\Delta dt\_r^{\text{C3,C2}} = \left(t^{\text{C3}} - BDT\right) - \left(t^{\text{C2}} - BDT\right) = t^{\text{C3}} - t^{\text{C2}} \tag{11}
$$

where *t <sup>C</sup>*<sup>3</sup> is the external time and frequency references when using BDS-3 observation and is equivalent to *t <sup>C</sup>*<sup>2</sup> when using a multimode BDS receiver. Considering the occurrence of unknown errors and unmodeled deviation in ISB estimation, we introduced parameter *τISB* to represent these errors. Equation (10) can be written further as:

$$ISB^{C3,C2} = \Delta b^{C3,C2}\_{r,IF} + \tau\_{ISB} \tag{12}$$

As shown by Equations (7)–(9) and (12), the ISB parameter is important when determining the receiver clock and further carrying out the time and frequency transfer.

#### **3. ISB Stochastic Models in Multi-GNSS Time Transfer**

With respect to the ISB parameter, it usually performs three stochastic models, namely white noise, random constant, and random walk process. Although the previous research shows that the stochastic model is closely related to the ISB performance in the multi-GNSS positioning and time transfer [32,33], the performance in BDS-2/BDS-3 is still unclear.

Regarding the white noise process, the ISB is assumed to be uncorrelated between the different epochs. The white noise process is applied widely when the characteristics of a parameter are not known. The model is expressed as

$$Q\_{ISB}(k) \sim N\left(0, \sigma^2\right) \tag{13}$$

where *Q* denotes a covariance; *k* is the epoch index.

For the random constant process, it is estimated as a piecewise mode. As for the entire data block, it is usually divided into sub-blocks; the mathematical model in the sub-blocks is expressed as

$$Q\_{ISB}(k+1) = Q\_{ISB}(k)\tag{14}$$

where *k* is the epoch index.

The random walk process can be formulated as follows:

$$Q\_{ISB}(k+1) = Q\_{ISB}(k) + \omega\_{ISB\prime} \qquad \omega\_{ISB} \sim N\left(0, \sigma\_{\omega\_{ISB}}^2\right) \tag{15}$$

where *ω* is the process noise of a random walk.

## **4. Results**

To explicitly investigate the temporal characteristics of ISBs in BDS-2/BDS-3 time transfer, we collected data from the Multi-GNSS Experiment (MGEX), which is piloted by the IGS to collect all available GNSS observations of new signals since 2013. The MGEX network has expanded to more than 300 stations, providing an excellent opportunity to track multi-GNSS constellations and to conduct tracking data analysis. More than 200 stations are tracking BDS satellites. However, most stations only track the dual-frequency BDS-2 signals, but single-frequency BDS-3 observations or the data received are only for a short period during a day. Moreover, the external time and frequency of atomic clocks are not equipped in most stations. Consequently, a limited number of available BDS-2 and BDS-3 stations are available for analyzing ISB variation. Ten stations that track common dualfrequency (B1I and B3I) signals for BDS-2 and BDS-2, equipped with atomic clocks, and that have a relatively complete receiver period in a day, are collected. The geographical distribution of collected stations is shown in Figure 1. The experiment was conducted from day of year (DOY) 100–120, 2021. Detailed information on these stations, e.g., type of receiver, antenna, and frequency standard, is presented in Table 1.

**Table 1.** Information of utilized BDS stations in the experiment.


**Figure 1.** Geographical distribution of collected stations in the experiment.

During the data processing, the precise time-transfer solution (PTTSol) was used [34]. The BDS satellite orbit and clock products from three IGS analysis centers, WUM, GFZ, and CODE, were employed for further research of the stability of ISBs when combining BDS-2 and BDS-3 for time transfer. The Wuhan University analysis center has provided BDS-3 satellite orbit and clock products since GPS week 2034, with 15 min and 5 min updates, using "wum" ID. GeoForschungsZentrum Potsdam has provided them since GPS week 2081, with 5 min and 30 s intervals, using "gbm" ID. Fortunately, starting from GPS week 2148, the CODE satellite solution has included BDS-3 (apart from GEO satellites), with a 5 min orbit and 30 s clock, using "com" ID. The pseudorange and CP measurements for B1, B3 of BDS-2, and BDS-3 dual-frequency observations are used to alleviate the ionosphere effect. In preprocessing, both the geometry-free combination and the Melbourne–Wübbena combination were used for cycle slip detection. Tropospheric wet delay is typically modeled as the sum of the Saastamoinen model and a random walk process. The receiver clock offset parameter was estimated as a white stochastic noise process. We also considered the wind-up effects on phase measurements. The data-processing strategies in this study are summarized in Table 2.



#### *4.1. Characterization of ISB over Different Time Periods*

As the research was conducted over numerous days, we randomly selected the result on DOY 100, 110, and 120 of 2021 to analyze ISB variation. Figure 2 shows that the ISBs are stable with the different analysis center products on the three days. The average variation over the experimental period is within 0.08m for ISB\_com and 0.07 m for ISB\_gbm and ISB\_wum. The variations on DOY 100 and 120 are within 0.06 m; however, the range on DOY 110 at 0.09 m is much larger than that for the other days. The difference between the minimum and maximum values is approximately 0.28 m for the ISB\_com solution, 0.26 m for ISB\_wum, and 0.25 m for ISB\_gbm. Notably, the mean values of the results from three analysis centers differ markedly for one daily result. Figure 3 shows the mean value character of the results of different analysis centers of 10 stations. Although obvious systematic bias does exist among the com, gbm, and wum results, it is relatively stable among the ten stations. The systematic bias values between the com and wum results are −0.46 m, −0.50 m, and −0.33 m for DOY 100, 110, and 120, respectively, whereas for com and gbm the values are −1.65 m, −1.75 m, and −0.78 m, respectively, i.e., larger than the former. Systematic bias is caused mainly by the different data-processing strategies of the three analysis centers when determining the BDS-2 and BDS-3 satellite orbit and clock products. Figures 2 and 3 show that the ISB trend is generally stable. The stability of ISB\_gbm and ISB\_wum is slightly superior to that of ISB\_com. Remarkably, the mean value of the former four stations equipped with SEPT POLARX5 receivers is not more stable than that of the middle three stations (SEPT POLARX5TR) or the latter three stations (TRIMBLE NETR9). Relevant details on this aspect are presented in the following sections.

**Figure 2.** Daily variations in ISB between BDS-3 and BDS-2 for 10 stations. The panels from top to bottom show the results of the com, wum, and gbm precise products.

**Figure 3.** Average of ISB between BDS-3 and BDS-2 for 10 stations for the com, wum, and gbm solutions.
