**1. Introduction**

With the development and application of Global Navigation Satellite Systems (GNSSs), GNSSs are now gradually moving from dual-frequency to multi-frequency, modes and likewise, multi-frequency precise point positioning (PPP) is widely studied by many scholars [1–4]. Among them, the MEO and ISGO satellites of the BeiDou-3 (BDS-3) can provide data in five frequencies (B1C, B1I, B2a, B3I, and B2b), GPS Block IIF and Block III satellites can provide triple-frequency observation data, and the Galileo system currently has 26 satellites providing five-frequency data [5–7]. Information on the available multifrequency GPS, Galileo, and BDS-3 satellites is provided in Table 1. It is well-known that accurate satellite orbits and clocks are important prerequisites for PPP. Zhou et al. [8]. performed PPP analysis using the orbit and clock products of iGMAS and obtained GNSS kinematic PPPs of 1.4, 1.2, and 2.9 cm in the E, N, and U directions, respectively, along with orbit/clock agreement of 1.5 cm and 60 s, respectively, compared to the orbit/clock of the

**Citation:** Chen, Y.; Mi, J.; Gu, S.; Li, B.; Li, H.; Yang, L.; Pang, Y. GPS, BDS-3, and Galileo Inter-Frequency Clock Bias Deviation Time-Varying Characteristics and Positioning Performance Analysis. *Remote Sens.* **2022**, *14*, 3991. https://doi.org/ 10.3390/rs14163991

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

Received: 13 July 2022 Accepted: 10 August 2022 Published: 16 August 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/).

IGS. Yang et al. [9] used a triple-frequency ambiguity solution based on undifferentiated observations for satellite clock estimation and compared it with ambiguity floating-point clock and found that the ambiguity-fixed clock was solution improved by 32% and 42.9% in the horizontal and vertical directions. However, with the widespread use of triplefrequency observations, the impact of periodic variation in satellite phase hardware delay on the triple-frequency data is becoming significant [10,11]. Montenbruck et al. [12] found that the carrier phase observations of L1, L2, and L5 of the GPS have an inconsistency of 20 cm, which was labeled as the inter-frequency clock bias (IFCB).

**Table 1.** Information on the available multi-frequency GPS, Galileo, and BDS-3 satellites.


In order to make better use of the multi-frequency observations, numerous scholars have investigated the IFCB. Montenbruck et al. [13] employed prior correction of the satellite IFCB to weaken the impact of IFCB on PPP. Pan et al. [14] found that the east, north, and up accuracy improved from 3.1 cm, 1.1 cm, and 3.3 cm to 2.1 cm, 0.7 cm, and 2.3 cm, respectively, after taking into account the triple-frequency PPP with IFCB compared to the PPP localization with uncorrected IFCB. In another work, Li et al. [15] observed that correcting the IFCB while performing the triple-frequency uncalibrated phase delay (UPD) estimation can significantly improve the quality of extra-wide-line UPD. Fan et al. [16] analyzed the IFCB of the GPS Block IIF satellite via eight months of observations, and concluded that the inter-peak amplitude could reach up to 10–40 cm and that the IFCB varied more during the eclipse than during the other periods. Furthermore, Zhao and Montenbruck et al. [17,18] analyzed the IFCB of BDS-2 satellites, which are affected by 2–4 cm of IFCB. In contrast, Steigenberger and Zhao [19,20] showed good consistency among the triple frequencies of Galileo and QZSS. In addition to these works, a better understanding of the characteristics of the BDS-3 IFCB is needed. Furthermore, the compatibility of ionosphere-free combination with non-combination estimated IFCB was also focused on in one of the reported studies [21].

Additionally, some scholars further analyzed the cycle variation in IFCB. In this regard, Gong [22] pointed out that due to solar illumination variations, such as the relative Sun– satellite–Earth geometry changes, the internal temperature of the satellite also changes, leading to periodic changes in the satellite phase delay, thereby resulting in the periodic changes in IFCB. Therefore, Li et al. [23] modeled the estimated IFCB using linear and fourth-order harmonic functions and reported more than 89% correction of the IFCB, with an average fitted RMS of 1.35 cm for the GPS IFCB. Moreover, Zhang [24] found that the periodic variation in IFCB is to some extent related to the orbital plane in which the satellite is located, and for the two satellites distributed in same orbital plane, IFCB shows similar amplitudes and waveforms.

Considering the current status of the existing research on IFCB, this paper utilized 117 MGXE (muti-GNSS experiment) stations worldwide for IFCB estimation, and analyzed the intra-day and inter-day time-varying characteristics of IFCB. Furthermore, this work analyzed the impact of IFCB on the GPS, BDS-3, and Galileo multi-frequency precise point positioning in terms of IFCB amplitude, PPP positioning accuracy, and post-test residuals, and investigated the short-term stability of IFCB.

#### **2. Methods**

*2.1. Triple-Frequency Uncombined PPP Model with IFCB*

The pseudorange *P<sup>s</sup> <sup>r</sup>*,*<sup>i</sup>* and carrier phase *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*<sup>i</sup>* observation equations of GNSSs are given as [15]:

$$\begin{array}{l}P\_{r,i}^{s} = \rho\_r^s + dt\_r - dt^s + T\_r^s + \mu\_i I\_{r,1}^s + d\_{r,i} + d\_i^s + \varepsilon\_{r,i}^s\\L\_{r,i}^s = \rho\_r^s + dt\_r - dt^s + T\_r^s - \mu\_i I\_{r,1}^s + \lambda\_i N\_{r,i}^s + b\_{r,i} + b\_i^s + \zeta\_{r,i}^s\end{array} \tag{1}$$

where *s* is the satellite and *r* is the receiver; *i* (*i* = 1, 2, 3) is the carrier frequency; *ρ<sup>s</sup> <sup>r</sup>* is the geometric distance between satellite and receiver; *dtr* and *dt<sup>s</sup>* are the receiver and satellite clock errors, respectively; *T<sup>s</sup> <sup>r</sup>* is the tropospheric delay; *I<sup>s</sup> <sup>r</sup>*,1 is the slant ionospheric delay at frequency *<sup>f</sup>*1; *<sup>μ</sup><sup>i</sup>* <sup>=</sup> *<sup>f</sup>* <sup>2</sup> 1 *f* 2 *<sup>i</sup>* is the frequency-dependent ionospheric delay amplification factor; *f* is the frequency; *λ<sup>i</sup>* is the carrier wavelength; *N<sup>s</sup> <sup>r</sup>*,*<sup>i</sup>* is the carrier phase ambiguity; *dr*,*<sup>i</sup>* and *d<sup>s</sup> <sup>i</sup>* are the code hardware delays at the receiver and satellite, respectively; *br*,*<sup>i</sup>* and *bs <sup>i</sup>* are the phase hardware delays from the receiver and satellite, respectively; and *<sup>ε</sup><sup>s</sup> r*,*i* and *ζ<sup>s</sup> <sup>r</sup>*,*<sup>i</sup>* are the unmodeled error and the observation noise of the code and carrier phase observations for each frequency.

The phase hardware deviation has obvious time-varying characteristics, and accordingly, it can be decomposed into a constant part and a time-varying part [7,13], as elaborated in Equation (2).

$$\begin{cases} \begin{array}{c} b\_{r,i} = \overline{b}\_{r,i} + \delta b\_{r,i} \\ b\_i^s = \overline{b}\_i^s + \delta b\_i^s \end{array} \end{cases} \tag{2}$$

where *br*,*<sup>i</sup>* and *b s <sup>i</sup>* are the constant parts of the receiver and satellite phase hardware delays, respectively, while *δbr*,*<sup>i</sup>* and *δb<sup>s</sup> <sup>i</sup>* are the time-varying parts of the receiver and satellite phase hardware delays. Moreover, the following variables are defined herein for the ease of expression:

$$\begin{cases} \begin{aligned} a\_{12} &= \frac{f\_1^2}{f\_1^2 - f\_2^2} \\ \beta\_{12} &= -\frac{f\_2^2}{f\_1^2 - f\_2^2} \\ DB\_{12}^s &= d\_1^s - d\_2^s \\ DCB\_{r,12} &= d\_{r,1} - d\_{r,2} \\ \delta DPB\_{12}^s &= \delta b\_1^s - \delta b\_2^s \\ \delta DB\_{r,12} &= \delta b\_{r,1} - \delta b\_{r,2} \\ \delta b\_{1F12}^s &= a\_{12}\delta b\_1^s + \beta\_{12}\delta b\_2^s \\ \delta b\_{r,IF12} &= a\_{12}\delta b\_{r,1} + \beta\_{12}\delta b\_{r,2} \\ d\_{IF12}^s &= a\_{12}d\_1^s + \beta\_{12}d\_2^s \\ d\_{r,IF12} &= a\_{12}d\_{r,1} + \beta\_{12}d\_{r,2} \end{aligned} \end{cases} \tag{3}$$

where *α*<sup>12</sup> and *β*<sup>12</sup> are the frequency factors' ionosphere-free combinations; *DCB<sup>s</sup>* <sup>12</sup> and *DCBr*,12 are the satellite and receiver differential code bias values, respectively; *δDPB<sup>s</sup>* 12 and *δDPBr*,12 are the satellite and receiver time-variant parts of differential phase bias, respectively; *δb<sup>s</sup> IF*<sup>12</sup> and *δbr*,*IF*<sup>12</sup> are the ionosphere-free combination time-variant parts of receiver and satellite phase hardware delays, respectively; and *d<sup>s</sup> IF*<sup>12</sup> and *dr*,*IF*<sup>12</sup> are the IF pseudorange hardware delays at the receiver and satellite, respectively.

After applying the precise satellite clock, track, and DCB product corrections, a triplefrequency uncombined PPP model with IFCB is expressed as:

$$\begin{cases} P\_{r,1}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r + \overline{T}\_{r,1}^s + \delta b\_{r,1}^s + \varepsilon\_{r,1}^s \\ P\_{r,2}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r + \mu\_2 \overline{T}\_{r,1} + \delta b\_{r,2}^s + \varepsilon\_{r,2}^s \\ P\_{r,3}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r + \mu\_3 \overline{T}\_{r,1}^s + IFB\_r + \delta b\_{r,3}^s + \varepsilon\_{r,3}^s \\ L\_{r,1}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r - \overline{t}\_{r,1}^s + \lambda\_1 \overline{N}\_{r,1}^s + \zeta\_{r,1}^s \\ L\_{r,2}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r - \mu\_2 \overline{T}\_{r,1}^s + \lambda\_2 \overline{N}\_{r,2}^s + \zeta\_{r,2}^s \\ L\_{r,3}^s = u\_r^s x + d\overline{t}\_r + m\_r^s Z\_r - \mu\_3 \overline{T}\_{r,1}^s + \lambda\_3 \overline{N}\_{r,3}^s + IFCB + \zeta\_{r,3}^s \end{cases} (4)$$

with

$$\begin{cases} \begin{aligned} \frac{d\overline{I}\_{r}}{dt} &= dt\_{r} + d\_{r,IF12} + \delta b\_{r,IF12} \\ \overline{I}\_{r,1}^{s} &= I\_{r,1}^{s} + \beta\_{12} DCB\_{r,12} - \beta\_{12} (\delta DPB\_{12}^{s} + \delta DPB\_{r,12}) \\ IFB\_{r} &= d\_{r,3} - d\_{r,IF12} - \mu\_{3} \beta\_{12} DCB\_{r,12} \\ \delta b\_{r,j}^{s} &= \mu\_{j} \beta\_{12} (\delta DPB\_{12}^{s} + \delta DPB\_{r,12}) - (\delta b\_{Ir}^{s} \mathbf{1}\_{2} + \delta b\_{r,IF12}) \\ IFCB &= (\delta b\_{3}^{s} - \delta b\_{Ir12}^{s} - \mu\_{3} \beta\_{12} DDB\_{12}^{s}) + (\delta b\_{r,3} - \delta b\_{r,IF12} - \mu\_{3} \beta\_{12} \delta DPB\_{r,12}) \\ \lambda\_{i} \overline{N}\_{r,i}^{s} &= \lambda\_{i} N\_{r,i}^{s} + \overline{b}\_{r,i} + \overline{b}\_{i}^{s} - d\_{i,IF12}^{s} - d\_{r,IF12} + \mu\_{i} \beta\_{12} DCB\_{r,12} \end{aligned} \tag{5}$$

where *u<sup>s</sup> <sup>r</sup>* is the directional cosine of the receiver–satellite linkage and *x* is the 3D coordinate correction value; *Zr* is the wet troposphere delay at the zenith path with a mapping function *m<sup>s</sup> <sup>r</sup>*; *IFBr* is the inter-frequency bias; and *δb<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the combined time-varying part of the unparameterized satellite and receiver-side phase hardware deviations, the effect of which can be ignored owing to its small magnitude [7]. Moreover, *dtr* is the estimated receiver clock error; *I s <sup>r</sup>*,1 is the estimated slant ionospheric delay at frequency *<sup>f</sup>*1; and *<sup>N</sup><sup>s</sup> <sup>r</sup>*,*<sup>i</sup>* is the estimated carrier phase ambiguity. ALL estimated parameters in our PPP models with IFCB are listed as:

$$X = \left[ \mathbf{x}, d\overline{\mathbf{f}}\_r, Z\_{r\prime} \overline{\mathbf{f}}\_{r,1\prime}^s, IFB\_{r\prime} \overline{\mathbf{N}}\_{r,1\prime}^s, \overline{\mathbf{N}}\_{r,2\prime}^s, \overline{\mathbf{N}}\_{r,3}^s \right] \tag{6}$$

#### *2.2. IFCB Estimation Method*

In this work, IFCB was estimated by using a difference for two ionosphere-free combinations, where the ionosphere-free combined carrier observations of L1 and L2, and L1 and L5 for *L<sup>s</sup> <sup>r</sup>*,*IF*<sup>12</sup> and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*IF*13, respectively, are

$$\begin{array}{c} L\_{r,IF12}^{s} = \rho\_r^s + dt\_r - dt^s + T\_r^s + \lambda\_{IF12} N\_{r,IF12}^s + b\_{r,IF12} + b\_{IF12}^s + \zeta\_{r,IF12}^s\\ L\_{r,IF13}^{s} = \rho\_r^s + dt\_r - dt^s + T\_r^s + \lambda\_{IF13} N\_{r,IF13}^s + b\_{r,IF13} + b\_{IF13}^s + \zeta\_{r,IF13}^s \end{array} \tag{7}$$

Note that the errors in the above equation between the receiver and the satellite antenna, phase winding, etc., have been corrected by the proposed model. Essentially, the difference between *L<sup>s</sup> <sup>r</sup>*,*IF*<sup>12</sup> and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*IF*<sup>13</sup> yields a combination of triple-frequency geometryfree and ionosphere-free phase observations (GFIF).

$$\begin{array}{lcl}GFIF & \equiv L\_{r,IF12}^s - L\_{r,IF13}^s \\ & = N\_{GFIF} + B\_{r,GFIF}^s + \delta B \end{array} \tag{8}$$

where *NGFIF* is the carrier phase ambiguity of the GFIF combination, and *B<sup>s</sup> <sup>r</sup>*,*GFIF* is the combination of the constant parts of pseudorange and phase hardware delays for GFIF combination. Meanwhile, *δB* is the IFCB value of the ionosphere-free combination, and its relationship with the non-combined *IFCB* is given as:

$$
\delta B = \frac{f\_3^2}{f\_1^2 - f\_3^2} IFCB \tag{9}
$$

In order to eliminate the ambiguity *NGFIF* and the constant term *B<sup>s</sup> <sup>r</sup>*,*GFIF* of phase hardware delay, the difference between epoch elements is calculated in the continuous observation arc without cycle slip:

$$
\Delta \delta B\_{(t, t-1)} = GFIF\_t - GFIF\_{t-1} \tag{10}
$$

where <sup>Δ</sup>*δB*(*t*,*t*−1) is the variation between the epochs at time *<sup>t</sup>* and *<sup>t</sup>* − 1. A weighted average of <sup>Δ</sup>*δB*(*t*,*t*−1) is usually calculated for the multiple stations of same epoch to improve the data stability and avoid the occasionality of the solution. For *n* stations, the following expression can be written:

$$
\Delta\delta\overline{B}\_{(t,t-1)} = \frac{\sum\_{r=1}^{n} \Delta\delta B\_{(t,t-1)} \cdot \omega\_{r,(t,t-1)}}{\sum\_{r=1}^{n} \omega\_{r,(t,t-1)}} \tag{11}
$$

where <sup>Δ</sup>*δB*(*t*,*t*−1) is the weighted average of the variation, and *<sup>ω</sup>r*,(*t*,*t*−1) is the corresponding weight of each station, as expressed in Equation (12).

$$
\omega\_{r,(t,t-1)} = \begin{cases}
0 & E\_{r(t,t-1)}^s < 10^\circ \\
\sin E\_{r(t,t-1)}^s & 10^\circ \le E\_{r(t,t-1)}^s < 30^\circ \\
1 & 30^\circ \le E\_{r(t,t-1)}^s
\end{cases}
\tag{12}
$$

where *E<sup>s</sup> <sup>r</sup>*(*t*,*t*−1) <sup>=</sup> *<sup>E</sup><sup>s</sup> r*(*t*)+*E<sup>s</sup> <sup>r</sup>*(*t*−1) <sup>2</sup> and *<sup>E</sup><sup>s</sup> <sup>r</sup>*(*t*) are the elevation angles formed by the station *r* and the satellite *s* at moment *t*. It is worthwhile to note that the use of a segmented elevation angle weighting method also reduces the impact of errors caused by the low elevation angles, while avoiding the chance of high elevation angles caused by the fewer moments of measuring stations. Therefore, the value of *δBt* at moment *t* is:

$$
\delta B\_t = \delta \overline{B}\_0 + \sum\_{k=1}^{k=t} \Delta \delta \overline{B}\_{(k,k-1)} \tag{13}
$$

In this study, *δB*<sup>0</sup> was set to 0 to introduce a common deviation for all epoch elements, which has no effect on the PPP floating solution, and the common deviation will be absorbed into the fuzziness parameter. Nevertheless, while conducting PPP-AR, the same benchmark needs to be added to the deviation products of the third frequency to avoid *δB*<sup>0</sup> effect [14].

#### **3. Experiment and Analysis**
