BDS-3/Galileo Time and Frequency Transfer with Quad-Frequency Precise Point Positioning

: In this work, quad-frequency precise point positioning (PPP) time and frequency transfer methods using Galileo E1/E5a/E5b/E5 and BDS-3 B1I/B3I/B1C/B2a observations were proposed with corresponding mathematical models. In addition, the traditional dual-frequency (BDS-3 B1I/B3I and Galileo E1/E5a) ionospheric-free (IF) model was also described and tested for comparison. To assess the proposed method for time transfer, datasets selected from timing labs were utilized and tested. Moreover, the number of Galileo or BDS-3 satellites, pseudorange residuals, positioning accuracy and tropospheric delay at receiver end were all analyzed. The results showed that the proposed quad-frequency BDS-3 or Galileo PPP models could be used to time transfer, due to stability and accuracy identical to that of dual-frequency IF model. Furthermore, the quad-frequency models can provide potential for enhancing the reliability and redundancy compared to the dual-frequency time transfer method.

Currently, numerous studies focus on Galileo [10] and BDS-3 positioning performance [11][12][13]. The positioning accuracy using pseudorange observations is better 10 m for Galileo-only [10] and BDS-3-only [11,12]. For high-precise positioning, for example, PPP and real-time kinematic (RTK) positioning, the accuracy all can reach the centimeter level [10,[14][15][16]. By the end of 2018, the timing accuracy of BDS-3 has officially announced to reach better than 20 ns [11]. For Galileo timing service, the maximum tolerable error (MTE) of service levels 3 is about 10 ns for the current and needs of users [17].

Methods
In this section, we shall start with general observations models for BDS-3 or Galileo quad-frequency signals. Following is that four BDS-3 or Galileo PPP time transfer models are then proposed. We end with detailed analysis of four different models.

General Observations
The quad-frequency un-combined pseudorange and carrier phase observations can be written as [27,28]          p s r,1 = g s r · x + cdt r + m s r · Z w + I s r,1 + β 12 · cDCB s 12 + cd r,1 + ε s r,1 p s r,2 = g s r · x + cdt r + m s r · Z w + γ 2 · I s r,1 − α 12 · cDCB s 12 + cd r,2 + ε s r,2 p s r,3 = g s r · x + cdt r + m s r · Z w + γ 3 · I s r,1 − cd s IF 12 + c(d r,3 + d s 3 ) + ε s r,3 p s r,4 = g s r · x + cdt r + m s r · Z w + γ 4 · I s r,1 − cd s IF 12 + c(d r,4 + d s 4 ) + ε s r,4 (1) l s r,1 = g s r · x + cdt r + m s r · Z w − I s r,1 − cd s IF 12 + λ 1 (N s r,1 + b r,1 + b s 1 ) + ξ s r,1 l s r,2 = g s r · x + cdt r + m s r · Z w − γ 2 · I s r,1 − cd s IF 12 + λ 2 (N s r,2 + b r,2 + b s 2 ) + ξ s r,2 l s r,3 = g s r · x + cdt r + m s r · Z w − γ 3 · I s r,1 − cd s IF 12 + λ 3 (N s r,3 + b r,3 + b s 3 ) + ξ s r,3 l s r,4 = g s r · x + cdt r + m s r · Z w − γ 4 · I s r,1 − cd s IF 12 + λ 4 (N s r,4 + b r,4 + b s 4 ) + ξ s r,4 (2) where S and r represent satellite and receiver, respectively; 1, 2, 3, and 4 are BDS-3 (B1I/B3I/B1C/B2a) or Galileo (E1/E5a/E5b/E5) quad-frequency signals. l and p are the observed-minus-computed (OMC) values of carrier-phase and pseudorange observations, respectively. g s r is the unit vector; x refers to the vector of the receiver coordinate increments in three components. c indicates the speed of light; dt r is the receiver clock offset; m s r indicates the wet mapping function; Z w demonstrates the zenith wet delay; γ j is the frequency-dependent multiplier factor (γ j = ( f 1 / f j ) 2 , j = 2, 3, 4). I s r,1 is the slant ionospheric delay on the first frequency f s 1 . α mn and β mn are the frequency factors (m = n). d s m and d r,m are the uncalibrated code bias (UCD) (m = 1, 2, 3, 4) at satellite and receiver end, respectively; d s IF mn is the ionospheric-free UCD at satellite s. DCB s mn and DCB s r,mn are the differential code bias at satellite and receiver end. λ j refers to the wavelength; N s r,j indicates the integer ambiguity; b r,j and b s j are phase delay at the receiver and satellite end, respectively. ε s r,j and ξ s r,j are the measurement noise of pseudorange and carrier phase observations. Here, m and n are the different frequency. In addition, in order to keep the receiver offset in multi-frequency PPP model consistent with that of dual-frequency, we use the receiver offset in dual-frequency ionospheric-free (IF) model as a reference.

Dual-Frequency IF PPP Model
The dual-frequency IF PPP is commonly employed for PPP model. Here, BDS-3 (B1I and B3I) or Galileo (E1 and E5a) were utilized to generate the dual-frequency IF PPP model, which is called IF0 in our study. The IF0 model can be expressed as [27,29]: with where p s r,IF12 and l s r,IF12 are the IF combination pseudorange and carrier phase OMC values, respectively; cdt s r,IF12 is the receiver clock offset, which has absorbed the ionospheric-free combination UCD at receiver end; N s r,IF12 indicates the float ambiguity; ε s r,IF12 and ξ s r,IF12 are the ionospheric-free measurement noise. Here, 1 and 2 represent BDS-3 (B1I, B1I) and Galileo (E1, E5a). The estimated parameters of IF0 model can be expressed as
Unlike IF0, three dual-frequency IF model will generate three receiver clock offsets. Note that the receiver clock offset has absorbed the UCD at receiver end (see Equation (3)). In order to separate cdt r,IF12 , the receiver clock of E1/E5b and E1/E5 combination will be divided into cdt r,IF12 and an inter-frequency (IFB) parameter. The method will simplify measurement of UCDs for different frequency at receiver end. User just calibrates the UCD of IF0 model at receiver end. Then, the linearized observation equations can be expressed as [30]:      p s r,IF12 = g s r · x + cdt r,IF12 + m s r · Z w + ε s r,IF12 p s r,IF13 = g s r · x + cdt r,IF12 + m s r · Z w + Ω s r,IF13 + ε s r,IF13 p s r,IF14 = g s r · x + cdt r,IF12 + m s r · Z w + Ω s r,IF14 + ε s r,IF14 For BDS-3, two dual-frequency combination (B1I/B3I and B1C/B2a) are recommended by the Interface Control Document (ICD) [31]. As we point out, an IFB parameter needs to be added. Thus, quad-frequency BDS-3 IF model, called IF1 (BDS-3) in our study, will generate by two dual-frequency IF combinations, and can be expressed as l s r,IF12 = g s r · x + cdt r,IF12 + m s r · Z w + N s r,IF12 + ξ s r,IF12 l s r,IF34 = g s r · x + cdt r,IF12 + m s r · Z w + N s r,IF34 + ξ s r,IF34 where p s r,IF34 and l s r,IF34 are the B1C/B2a IF combination pseudorange and carrier phase OMC values, respectively; Ω s r,IF34 is the IFB for B1C/B2a IF combination; N s r,IF34 is the redefinition ambiguity for B1C/B2a IF combination; note further that 1, 2, 3, and 4 represent B1I, B3I, B1C, and B2a signals for BDS-3 here.
Then, the estimated parameters of IF1 (BDS-3) can be written as Quad-frequency IF model, namely IF2, can also be generated by two triple-frequency combinations. The three combination coefficients need to meet three conditions (see Equation (5)), consists of unchanged geometric range, elimination of first-order ionospheric delays, and minimum noise [32]. where e 1 , e 2 , and e 3 are the combination coefficients. e 1 , e 2 , and e 3 can be calculated by the Lagrange Equation [18], the results can be described as [33]: An IFB parameter should be added in IF2 model. Then, the IF2 model for BDS-3 (B1I/B3I/B1C and B1I/B3I/B1C) or Galileo (E1/E5a/E5b and E1/E5a/E5) can be written as [32]: p s r,IF123 = e 1 p s r,1 + e 2 p s r,2 + e 3 p s r,3 = g s r · x + cdt r,IF123 + m s r · Z w + ε s r,IF123 p s r,IF124 = e 11 p s r,1 + e 22 p s r,2 + e 44 p s r,4 = g s r · x + cdt r,IF123 + m s l s r,IF123 = e 1 l s r,1 + e 2 l s r,2 + e 3 l s r,3 = g s r · x + cdt r,IF123 + m s r · Z w + N s r,IF123 + ξ s r,IF123 l s r,IF124 = e 11 l s r,1 + e 22 l s r,2 + e 44 l s r,4 = g s r · x + cdt r,IF123 + m s r · Z w + N s r,IF124 + ξ s r,IF124 where cdt r,IF123 refers to the receiver clock; Ω s r,IF124 is the IFB parameter; N s r,IF123 and N s r,IF124 are the redefinition ambiguity for two triple-frequency IF models, respectively. e 11 , e 22 , and e 44 is the combination coefficients. The estimates parameters of IF2 can be obtained as Not that 1, 2, 3, and 4 represent B1I, B3I, B1C, and B2a signals for BDS-3 or E1, E5a, E5b, and E5 for Galileo here.

Quad-Frequency Uncombined PPP Model
For quad-frequency uncombined PPP model, namely UC in this work, the slant ionospheric delay will be generally estimated as a parameter. Combined Equations (1) and (2), the precise quad-frequency uncombined PPP model can be described as [29,30]: Remote Sens. 2021, 13, 2704 where Ω s r, 3 and Ω s r,4 are the corresponding IFB parameters; I s r,UC1234 and N s r,j are the redefined slant ionospheric delay and float ambiguity parameters, respectively. Then, the unknown parameters of UC model can be expressed as:

Characteristic of Quad-Frequency GNSS PPP Time Transfer Models
To compare the IF0, IF1, IF2, and UC models with BDS-3 or Galileo observations, their major characteristics are concluded in Tables 1 and 2, respectively, including the observations used, the combination coefficients, and the noise amplification factor.
From Table 1, we can see that the B1C/B2a presents smaller noise amplification than B1I/B3I or triple-frequency combination. In addition, the noise amplification of triple-frequency combination is better than that of B1I/B3I. Interestingly, B1I/B3I/B2a combination is analogous to B1I/B2a combination due to low contribution of B3I. For Galileo quad-frequency PPP time transfer models (see, Table 2), we can see that the noise amplification of E1/E5a, E1/E5b, and E1/E5 show the similar performance, while is slightly larger than E1/E5a/E5b and E1/E5a/E5. More interestingly, E1 presents the Remote Sens. 2021, 13, 2704 7 of 26 greatest contribution for E1/E5a/E5b and E1/E5a/E5 combination. Note further that the combination coefficient retains at least 8 decimal points, otherwise it will result in an abnormal combined observation, although we keep only 4 decimal points in Table 2. The reader can calculate the combination coefficient by Equation (16).
On the other hand, the IF1 and UC models exhibit more flexible than the IF2 model. We do not recommend users to apply IF2 model for PPP time transfer. That can be explained by two reasons. For one thing, the IF2 model will be not used due to the absence of particular frequency. For another thing, users need to calibrate UCDs at three frequencies. As we know, the calibration of UCD is complex. For IF1 and UC model, UCD at two frequencies are only calibrated.

Results and Discussion
Processing strategies and datasets are first introduced. Subsequently, the quadfrequency BDS-3 and Galileo PPP time transfer solutions with different models are provided. Finally, the observation residuals, positioning, tropospheric delay, and the IFB are presented.  Table 3. Note that high-performance H-master are available at all stations. The precise orbit and clock products (WUM) with intervals of 15 min and 30 s, respectively, are released by Wuhan University, China. The BDS-3-and Galileo-only quad-frequency PPP time transfer models are studied and analyzed. The detailed processing strategies for BDS-3-or Galileo-only quad-frequency PPP time transfer models are described in Table 4. Note that the GAMP [34] software was developed to meet the requirements of BDS-3 and Galileo quad-frequency PPP-derived time transfer experiments in our work.   7.93, 6.32) for GPS, GLONASS, Galileo, and BDS-3, respectively. The number of BDS-3 satellites is a little less than that of Galileo at current state.    Figure 3 presents the time difference of XIA3-BRCH with BDS-3 IF0, IF1, IF2, and UC PPP models. Note that the difference has absorbed the hardware delay. The time series of IF0, IF1, and UC match each other very well. However, there is a clear system bias between IF2 and other models, which reflects the UCD at receiver end (see Equation (19)). In addition, the trends of time series obtained by different models are nearly consistent. For a clearer surface of our findings, we enlarge some of results in Figure 3 and present them in The frequency stability can be used to verify precise time transfer performance in another way. Here, the modified Allan deviation (MDEV) method is employed to assess the frequency stability. Figure 5 Figure 3 presents the time difference of XIA3-BRCH with BDS-3 IF0, IF1, IF2, and UC PPP models. Note that the difference has absorbed the hardware delay. The time series of IF0, IF1, and UC match each other very well. However, there is a clear system bias between IF2 and other models, which reflects the UCD at receiver end (see Equation (19)). In addition, the trends of time series obtained by different models are nearly consistent. For a clearer surface of our findings, we enlarge some of results in Figure 3 and present them in that the performance of quad-frequency PPP is equal to or better than that of dual-frequency PPP both for short-and long-term frequency stability. This finding can also be proved by Figure 6 (enlarged figure), which exhibits MDEV between different BDS−3 PPP models (DOY 15) on XIA3-BRCH time-link.   that the performance of quad-frequency PPP is equal to or better than that of dual-frequency PPP both for short-and long-term frequency stability. This finding can also be proved by Figure 6 (enlarged figure), which exhibits MDEV between different BDS−3 PPP models (DOY 15) on XIA3-BRCH time-link.   The frequency stability can be used to verify precise time transfer performance in another way. Here, the modified Allan deviation (MDEV) method is employed to assess the frequency stability. Figure 5 exhibits the modified Allan deviation (MDEV) values of the time differences, obtained from four kinds of BDS-3 precise time transfer models between BRCH and XIA3 stations at different time intervals. As shown in Figure 5 We can conclude that the performance of quad-frequency PPP is equal to or better than that of dual-frequency PPP both for short-and long-term frequency stability. This finding can also be proved by Figure 6 (enlarged figure), which exhibits MDEV between different BDS−3 PPP models (DOY 15) on XIA3-BRCH time-link.       Figures 9 and 10, respectively, for clearly presentation. Note that the receivers at PT11 and PTBB are connected to the common clock at PTB. We can conclude two findings here. First, similar to the above conclusions, Galileo IF0, IF1, and UC agree well with each other. Second, IF2 time series is different from the other models, which is relative to the UCD. To further study the characteristic of Galileo precise time transfer model, we calculated mean and STD values of time series obtained from IF0, IF1, IF2, and UC models and listed in Tables 5 and 6 , and UC models, respectively. Therefore, we obtain the second conclusion that quad-frequency Galileo precise time transfer models show the same performance. In addition, the quad-frequency Galileo PPP time transfer performance present a slightly better than that of BDS−3. There are several possible explanations for this result. First, the accuracy of BDS−3 orbit and clock products need to be further improved. Second, the TDOP values of Galileo satellite system is less than that of BDS-3. Third, the DCB values of BDS-3 we are currently used is just a preliminary file and has not yet been officially released. Besides, an unknown system bias exists in the IF2 model, so we do not recommend users to use this model for time transfer (See the mean values).

Quad-Frequency Galileo PPP Time Transfer Solutions
ble explanations for this result. First, the accuracy of BDS−3 orbit and clock products need to be further improved. Second, the TDOP values of Galileo satellite system is less than that of BDS-3. Third, the DCB values of BDS-3 we are currently used is just a preliminary file and has not yet been officially released. Besides, an unknown system bias exists in the IF2 model, so we do not recommend users to use this model for time transfer (See the mean values).                  Figure 15 illustrates the residuals of pseudorange observations at BRCH and XIA3 stations of BDS-3 IF0, IF1, IF2, and UC PPP models on DOY 16. Figure 16 demonstrates the residuals of pseudorange observations at BRUX and PT11 stations of Galileo IF0, IF1, IF2, and UC PPP models on DOY 278. Four findings can be concluded here. First, the pseudorange residuals are basically close to zero. It suggests that our DCB estimation and correction method is right. Second, the residuals of UC model are smaller than ionospheric-free models. It seems possible that these results are due to low noise amplification (see Tables 1 and 2). Third, for BDS-3 pseudorange residuals may be affected by the preliminary DCB file. Fourth, the pseudorange residuals of Galileo quad-frequency precise time transfer is obviously smaller than that of BDS-3. These differences can be explained  Figure 15 illustrates the residuals of pseudorange observations at BRCH and XIA3 stations of BDS-3 IF0, IF1, IF2, and UC PPP models on DOY 16. Figure 16 demonstrates the residuals of pseudorange observations at BRUX and PT11 stations of Galileo IF0, IF1, IF2, and UC PPP models on DOY 278. Four findings can be concluded here. First, the pseudorange residuals are basically close to zero. It suggests that our DCB estimation and correction method is right. Second, the residuals of UC model are smaller than ionosphericfree models. It seems possible that these results are due to low noise amplification (see Tables 1 and 2). Third, for BDS-3 pseudorange residuals may be affected by the preliminary DCB file. Fourth, the pseudorange residuals of Galileo quad-frequency precise time transfer is obviously smaller than that of BDS-3. These differences can be explained by the fact that the orbit and clock products, DCB values and number of BDS-3 satellites need to be further improved.

Positioning, Tropospheric Delay, and IFB Estimates
For time users, the station coordinates are commonly estimated as a constant or constrained to improve receiver clock offset estimation. To fully study the quad-frequency precise time transfer method, the selected station coordinates were estimated as unknow parameters. Figures 17 and 18 exhibit the positioning error for different BDS-3 or Galileo PPP models on XIA3 and BRUX stations, respectively. Note that the precise receiver coordinates of PT11, BRCH, XIA3 are calculated by Bernese 5.2 software with IGS final products to assess the positioning accuracy. To further quantify positioning accuracy, the root mean square (RMS) of positioning error for BRCH and XIA3 stations with BDS-3 different PPP models and for BRUX, PT11, PTBB with Galileo different PPP models are listed in Tables 7 and 8 (7)-(22)), which are related to the receiver after satellite DCB correction.
Here, we will just show IFB parameters obtained from UC model at XIA3 and BRCH, which depicts at Figure 21. The mean STD values of IFB values are (0.018, 0.006) ns for BDS-3 UC models at BRCH and XIA3 stations, respectively. Hence, we can conclude that the IFB parameters are very stable during the whole day.   For Galileo IF1, IF2, and UC models, two, one and two additional parameters also need to be estimated. As previous mentioned, Figure 22 exhibits the IFB parameters time series obtained from Galileo UC models. The mean STD values of IFB parameters are (0.014, 0.004) ns for UC models at BRUX and PT11 stations. Similar to the IFB parameters of BDS-3 PPP model, the IFB values vary stable, with the small STD values.

Conclusions
PPP time transfer with Galileo and BDS-3 quad-frequency observations was proposed in this work. Three quad-frequency PPP time and frequency transfer models, namely IF1, IF2, and UC, were studied and developed. The mathematical models of quad-frequency PPP were first described in detailed. Then, datasets from timing lab were selected and employed to evaluate reliability and feasibility of the proposed models. With the experimental results, three findings can be concluded.
First, the proposed three quad-frequency PPP models are still utilized to precise time transfer with BDS-3-only and Galileo-only quad-frequency observations. The reason for this is the fact that the frequency stability and time transfer accuracy are identical to or better than the corresponding IF0 PPP solutions.
Second, BDS-3-only or Galileo-only IF1, IF2, and UC models agree well with each other due to the same performance of positioning accuracy and tropospheric delay. However, we do not recommend time users to use IF2 models. This result may be explained by the fact that an unknown system bias exists in the IF2 model. Besides, IF2 model will be not used due to the absence of particular frequency. In addition, the receiver DCB values obtained from BDS-3 or Galileo quad-frequency PPP models vary stable.
Third, the Galileo quad-frequency PPP solutions perform slightly better than that of BDS-3. This result is explained by the fact that the orbit and clock products need to be further improved.
Currently, a few stations can receive all BDS-3 and Galileo quad-frequency signals. With the development of Galileo, BDS-3, system and GNSS tracking station, BDS-3 + Galileo combination PPP time and frequency with quad-frequency observations will be studied. In addition, PPP time and frequency transfer with ambiguity resolution of quadfrequency observations is of interest for future study. Unlike IF models, UC model provides not only station coordinates, receiver clock offset, troposphere delay, but also ionosphere information. Both BDS-3 or Galileo IF and UC models can be employed not only for positioning, troposphere delay, but also for precise time transfer.
Author Contributions: Y.G. and X.C. conceived and designed the experiments; Y.G. performed the experiments, analyzed the data, and wrote the paper; F.S., X.Y. and S.W. helped in the discussion and revision. All authors have read and agreed to the published version of the manuscript.