Galileo Time Transfer with Five-Frequency Uncombined PPP: A Posteriori Weighting, Inter-Frequency Bias, Precise Products and Multi-Frequency Contribution

Abstract

Galileo satellites can broadcast signals on five frequencies, namely E1, E5A, E5B, E5 (A+B), and E6. The multi-frequency integration has become an emerging trend in Global Navigation Satellite System (GNSS) data processing. This study focused on the precise time transfer based on Galileo five-frequency uncombined precise point positioning (PPP), including the performance comparison of PPP time transfer with a priori and a posteriori weighting strategies, with different inter-frequency bias (IFB) dynamic models, and with the precise satellite products from different analysis centers, as well as the contribution of multi-frequency observations for time transfer. Compared with the a priori weighting strategy, the short-term frequency stability of time transfer adopting the Helmert variance component estimation method can be improved by 28.9–37.6% when the average time is shorter than 100 s. The effect of IFB dynamic models on Galileo five-frequency PPP time transfer is not obvious. When employing the post-processed precise satellite products from seven analysis centers, the accuracy of time transfer can be better than 0.1 ns, while an accuracy of 0.253 ns can be obtained in the real-time mode. At an average time of approximately 10,000 s, the post-processed time transfer with Galileo five-frequency PPP can provide a frequency stability of 3.283 × 10⁻¹⁴ to 3.459 × 10⁻¹⁴, while that in real-time mode can be 3.541 × 10⁻¹⁴. Compared with dual-frequency PPP results, the contribution of multi-frequency combination to the accuracy and frequency stability of time transfer is not significant, but multi-frequency PPP can achieve more reliable time transfer results when the signal quality is poor.

1. Introduction

Global Navigation Satellite System (GNSS) is one of the main tools for precise time transfer. In comparison with traditional tools, GNSS time transfer has the advantages of low cost, small size, all-weather availability, and easy maintenance. Since 1980, many researchers have focused on the GNSS technology for time transfer, and several typical technologies have been proposed, including common view (CV), all-in view (AV), two-way satellite time and frequency transfer (TWSTFT), and carrier phase (CP) methods [1]. The basic idea, advantage, and disadvantage of the four GNSS time transfer methods are listed in

Table 1. Basic idea, advantage, and disadvantage of four typical GNSS time transfer technologies.

Technology	Basic Idea	Accuracy	Advantage	Disadvantage
CV	Utilizing observations of enough common-view satellites from two stations, and cancelling out or weakening many errors through single-difference between stations	1–10 ns	Low cost and simple data processing	The time transfer accuracy degrades with the increasing distance between the two stations, and is also limited by the precision of code observations
AV	Measuring the time differences between two stations with the code observations from all observed satellites as well as the precise satellite orbit and clock offset products	Better than CV method	Low cost and free of the negative effects of the distance between stations	The implementation needs to rely on the precise satellite orbit and clock offset products, and the time transfer accuracy is limited by the precision of code observations
TWSTFT	Exchanging time signals between two ground stations through a geostationary satellite, and cancelling out many delay variations of the transmitted signals due to the reciprocity of the propagation path	Higher than AV method	High accuracy, and high long-term frequency stability	This method needs the support of communication satellites and the employed equipment must have the ability to communicate with satellites, which limits its applications
CP	Employing the precise carrier phase observations as well as the precise satellite orbit and clock offset products to conduct the time transfer through the precise point positioning technology	A sub-nanosecond level	Low cost, free of the negative effects of the distance between stations, high accuracy, and high short-term frequency stability	The implementation needs to rely on the precise satellite orbit and clock offset products

. The CP method, which is also called the precise point positioning (PPP) method, combines the advantages of the other three methods. The short-term frequency stability of PPP time transfer is higher than that of TWSTFT method [2]. Currently, PPP has become the primary method to maintain the International Atomic Time (TAI) by Bureau International des Poids et Mesures (BIPM) and to keep the time synchronization by time laboratories in various countries [3,4].

For a long period, the PPP time transfer technology mainly relied on GPS. The frequency stability of GPS PPP time transfer with ambiguity resolution (AR) could reach a level of 1 × 10⁻¹⁶ and sub−10⁻¹⁶ at an average time of five days and seven days, respectively [5,6]. In recent years, with the overall deployment of Multi-GNSS Experiment (MGEX) stations, the PPP time transfer based on GLONASS-only, Galileo-only, BDS-only, and multi-GNSS integrated measurements has been a research focus. Ge et al. [7] investigated the effects of three different inter-frequency bias (IFB)-handling strategies on the time transfer performance with GLONASS-only dual-frequency PPP, and the results indicated that the standard deviation (STD) statistics of the differences between receiver clock estimates derived from GLONASS-only and GPS-only PPP could be 0.2 ns when estimating IFB for each GLONASS satellite. Xu et al. [8] evaluated the contribution of Galileo PPP AR on the frequency stability of time transfer. The results indicated that the ambiguity-fixed solutions had slightly reduced short-term stability compared with that of the ambiguity-float solutions, but the long-term stability of the ambiguity-fixed solutions was significantly improved. Zhang et al. [9] compared the time transfer performance of GPS-only, Galileo-only, BDS-2-only, and GPS/Galileo/BDS-2 combined dual-frequency PPP, and it was indicated that the STD statistics of the short time link could be 0.059, 0.116, 0.089, and 0.055 ns for the four different constellation combination cases, respectively. Ge et al. [10,11] further jointly used the single- or dual-frequency observations from GPS, GLONASS, Galileo, and BDS-2 satellites to carry out the time transfer. The four-system integrated PPP could improve the frequency stability of time links by 9.61% and 20.3% for single- and dual-frequency cases over GPS-only PPP, respectively. Lyu et al. [12] analyzed the time transfer accuracy of GPS/GLONASS/Galileo/BDS-2 combined PPP with fixed ambiguities for GPS, and it was proved that the STD statistics of the time link were decreased by 38.1% over GPS-only ambiguity-float PPP.

As most of the modernized GNSS satellites start to broadcast multi-frequency signals, many researchers investigated the time transfer performance of multi-frequency combined PPP. Tu et al. [13] proved that the BDS-2 triple-frequency uncombined (UC) PPP model could be applied to the time transfer, and the obtained accuracy and frequency stability of time links were comparable to those of the traditional dual-frequency ionospheric-free (IF) model. With the use of BDS-3 observations, Su and Jin [14] derived similar conclusions. Zhang et al. [15] rigorously evaluated the time transfer performance of five Galileo four-frequency IF combined models. The average STD statistics of time link differences between four-frequency and triple-frequency PPP models, and between four-frequency and dual-frequency PPP models could be 0.041 ns, but they exhibited consistency in terms of frequency stability of time links. Ge et al. [16] compared the time transfer performance of Galileo and BDS-3 four-frequency PPP. The results indicated that the frequency stability of Galileo time transfer results was slightly better than that of BDS-3, and both the accuracy and frequency stability of them were comparable to those of the corresponding dual-frequency PPP model.

According to the above research status, the precise time transfer with multi-frequency PPP has attracted an increasing attention from the GNSS community, but most of the existing studies were confined to triple- or four-frequency observations. Actually, Galileo satellites can provide signals on five frequencies, namely E1, E5A, E5B, E5 (A+B), and E6, and BDS-3 satellites can transmit five open-service signals, namely B1I, B3I, B1C, B2a, and B2b. The current multi-frequency PPP time transfer models need to be extended so as to incorporate the five-frequency observations. Most importantly, the existing studies did not carefully consider the differences of multi-frequency observations. Taking Galileo five-frequency observations as an example, the signals on different frequencies differ in data quality, anti-jamming capability, hardware delay, and measurement noises because they employ different modulation modes. The E1 signal is modulated using Composite Binary Offset Carrier (CBOC) modulation, while the Alternative Binary Offset Carrier (AltBOC) modulation is used for the E5 signal (E5A/E5B/E5 (A+B)) [17,18,19]. Thus, the precision of observations differs among different frequencies. Most of the existing studies lost sight of this and set the same variances for the measurements on all the frequencies, which is obviously unreasonable.

This study focused on the precise time transfer based on the Galileo five-frequency UC PPP, including the effects of weighting strategies and IFB dynamic models as well as the contribution of multi-frequency combination. In addition, there are seven analysis centers (ACs) in support of post-processed or real-time Galileo precise satellite products at present. The systematic biases contained in the precise satellite products differ among different ACs, as they employ different strategies and models during precise orbit determination and precise clock estimation. These biases will be partially absorbed by the receiver clock parameters during PPP data processing, and thus affect the time transfer performance. Therefore, we also compared the time transfer performance of Galileo five-frequency UC PPP with the precise products from different ACs. In Section 2, the Galileo five-frequency UC PPP model for time transfer is first described, and then the Helmert variance component estimation method in Kalman filter is introduced. Subsequently, Section 3 presents the results, including the accuracy and frequency stability of time transfer derived from two different weighting strategies, three different IFB dynamic models, eight different precise satellite products, and four different inter-frequency combination cases. Next, the discussion is conducted in Section 4. Finally, we summarize the main points and the conclusions in Section 5.

2. Methods

2.1. Galileo Five-Frequency Uncombined PPP Model for Time Transfer

The single-frequency code and carrier phase observations from a Galileo receiver can be described below:

$$\{\begin{array}{l}{P}_i=\rho +c\mathrm{d}{t}_r-c\mathrm{d}t+I\cdot {\gamma }_i+T+d_{r,i}+d_i\hfill \\ {\Phi }_i=\rho +c\mathrm{d}{t}_r-c\mathrm{d}t-I\cdot {\gamma }_i+T+N_i+b_{r,i}+b_i\hfill \end{array}$$

(1)

with ${\gamma }_i=f_1^2/f_i^2$ where $i$ refers to the Galileo E1, E5A, E5B, E6, and E5 (A+B) frequency when it is equal to 1, 2, 3, 4, and 5, respectively, $P$ is the measured pseudorange, $\Phi$ is the measured carrier phase, $\rho$ is the geometric range, $c\mathrm{d}{t}_r$ is the receiver clock offset, $c\mathrm{d}t$ is the satellite clock offset, $I$ is the ionospheric path delay on the E1 frequency, $\gamma$ is the ionospheric delay factor, $f$ is the carrier frequency, $T$ is the tropospheric path delay, $N$ is the phase ambiguity, $d_r$ and $d$ are the receiver- and satellite-dependent code hardware delay, respectively, and $b_r$ and $b$ refer to the phase-specific hardware delay at the receiver and satellite, respectively. The four hardware delay terms are considered to be stable within a short period of time, for example, a day.

In precise clock estimation (PCE), the common clock is usually assumed for both code and carrier phase observation equations, but the code hardware delay is ignored. Consequently, the satellite clock parameters will absorb the time-invariant satellite code hardware delay. The current Galileo satellite clocks, which are computed with E1/E5A IF combined code and carrier phase observations, can be expressed as:

$$c\mathrm{d}{t}_{IF,12}=c\mathrm{d}t-(a_{12,1}\cdot d_1+a_{12,2}\cdot d_2)$$

(2)

with $a_{12,1}=f_1^2/(f_1^2-f_2^2)$ and $a_{12,2}=-f_2^2/(f_1^2-f_2^2)$.

When the PCE is carried out based on different observations, the estimated satellite clocks will be different. The satellite clock estimates are frequency dependent due to the frequency-dependent code hardware delays lumped into them. Thus, the current E1/E5A IF satellite clocks provided by the international GNSS service (IGS) cannot be directly applied to the Galileo multi-frequency PPP processing. The E1/E5A IF satellite clocks can be converted into the desired satellite clocks for the UC observables on a specific frequency by combining them with the differential code biases (DCBs), that is:

$$\{\begin{array}{l}c\mathrm{d}{t}_{UC,1}=c\mathrm{d}{t}_{IF,12}-a_{12,2}\cdot \mathrm{DCB}(P_1,P_2)=c\mathrm{d}t-d_1\hfill \\ c\mathrm{d}{t}_{UC,2}=c\mathrm{d}{t}_{IF,12}+a_{12,1}\cdot \mathrm{DCB}(P_1,P_2)=c\mathrm{d}t-d_2\hfill \\ c\mathrm{d}{t}_{UC,3}=c\mathrm{d}{t}_{IF,12}-a_{12,2}\cdot \mathrm{DCB}(P_1,P_2)+\mathrm{DCB}(P_1,P_3)=c\mathrm{d}t-d_3\hfill \\ c\mathrm{d}{t}_{UC,4}=c\mathrm{d}{t}_{IF,12}-a_{12,2}\cdot \mathrm{DCB}(P_1,P_2)+\mathrm{DCB}(P_1,P_4)=c\mathrm{d}t-d_4\hfill \\ c\mathrm{d}{t}_{UC,5}=c\mathrm{d}{t}_{IF,12}-a_{12,2}\cdot \mathrm{DCB}(P_1,P_2)+\mathrm{DCB}(P_1,P_5)=c\mathrm{d}t-d_5\hfill \end{array}$$

(3)

with

$$\{\begin{array}{l}\mathrm{DCB}(P_1,P_2)=d_1-d_2\hfill \\ \mathrm{DCB}(P_1,P_3)=d_1-d_3\hfill \\ \mathrm{DCB}(P_1,P_4)=d_1-d_4\hfill \\ \mathrm{DCB}(P_1,P_5)=d_1-d_5\hfill \end{array}$$

(4)

where $\mathrm{DCB}(P_1,P_2)$, $\mathrm{DCB}(P_1,P_3)$, $\mathrm{DCB}(P_1,P_4)$ and $\mathrm{DCB}(P_1,P_5)$ denote the DCBs between E1 and E5A, between E1 and E5B, between E1 and E6, and between E1 and E5 (A+B), respectively, and $c\mathrm{d}{t}_{UC,1}$, $c\mathrm{d}{t}_{UC,2}$, $c\mathrm{d}{t}_{UC,3}$, $c\mathrm{d}{t}_{UC,4}$ and $c\mathrm{d}{t}_{UC,5}$ denote the satellite clocks on a single frequency, namely the E1, E5A, E5B, E6, and E5 (A+B) UC satellite clocks, respectively. The four sets of DCB products have been provided by the Chinese Academy of Sciences (CAS) [20]. As shown in Equation (3), in addition to the frequency-independent physical clock error of the satellite (i.e., $c\mathrm{d}t$), the UC satellite clocks only contain the single-frequency satellite code hardware delays on the corresponding frequency. As such, the satellite code hardware delays in the UC code observation equation (see Equation (1)) can be removed by applying these UC satellite clocks.

The ionospheric delay in the UC-based PPP model is taken as unknown parameter. As both the hardware delay and ionospheric delay are frequency dependent, only their combination rather than the separated term is unbiasedly estimable because of rank deficiency. Due to the perfect correlation between the receiver code hardware delay and ionospheric delay for the dual-frequency PPP with E1 and E5A UC observations, we usually estimate them as a lumped term in the dual-frequency UC model. As to the Galileo five-frequency UC PPP, the ionospheric delay estimate cannot completely absorb the effect of receiver code hardware delay on the E5B, E6, and E5 (A+B) frequencies. Therefore, three additional IFB parameters need to be respectively introduced into the E5B, E6, and E5 (A+B) code observation equations. Regarding the tropospheric delays, we usually employ an a priori model to correct the dry parts, while the wet parts are estimated as unknowns from the observations. After applying the E1, E5A, E5B, E6, and E5 (A+B) UC satellite clocks as well as the satellite orbits, the linearized observation model for Galileo five-frequency UC PPP can be given below:

$$\{\begin{array}{l}{p}_1=\mu \cdot X+c\mathrm{d}{t}_{r,E}+I_E\cdot {\gamma }_1+m\cdot Z\hfill \\ p_2=\mu \cdot X+c\mathrm{d}{t}_{r,E}+I_E\cdot {\gamma }_2+m\cdot Z\hfill \\ p_3=\mu \cdot X+c\mathrm{d}{t}_{r,E}+I_E\cdot {\gamma }_3+m\cdot Z+g_{UC-3}\hfill \\ p_4=\mu \cdot X+c\mathrm{d}{t}_{r,E}+I_E\cdot {\gamma }_4+m\cdot Z+g_{UC-4}\hfill \\ p_5=\mu \cdot X+c\mathrm{d}{t}_{r,E}+I_E\cdot {\gamma }_5+m\cdot Z+g_{UC-5}\hfill \\ {\phi }_1=\mu \cdot X+c\mathrm{d}{t}_{r,E}-I_E\cdot {\gamma }_1+m\cdot Z+N_{1,E}\hfill \\ {\phi }_2=\mu \cdot X+c\mathrm{d}{t}_{r,E}-I_E\cdot {\gamma }_2+m\cdot Z+N_{2,E}\hfill \\ {\phi }_3=\mu \cdot X+c\mathrm{d}{t}_{r,E}-I_E\cdot {\gamma }_3+m\cdot Z+N_{3,E}\hfill \\ {\phi }_4=\mu \cdot X+c\mathrm{d}{t}_{r,E}-I_E\cdot {\gamma }_4+m\cdot Z+N_{4,E}\hfill \\ {\phi }_5=\mu \cdot X+c\mathrm{d}{t}_{r,E}-I_E\cdot {\gamma }_5+m\cdot Z+N_{5,E}\hfill \end{array}$$

(5)

with

$$\{\begin{array}{l}c\mathrm{d}{t}_{r,E}=c\mathrm{d}{t}_r+a_{12,1}\cdot d_{r,1}+a_{12,2}\cdot d_{r,2}\hfill \\ I_E=I-a_{12,2}\cdot (d_{r,2}-d_{r,1})\hfill \\ g_{UC-3}=(-{\gamma }_3\cdot a_{12,2}-a_{12,1})\cdot d_{r,1}-a_{12,2}\cdot (1-{\gamma }_3)\cdot d_{r,2}+d_{r,3}\hfill \\ g_{UC-4}=(-{\gamma }_4\cdot a_{12,2}-a_{12,1})\cdot d_{r,1}-a_{12,2}\cdot (1-{\gamma }_4)\cdot d_{r,2}+d_{r,4}\hfill \\ g_{UC-5}=(-{\gamma }_5\cdot a_{12,2}-a_{12,1})\cdot d_{r,1}-a_{12,2}\cdot (1-{\gamma }_5)\cdot d_{r,2}+d_{r,5}\hfill \\ N_{1,E}=N_1+b_{r,1}+b_1-d_1-(a_{12,1}-{\gamma }_1\cdot a_{12,2})\cdot d_{r,1}-a_{12,2}\cdot (1+{\gamma }_1)\cdot d_{r,2}\hfill \\ N_{2,E}=N_2+b_{r,2}+b_2-d_2-(a_{12,1}-{\gamma }_2\cdot a_{12,2})\cdot d_{r,1}-a_{12,2}\cdot (1+{\gamma }_2)\cdot d_{r,2}\hfill \\ N_{3,E}=N_3+b_{r,3}+b_3-d_3-(a_{12,1}-{\gamma }_3\cdot a_{12,2})\cdot d_{r,1}-a_{12,2}\cdot (1+{\gamma }_3)\cdot d_{r,2}\hfill \\ N_{4,E}=N_4+b_{r,4}+b_4-d_4-(a_{12,1}-{\gamma }_4\cdot a_{12,2})\cdot d_{r,1}-a_{12,2}\cdot (1+{\gamma }_4)\cdot d_{r,2}\hfill \\ N_{5,E}=N_5+b_{r,5}+b_5-d_5-(a_{12,1}-{\gamma }_5\cdot a_{12,2})\cdot d_{r,1}-a_{12,2}\cdot (1+{\gamma }_5)\cdot d_{r,2}\hfill \end{array}$$

(6)

where $p$ and $\phi$ are the observed-minus-computed (OMC) code and carrier phase observables, respectively, $\mu$ is the unit vector in the direction from the receiver to the satellite, $X$ is the three-dimensional (3D) receiver positions, $c\mathrm{d}{t}_{r,E}$ is the receiver clock estimates, which are a lumped term of frequency-independent receiver physical clocks and receiver code hardware delays on E1 and E5A frequencies combined within the E1/E5A IF combination, $I_E$ is the slant ionospheric delay estimate on the E1 frequency, which is not the real ionospheric delay because of the lumped receiver code hardware delays on E1 and E5A frequencies, $m$ is the wet mapping function, $Z$ is the tropospheric zenith wet delay (ZWD), $g_{UC-3}$, $g_{UC-4}$, and $g_{UC-5}$ are the receiver IFBs on the E5B, E6, and E5 (A+B) frequencies, respectively, all of which are a linear combination of receiver code hardware delays on E1 and E5A frequencies as well as those on a third frequency (E5B, E6, or E5 (A+B)), and $N_E$ is the phase ambiguity estimates, which are non-integers because of the integrated code and phase hardware delays of the satellite and receiver.

The parameters which need estimation for Galileo five-frequency UC PPP include the receiver 3D coordinates, receiver clocks, slant ionospheric delays on the E1 frequency, ZWD, receiver IFBs on the E5B, E6, and E5 (A+B) frequencies, and E1, E5A, E5B, E6, and E5 (A+B) UC float ambiguities, that is:

$$S_{UC}={\left[X,c\mathrm{d}{t}_{r,E},I_E,Z,g_{UC-3},g_{UC-4},g_{UC-5},N_{1,E},N_{2,E},N_{3,E},N_{4,E},N_{5,E}\right]}^{\mathrm{T}}$$

(7)

where $S_{UC}$ is the state vector for the Galileo five-frequency UC PPP.

To optimize the time transfer performance with Galileo five-frequency UC PPP, the proper stochastic models should also be determined. Following Satirapod and Luansang [21], the satellite elevations with a mapping form of trigonometric function or exponential function can be used to quantify the precision of satellite observations. The elevation-dependent weighting scheme with a Sine mapping function was adopted in this study. It was assumed that the observations in the E1, E5A, E5B, E6, and E5 (A+B) carriers are uncorrelated. In addition, they have the same a priori noises. The covariance matrix of UC code and phase observations can be computed as:

$$R_{UC}=\mathrm{diag}\left({\sigma }_P^2,{\sigma }_P^2,{\sigma }_P^2,{\sigma }_P^2,{\sigma }_P^2,{\sigma }_{\Phi }^2,{\sigma }_{\Phi }^2,{\sigma }_{\Phi }^2,{\sigma }_{\Phi }^2,{\sigma }_{\Phi }^2\right)$$

(8)

with ${\sigma }_P=k_1/\sin (el)$ and ${\sigma }_{\Phi }=k_2/\sin (el)$ where both $k_1$ and $k_2$ are constants, which are set to 0.3 m for code observations and 0.003 m for phase observations, respectively, and $el$ is the satellite elevation angle. For the Galileo five-frequency UC PPP model, the covariance matrix is a diagonal matrix.

Based on the established observation model and stochastic model, the Kalman filter is employed to estimate the unknown parameters (state vector) by the following equation:

$$\{\begin{array}{l}{S}_k=S_{k,k-1}+K_k(L_k-B_k{S}_{k,k-1})\\ D_k=(I-K_k{B}_k)D_{k,k-1}\\ K_k=D_{k,k-1}{B}_k^{\mathrm{T}}{(B_k{D}_{k,k-1}{B}_k^{\mathrm{T}}+R_k)}^{-1}\end{array}$$

(9)

where $S_{k,k-1}$ and $D_{k,k-1}$ are the predicted state vector and its covariance matrix at epoch k, respectively, and $S_k$ and $D_k$ are the corresponding filtering results. $L_k$ and $R_k$ are the observation vector and its covariance matrix, respectively, $I$ is the unit matrix, $B_k$ is the design matrix, and $K_k$ is the Kalman gain matrix.

2.2. Helmert Variance Component Estimation Method for Kalman Filter

Obviously, the a priori weighting method (see Equation (8)) is not rigorous. In the Helmert variance component estimation (HVCE) method, the a posteriori observation residuals are used to determine the variance components of various observations so that the influence of different types of GNSS observations on the solutions can be effectively balanced [22,23]. In this study, based on the a priori weighting method, the HVCE method was further introduced into the Galileo five-frequency UC PPP data processing to determine the reasonable variances of observations on different frequencies.

Considering that the noise level is different between code and carrier phase observations as well as among the observations on different frequencies, the code and carrier phase observations on five frequencies are grouped into ten types to conduct the HVCE. To obtain the reliable HVCE time transfer solutions, there must be enough observations for each type and the difference of the number of observations among the ten types should not be significant. Based on this classification principle, the residuals and covariance matrix of various observations after carrying out the Kalman filter can be expressed as:

$$\{\begin{array}{c}{V}_{P_1}=B_{P_1}\widehat{S}-L_{P_1}\\ ⋮\\ V_{P_5}=B_{P_5}\widehat{S}-L_{P_5}\\ V_{{\Phi }_1}=B_{{\Phi }_1}\widehat{S}-L_{{\Phi }_1}\\ ⋮\\ V_{{\Phi }_5}=B_{{\Phi }_5}\widehat{S}-L_{{\Phi }_5}\end{array}\hspace{1em}$$

(10)

$$R=\left[\begin{array}{cccccc}{R}_{P_1}& \cdots & 0& 0& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & R_{P_5}& 0& \cdots & 0\\ 0& \cdots & 0& R_{{\Phi }_1}& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & R_{{\Phi }_5}\end{array}\right]$$

(11)

where $P_i$ and ${\Phi }_i$ denote the code and phase observations on the $i\mathrm{th}\hspace{0.17em}(i=1,2,\cdots ,5)$ frequency for all visible satellites, respectively, $R_{P_i}$ and $R_{{\Phi }_i}$ are the corresponding covariance matrix, $L$, $V$, and $B$ denote the observation vector, observation residual error vector, and design matrix, respectively, and $\widehat{S}$ denotes the estimated parameters with Kalman filter. As the HVCE method requires iterative processing, to omit the complex matrix operation and improve the calculation efficiency, the simplified HVCE calculation formula is used to estimate the variance components of all types of observations [24], that is:

$$\{\begin{array}{l}{\widehat{\sigma }}_{P_i}^2=\frac{V_{P_i}^{\mathrm{T}}{R}_{P_i}^{-1}{V}_{P_i}}{n_{P_i}-tr(N^{-1}{N}_{P_i})}\\ {\widehat{\sigma }}_{{\Phi }_i}^2=\frac{V_{{\Phi }_i}^{\mathrm{T}}{R}_{{\Phi }_i}^{-1}{V}_{{\Phi }_i}}{n_{{\Phi }_i}-tr(N^{-1}{N}_{{\Phi }_i})}\end{array}$$

(12)

with

$$\{\begin{array}{l}N={\displaystyle {\sum }_{i=1}^5(N_{P_i}+}{N}_{{\Phi }_i})+D_{\overline{S}}^{-1}\hfill \\ N_{P_i}=B_{P_i}^{\mathrm{T}}{R}_{P_i}^{-1}{B}_{P_i}\hfill \\ N_{{\Phi }_i}=B_{{\Phi }_i}^{\mathrm{T}}{R}_{{\Phi }_i}^{-1}{B}_{{\Phi }_i}\hfill \end{array}$$

(13)

where ${\widehat{\sigma }}_{P_i}^2$ and ${\widehat{\sigma }}_{{\Phi }_i}^2$ denote the variance components of code and phase observations on the $i\mathrm{th}$ frequency, respectively, $n_{P_i}$ and $n_{{\Phi }_i}$ denote the corresponding number of observations, and $D_{\overline{S}}^{}$ denotes the covariance matrix of the predicted state vector.

Based on the obtained variance components, the variances of observations can be re-computed as:

$$\{\begin{array}{l}{R}_{P_i,m}=\frac{{\widehat{\sigma }}_{P_i,m}^2}{c_{0,m}}\cdot R_{P_i,m-1}\\ R_{{\Phi }_i,m}=\frac{{\widehat{\sigma }}_{{\Phi }_i,m}^2}{c_{0,m}}\cdot R_{{\Phi }_i,m-1}\end{array}$$

(14)

where the subscripts $m$ and $m-1$ denote the times of iterations, and $c_0$ is an arbitrary constant, which is usually taken as one of the estimated variance components.

Figure 1 presents the flow chart of the HVCE method for Kalman filter. The computation procedure of the HVCE method includes the following steps:

Step 1: Obtain the state vector and observation residual error vector by performing Kalman filter with a priori covariance matrix at an epoch.

Step 2: Calculate the variance components of all types of observations based on Equation (12).

Step 3: Update the covariance matrix of observations based on Equation (14), and then rerun the Kalman filter to obtain a posteriori residuals once more.

Step 4: Repeat Step 2 and Step 3 until the difference of the maximum value and minimum value among the derived variance components of ten types of observations is smaller than 0.001, namely $\max \{{\widehat{\sigma }}_{P_i}^2,{\widehat{\sigma }}_{{\Phi }_i}^2\}-\min \{{\widehat{\sigma }}_{P_i}^2,{\widehat{\sigma }}_{{\Phi }_i}^2\}<0.001$.

3. Results

3.1. Datasets and Processing Strategies

To assess the time transfer performance of the Galileo five-frequency UC PPP model, the datasets from a week on day of year (DOY) 248 to 254, 2021, collected at 18 MGEX stations were used for static data processing. The sampling interval was 30 s, and the cut-off elevation angle was set to 10°. Following Prochniewicz and Grzymala [17], both the code and phase multipath errors (including measurement noises) of Galileo significantly increased with decreasing elevation angles when the elevation angles were lower than 20°. After considering a compromise between signal quality and satellite availability (note that there were only 24 operational Galileo satellites during the analysis period), the widely used threshold of 10° was taken as the cut-off elevation here. The geographical distribution of 18 MGEX tracking stations is shown in Figure 2. All stations were equipped with high-performance hydrogen clocks and could support the Galileo five-frequency observations. The BRUX station operated by a time laboratory at the Royal Observatory of Belgium was chosen as the central node for the time transfer, and a total of 17 time links were designed. Considering that the accuracy of the post-processed receiver clock offset products provided by IGS can reach 75 ps, they were taken as the reference values when evaluating the accuracy of time transfer results. Since the PPP time transfer technology has a significant advantage in short-term (shorter than one day) frequency stability, the data processing was carried out day-by-day. To analyze the effects of different precise products on PPP time transfer, eight precise satellite products in support of the Galileo system, which are identified by “COD”, “GBM”, “WUM”, “SHA”, “IAC”, “ESA”, “GRG”, and “CNT”, over the analysis period were collected. It should be noted that the Galileo precise satellite products were not provided by “SHA” on DOY 250, 2021. These precise products were generated by the analysis centers Center for Orbit Determination in Europe (CODE), GeoForschungsZentrum Potsdam (GFZ), Wuhan University (WHU), Shanghai Observatory (SHAO), Information and Analysis Center (IAC), European Space Agency (ESA), and Centre National d’Etudes Spatiales (CNES), and their information is detailed in

Table 2. Galileo precise satellite products from different ACs.

Institution	ID	Orbit	Clock
CODE	COD0MGXFIN	5 min	30 s
GFZ	GBM0MGXRAP	5 min	30 s
WHU	WUM0MGXFIN	15 min	30 s
SHAO	SHA0MGXRAP	5 min	5 min
IAC	IAC0MGXFIN	5 min	30 s
ESA	ESA0MGNFIN	5 min	30 s
CNES	GRG0MGXFIN	5 min	30 s
	CNT	5 min	5 s

. It is worth noting that CNES provides two sets of precise products, namely “GRG” and “CNT”, among which “GRG” is the post-processed precise products, and “CNT” is the real-time precise products retrieved based on real-time data stream CLK93 and broadcast ephemeris.

The sufficient visible satellites and the smaller time dilution of precision (TDOP) values are beneficial for the high-accuracy time transfer. As of the end of 2021, the Galileo constellation had a total of 24 operational satellites, including three in-orbit validation (IOV) satellites and 21 Full Operational Capability (FOC) satellites [25]. Figure 3 and Figure 4 provide the geographical distribution of the average number of visible satellites and average TDOP values for Galileo over the entire analysis period, respectively. There are seven visible Galileo satellites (on average) in areas spanning the latitudes 15°N to 55°N and 15°S to 55°S, while the satellite numbers are slightly increased to eight (on average) in other areas. Due to the evenly distributed visible satellites, the change of average TDOP values on a global scale is not significant, with a varying range of 1.1 to 1.4. Therefore, the current number of Galileo satellites is sufficient and the TDOP value is satisfactory, which can meet the requirements of PPP time transfer.

3.2. Time Transfer Results with HVCE Method

Four time links (i.e., the difference of receiver clock estimates at two different stations) between BRUX and KOKB, between BRUX and ONSA, between BRUX and SPT0, and between BRUX and PTBB on DOY 251, 2021, were selected to compare the time transfer performance of the Galileo five-frequency PPP with a priori weighting method and the HVCE a posterior weighting method, and the results are illustrated in Figure 5. The post-processed precise products “GBM” were used for PPP processing. It can be seen that the estimates of time links showed a significant difference between the two different weighting schemes in the first 8 h. After the filter fully converged, the corresponding difference of the time link results was not obvious. Most importantly, compared with the a priori weighting method, the time links of the HVCE method converged faster. The results indicate that the adjustment of weights of observations with HVCE method is helpful to the filter convergence, and thus is expected to improve the accuracy of PPP time transfer within the first several hours.

The frequency stability, which can be expressed by Modified Allan Deviation (MDEV), is an important index to evaluate the time transfer performance. Figure 6 displays the MDEV of four time links with two different weighting schemes on DOY 251, 2021, and

Table 3. MDEV values at several specific average times of four time links with the a priori weighting method and HVCE method on DOY 251, 2021.

Average Time (s)	Modified Allan Deviation
	BRUX-KOKB		BRUX-ONSA		BRUX-SPT0		BRUX-PTBB
	A Priori	HVCE	A Priori	HVCE	A Priori	HVCE	A Priori	HVCE
30	9.650 × 10−13	7.006 × 10−13	1.061 × 10−12	5.090 × 10−13	9.162 × 10−13	7.169 × 10−13	1.263 × 10−12	6.434 × 10−13
60	8.007 × 10−13	5.006 × 10−13	5.588 × 10−13	2.888 × 10−13	4.422 × 10−13	3.540 × 10−13	5.899 × 10−13	4.074 × 10−13
120	5.182 × 10−13	4.828 × 10−13	3.386 × 10−13	2.789 × 10−13	3.640 × 10−13	2.356 × 10−13	4.455 × 10−13	3.470 × 10−13
240	2.050 × 10−13	1.233 × 10−13	1.353 × 10−13	1.202 × 10−13	1.048 × 10−13	6.975 × 10−14	2.137 × 10−13	1.712 × 10−13
480	6.794 × 10−14	4.809 × 10−14	8.353 × 10−14	5.872 × 10−14	6.722 × 10−14	3.617 × 10−14	1.008 × 10−13	9.569 × 10−14
960	3.444 × 10−14	3.174 × 10−14	3.901 × 10−14	2.849 × 10−14	4.491 × 10−14	3.296 × 10−14	5.091 × 10−14	4.818 × 10−14
1920	1.969 × 10−14	2.080 × 10−14	1.677 × 10−14	1.419 × 10−14	2.201 × 10−14	1.957 × 10−14	1.796 × 10−14	1.550 × 10−14
3840	1.239 × 10−14	1.369 × 10−14	9.514 × 10−15	8.338 × 10−15	9.108 × 10−15	9.473 × 10−15	8.254 × 10−15	8.604 × 10−15
7680	7.280 × 10−15	7.627 × 10−15	5.543 × 10−15	5.028 × 10−15	4.300 × 10−15	5.054 × 10−15	6.399 × 10−15	5.832 × 10−15
15,360	1.110 × 10−14	1.114 × 10−14	4.651 × 10−15	4.398 × 10−15	3.835 × 10−15	3.738 × 10−15	3.901 × 10−15	3.716 × 10−15
28,800	1.079 × 10−14	9.842 × 10−15	3.168 × 10−15	2.016 × 10−15	1.978 × 10−15	3.568 × 10−16	1.730 × 10−15	2.396 × 10−15

lists the MDEV values at several specific average times. It is clearly seen that the frequency stability of the HVCE scheme was better than that of the a priori weighting scheme for all the four time links, especially for the short-term frequency stability. Based on the average results over the four time links, the HVCE scheme can improve the frequency stability by 28.9%–37.6% when the average time is shorter than 100 s (including 30, 60, and 90 s). With the increasing average time, the corresponding improvement rate of HVCE scheme decreased. When the average time fell between 1000 and 10,000 s, the frequency stability of the HVCE scheme was even slightly lower than that of the a priori weighting scheme sometimes for the time links between BRUX and KOKB, between BRUX and SPT0, and between BRUX and PTBB. For the long-term (with an average time longer than 10,000 s) frequency stability, the HVCE scheme outperformed the a priori weighting scheme, but the difference between them was confined to a level of 10⁻¹⁶ to 10⁻¹⁵. With the increment of average time, the number of samples for MDEV calculation was continuously reduced, which may have resulted in the increasing MDEV errors [26]. In addition, the difference of time link results between the HVCE scheme and a priori weighting scheme was not significant after the filter converged. Thus, the frequency stability at an average time longer than 1000 s seems to be at a similar level for the two different weighting schemes. The results indicate that the HVCE scheme can also be beneficial to improve the short-term frequency stability of time transfer compared with the a priori weighting scheme.

3.3. Time Transfer Results with Different IFB Dynamic Models

There are three IFB parameters in the Galileo five-frequency UC PPP observation model, and thus a reasonable IFB dynamic model should be taken to ensure the reliability of the receiver clock offsets estimated through Kalman filter. The post-processed precise products “GBM” are used for analysis. To investigate the time correlation of IFB parameters among the epochs, the IFBs are modeled as constants, random walk process, and white noise process for comparison to conduct the time transfer with Galileo five-frequency PPP. If IFBs are very stable and will not change over time, they should be estimated as constants. In the case of IFBs with a high time correlation, they should be estimated as random walk process. If the IFBs are independent among different epochs, they should be estimated as white noise process. Figure 7 illustrates the time series of receiver clock offset estimates with three different IFB dynamic models at BRUX and SPT0 stations on DOY 248, 2021. It can be seen that the receiver clock estimates were consistent for the three IFB dynamic models (with a maximum difference smaller than 0.01 ns). With the use of receiver clock estimates at the two stations, a time link can be generated. We can also obtain the corresponding time link using the IGS final products. Taking the time link derived from IGS final products as references, the epoch-wise time transfer errors are shown in Figure 8. The effects of different time datums are eliminated by subtracting the average value. The time transfer errors with a varying range of 0.25 ns (excluding the results within the first several minutes) were also consistent for the three IFB dynamic models (with a maximum difference smaller than 0.01 ns). The STD statistics can be used to assess the time transfer accuracy. The STD values of time link differences were 0.041 ns for all the three IFB dynamic models (excluding the results within the first four hours at a converging stage). To investigate the frequency stability, Figure 9 further provides the MDEV results of the time link BRUX-SPT0. Similarly, both the short-term and long-term frequency stability were basically consistent for the three IFB dynamic models. When the IFBs were modeled as constants, random walk process, and white noise process, the frequency stability at an average time of 90 s was 2.777 × 10⁻¹³, 2.778 × 10⁻¹³, and 2.769 × 10⁻¹³, respectively, while the corresponding stability at an average time of 9990 s was 4.066 × 10⁻¹⁵, 4.064 × 10⁻¹⁵, and 4.050 × 10⁻¹⁵, respectively. Therefore, the effects of IFB dynamic models on the precise time transfer with Galileo five-frequency UC PPP were marginal. For simplicity, the random walk process was employed for IFB estimation in the following data processing.

3.4. Time Transfer Results with Precise Products from Different Analysis Centers

In order to evaluate the time transfer performance of different precise products, the Galileo five-frequency PPP processing with eight different precise products was carried out at 18 selected MGEX stations to generate 17 time links. Figure 10 illustrates the time transfer results of three time links BRUX-SPT0, BRUX-KOKB, and BRUX-USN7 on DOY 248, 2021. The time links of seven post-processed precise products were consistent with each other, which indicates that all the post-processed precise products in support of Galileo can be used for time transfer. The time transfer results of the real-time product “CNT” were significantly different from those of the post-processed ones, especially for the time link BRUX-USN7, in which the maximum differences could reach 0.5 ns. Figure 11 shows the differences of the time links between Galileo five-frequency PPP estimates with different precise products and IGS final products on DOY 248, 2021, and

Table 4. STD values of differences of the time links between PPP estimates and IGS final products on DOY 248, 2021.

Products	STD (ns)
	BRUX-SPT0	BRUX-KOKB	BRUX-USN7	Mean
CNT	0.054	0.388	0.317	0.253
COD	0.071	0.136	0.056	0.088
ESA	0.066	0.126	0.046	0.079
GBM	0.066	0.095	0.076	0.079
GRG	0.069	0.103	0.045	0.072
IAC	0.074	0.154	0.070	0.099
SHA	0.070	0.114	0.054	0.079
WUM	0.068	0.139	0.041	0.083

provides the corresponding STD statistics. To ensure that only the fully converged results are involved in the analysis, the results in the first four hours were excluded. After removing the effects of different time datums, the time transfer accuracy of post-processed products (usually better than 0.2 ns) was significantly higher than that of real-time products (usually better than 0.6 ns). According to the average STDs over the three time links, the product “GRG” achieved the highest time transfer accuracy with a statistic of 0.072 ns, while that of the product “IAC” was the worst among the post-processed products, which was 0.099 ns. As to the real-time product “CNT”, the time transfer accuracy degraded to only 0.253 ns.

The frequency stability of time transfer with various precise products was further analyzed. Figure 12 illustrates the MDEV of the time link BRUX-KOKB derived with eight different precise products on DOY 248, 2021. In terms of the short-term frequency stability, the MDEV of post-processed product “SHA” was 5.802 × 10⁻¹³ at an average time of 90 s, while that of real-time product “CNT” was 4.844 × 10⁻¹³, which is better than the post-processed product “SHA” but slightly worse than the other six post-processed products. This may be because the sampling interval of “SHA” precise satellite clock product was as long as 5 min, which is much longer than that for other post-processed products (30 s) and real-time products (5 s). When employing the product “SHA” for PPP processing, the interpolation errors of satellite clock offsets will be larger. Due to the difference in accuracy between real-time products and post-processed products, the frequency stability of “CNT” results gradually deteriorated with the increase of average time, while that of “SHA” results became close to other post-processed products. In terms of the long-term frequency stability (e.g., at an average time of 10,000 s), the “CNT” results were obviously worse than those of the seven post-processed products.

Figure 13 exhibits the mean MDEV values for 17 time links derived from Galileo five-frequency PPP with eight different precise products over a week from DOY 248 to 254, 2021. For simplicity, the numbers of 10², 10³, and 10⁴ in the figure refer to the average time of 90, 990, and 9990 s, respectively. When the average time was less than 1000 s, the frequency stability of “SHA” post-processed results was significantly lower than that of other precise products for most of the cases, even in comparison with “CNT” real-time results, while that of “CNT” results was slightly worse than the other six post-processed products. In contrast, regarding the long-term frequency stability, the results of real-time product “CNT” were slightly worse than those of all the seven post-processed products.

Table 5. Mean MDEV values over all the time links with eight different precise products at an average time of 90, 990, and 9990 s.

Products	Mean Modified Allan Deviation
	90 s	990 s	9990 s
CNT	6.648 × 10−13	8.460 × 10−14	3.541 × 10−14
COD	6.125 × 10−13	7.630 × 10−14	3.367 × 10−14
ESA	6.257 × 10−13	7.529 × 10−14	3.283 × 10−14
GBM	6.235 × 10−13	7.661 × 10−14	3.301 × 10−14
GRG	6.185 × 10−13	7.562 × 10−14	3.296 × 10−14
IAC	6.212 × 10−13	7.814 × 10−14	3.406 × 10−14
SHA	8.539 × 10−13	1.152 × 10−13	3.459 × 10−14
WUM	6.207 × 10−13	7.881 × 10−14	3.355 × 10−14

further lists the mean MDEV values over all the time links with eight different precise products at an average time of 90, 990, and 9990 s. When the average time was 90 s, the average MDEV of “SHA” post-processed results was 8.539 × 10⁻¹³, which is much larger than that of other precise products, while that of “CNT” real-time results was 6.648 × 10⁻¹³. Except for the above two products, the difference of frequency stability among the other six post-processed precise products was confined to 1.4 × 10⁻¹⁴. The average MDEV of “COD” results was the smallest with a value of 6.125 × 10⁻¹³, while that of “ESA” results was the largest, which was 6.257 × 10⁻¹³. When the average time was 990 s, there was still a big gap between “SHA” and other precise products in terms of MDEV statistics. When the average time was further increased to 9990 s, the frequency stability of “ESA”, “GBM”, and “GRG” results was slightly better than that of the other four post-processed precise products, which was smaller than 3.301 × 10⁻¹⁴, while the frequency stability of real-time product “CNT” was the worst, with an average MDEV of 3.541 × 10⁻¹⁴.

According to the above analysis, in terms of time transfer accuracy, the results of real-time product “CNT” showed a worse performance compared with those of post-processed products, while the corresponding accuracy difference among the seven post-processed products was not obvious. As to frequency stability, both the short-term and long-term performance of “CNT” real-time results were slightly worse than those of the post-processed products (except for “SHA”). Due to the lower sampling rate of satellite clock offset products, the “SHA” post-processed results had poor short-term frequency stability but satisfied long-term frequency stability. The short-term frequency stability of “COD”, “GRG”, and “WUM” products was slightly better than that of other post-processed products, while the “ESA”, “GBM”, and “GRG” products slightly outperformed other post-processed products in terms of the long-term frequency stability. Therefore, when the Galileo five-frequency PPP is used for short-term or long-term time transfer, these precise products can be given priority.

3.5. Time Transfer Results with Galileo Dual-, Triple-, Four- and Five-Frequency UC PPP

To investigate the contribution of Galileo multi-frequency observations to PPP time transfer, the E1/E5A dual-frequency (UC-2), E1/E5A/E5B triple-frequency (UC-3), E1/E5A/E5B/E6 four-frequency (UC-4), and E1/E5A/E5B/E6/E5 (A+B) five-frequency (UC-5) integration were employed to conduct the Galileo UC PPP time transfer. The post-processed precise product “GBM” was used for data processing. Figure 14 illustrates three time links (BRUX-SPT0, BRUX-KOKB, and BRUX-USN7) derived from Galileo dual-, triple-, four-, and five-frequency UC PPP solutions on DOY 248, 2021. To ensure that only the fully converged results were involved in the analysis, the results in the first four hours were excluded. It was seen that there was no significant difference in the time transfer results among the four different frequency combination cases. Figure 15 shows the differences of the time links between Galileo UC PPP estimates and IGS final products on DOY 248, 2021, and

Table 6. STD value of differences of the time links between Galileo UC PPP estimates and IGS final products on DOY 248, 2021.

Strategy	STD (ns)
	BRUX-SPT0	BRUX-KOKB	BRUX-USN7
UC-2	0.066	0.090	0.075
UC-3	0.067	0.092	0.077
UC-4	0.067	0.094	0.077
UC-5	0.066	0.095	0.076

provides the corresponding STD statistics. The joint use of more frequencies did not obviously improve the accuracy of time transfer compared with the results of dual-frequency UC PPP. Figure 16 presents the MDEV of the time link BRUX-KOKB. It was indicated that there was still no significant difference in the frequency stability of time transfer among the four different frequency combination cases.

Figure 17 shows the mean MDEV values for 17 time links generated with Galileo dual-, triple-, four-, and five-frequency UC PPP over a week from DOY 248 to 254, 2021. Except for the time link BRUX-MRC1, the results of dual-, triple-, four-, and five-frequency cases showed consistent levels in terms of MDEV values.

Table 7. Mean MDEV values over 16 time links (except for BRUX-MRC1) generated with Galileo dual-, triple-, four-, and five-frequency UC PPP at an average time of 90, 990, and 9990 s.

Strategy	Mean Modified Allan Deviation
	90 s	990 s	9990 s
UC-2	4.955 × 10−13	6.198 × 10−14	2.034 × 10−14
UC-3	4.944 × 10−13	6.132 × 10−14	2.026 × 10−14
UC-4	4.983 × 10−13	6.149 × 10−14	2.032 × 10−14
UC-5	4.885 × 10−13	6.155 × 10−14	2.033 × 10−14

further lists the mean MDEV values over 16 time links (except for BRUX-MRC1) for different frequency integration cases at an average time of 90, 990, and 9990 s. In terms of short-term frequency stability (90 s), the UC-5 case slightly outperformed other multi-frequency cases, with a statistic of 4.885 × 10⁻¹³, while those of other three cases varied within the range of 4.944 × 10⁻¹³ to 4.983 × 10⁻¹³. Regarding the long-term frequency stability (9990 s), there was no significant difference among the four different frequency combination cases, all of which were close to 2.03 × 10⁻¹⁴.

For the time link BRUX-MRC1, the MDEV results of time transfer with UC-2 were dramatically worse than those of the cases employing more frequencies. It was found that there was a significant difference between the dual-frequency PPP solutions and other multi-frequency PPP solutions on DOY 249, 2021. Figure 18 shows the MDEV of the time link BRUX-MRC1 derived from Galileo dual-, triple-, four-, and five-frequency UC PPP on DOY 249, 2021. It was seen that the MDEV results of UC-3, UC-4, and UC-5 were consistent, but far different from those of UC-2. The reason for this phenomenon is that the quality of observations at station MRC1 on DOY 249, 2021 was poor, with frequent cycle slips and loss of lock of satellite signals on the E1 and E5A frequencies. This leads to the unreliable time transfer results of dual-frequency PPP, and the frequency stability at an average time of 90 and 9990 s was only 2.092 × 10⁻¹¹ and 2.160 × 10⁻¹², respectively. By contrast, the corresponding frequency stability reached approximately 6.06 × 10⁻¹² and 1.14 × 10⁻¹³ for other multi-frequency cases, which was much better than that of dual-frequency PPP.

4. Discussion

The obtained results indicated that both the accuracy and frequency stability of time transfer were at a similar level for the dual-frequency and multi-frequency PPP under an open sky environment, namely under the observing conditions with no signal blockage around the receiver. However, in the constrained visibility environments (poor signal quality), the time transfer results of multi-frequency PPP were more reliable. This is because the multi-frequency observations can be regarded as the repeated measurement for the geometric range between the receiver and the satellite, and there is a strong correlation among them. When the observation conditions were ideal, the observations on all frequencies were of high quality, and the benefits from the repeated measurement were limited. Thus, the time transfer results of multi-frequency and dual-frequency PPP were almost the same. When the signal quality degrades due to various interference (surrounding the stations), the joint use of more available observations is beneficial to the cycle slip detection and data quality control, and thus the negative effects of cycle slips and gross errors on the PPP processing can be effectively weakened. Therefore, the time transfer results of multi-frequency PPP are expected to be much better than those of dual-frequency PPP in harsh environments.

5. Conclusions

Based on 18 MGEX stations equipped with high-performance hydrogen clocks, we rigorously investigated the precise time transfer with Galileo five-frequency UC PPP, including the effects of different weighting strategies (a priori method and a posterior method), of different IFB dynamic models (white noise process, random walk process and constants), and of precise products from different ACs, as well as the contribution of multi-frequency integration. After employing the HVCE a posterior weighting method, the accuracy of time transfer in the first several hours was significantly improved, and the short-term frequency stability (at an average time shorter than 100 s) could be improved by 28.9%–37.6%. The IFB dynamic models had little influence on time transfer results. The time transfer accuracy of Galileo five-frequency UC PPP with the real-time product “CNT” was 0.253 ns, while that of all the seven post-processed products was better than 0.1 ns. In terms of the short-term frequency stability, the PPP time transfer with the post-processed product “SHA” showed the worst performance, with the MDEV of 8.539 × 10⁻¹³ at an average time of 90 s, which was even worse than that of the real-time product “CNT” (6.648 × 10⁻¹³), while that of other post-processed products ranged from 6.125 × 10⁻¹³ to 6.257 × 10⁻¹³. The PPP time transfer with the real-time product “CNT” achieved the worst long-term frequency stability, with the MDEV of 3.541 × 10⁻¹⁴ at an average time of 9990 s, which was inferior to the post-processed products (3.283 × 10⁻¹⁴ to 3.459 × 10⁻¹⁴). Usually, the time transfer results of multi-frequency PPP and dual-frequency PPP are basically consistent in terms of accuracy and frequency stability. When the signal quality becomes poor (e.g., in harsh environments), the time transfer results of multi-frequency PPP are more reliable because of the higher redundancy of observations.