100 Picosecond/Sub-10⁻¹⁷ Level GPS Differential Precise Time and Frequency Transfer

Abstract

A Global Navigation Satellite System (GNSS) is a high-precision method for comparing clocks and time transfer. The GNSS carrier phase can provide more precise observable information than pseudorange. However, the carrier phase is ambiguous, and only pseudorange can provide the absolute time difference between two clocks. In our study, by taking full advantage of GNSS pseudorange and carrier-phase observables, a differential precise time transfer (DPT) method with a clustering constraint was employed to estimate the time difference between two clocks, aiming to achieve accurate and precise time and frequency transfer. Using this method, several time transfer results were analyzed for different baselines. For the common clock experiment, the time transfer results showed good consistency across different days, with an intra-day accuracy of within 20 ps. Furthermore, we evaluated the self-consistency of DPT using closure among three stations. DPT closure of the three stations had a peak-to-peak value of closure of about 25 ps. The closure did not change over time, indicating the self-consistency of the DPT processing in time transfer. Moreover, our results were compared to station clock solutions provided by the International GNSS Service (IGS), and the standard deviations (STDs) of the four baselines were all less than 100 ps within one month, confirming the time and frequency stability of the DPT method. In addition, we found that the time stability of DPT was less than 20 ps within one week. As for frequency stability, DPT achieved a 10⁻¹⁶ level of modified Allan deviation (MDEV) at an averaging time of about 1 day and a sub-10⁻¹⁷ level at an averaging time of one week.

1. Introduction

Time transfer technology is a prerequisite for time synchronization and timing, and high-precision time and frequency transfer technology is particularly important for achieving high-precision time and frequency services. Since the early 1980s, time and frequency transfer based on Global Navigation Satellite Systems (GNSSs) has been widely developed. Satellite one-way time dissemination, as the simplest GNSS time transfer method, provides clock solutions with an accuracy of 20–30 nanoseconds [1]. The Satellite Common View (CV) method provides a relatively high accuracy in short-baseline time comparisons [2]. However, the distance and accuracy of time transfer are limited by the availability of common satellites and the target range measurement [3,4]. In 2004, Jiang et al. [5] proposed the All-in-View (AV) satellite time transfer method as an improvement over the CV method. With AV, time transfer distance is no longer limited by the availability of common satellites [6,7]. However, both methods use pseudorange, and the pseudorange multipath and noise limit the time and frequency transfer performance, such that only nanosecond-level time transfer can be achieved [8].

Considering that GNSS carrier-phase measurements are two orders of magnitude more accurate than pseudorange measurements, it is possible to further improve time transfer performance by using phase measurements. Schildknecht et al. [9] first used GPS carrier-phase measurements to validate high-precision time and frequency transfer. As a continuation of the AV method, the Precise Point Positioning (PPP) method uses both pseudorange and carrier-phase measurements for high-precision un-difference time and frequency transfer, and it has been widely promoted in the time community, becoming the main means used by the Bureau International des Poids et Mesures (BIPM) for comparing clocks contributing to International Atomic Time (TAI) techniques [10,11,12]. In addition, using carrier-phase differential methods for high-precision time and frequency transfer has become a complementary method. In 2007, J. Delporte et al. [13] directly fixed the single-difference (SD) integer ambiguity by combining AV and CV technology with SD pseudorange and carrier-phase measurements, realizing SD carrier-phase time transfer. This method is suitable for post-processing analysis. In 2008, Lee et al. [14] utilized International GNSS Service (IGS) orbit and clock products to calculate carrier-phase AV and carrier-phase CV clock solutions and realized GPS carrier-phase CV and AV time transfer. In addition, the method based on real-time kinematics (RTK) double difference has also been used to recover high-precision inter-station clock information. Feng et al. [1] proposed a four-dimensional real-time kinematics (4D-RTK) timing method, which further recovered the inter-station SD integer carrier-phase ambiguity based on double-difference observations to obtain relative clock solutions and achieved a precision of better than 0.1 ns on a 21 km baseline for a duration of 1.5 h. What is more, Tu et al. [3] further compared the time difference estimates obtained from the RTK and PPP methods using short baseline links; the STD of the time difference between the two methods was about 0.1 ns. However, while using the carrier phase to improve the precision of clock comparison, the accuracy of time transfer is often overlooked.

If the primary focus is on time transfer rather than frequency transfer, the issue of an absolute time difference in carrier-phase time transfer solutions also needs to be solved. Since the carrier phase is ambiguous, only pseudorange can provide the absolute time difference between two stations, so the method using only the carrier phase can only provide frequency transfer results. In order to provide accurate time difference results between two stations, it is necessary to introduce pseudorange to determine receiver clock offset in carrier-phase time transfer. However, as it is affected by the noise and multipath of pseudorange, it is difficult to fix the SD integer ambiguity.

For the problem of integer ambiguity resolution, Teunissen et al. [15] proposed a least-squares ambiguity decorrelation adjustment (LAMBDA) method based on integer Gaussian transformation in 1995. Since then, the LAMBDA method has been widely used for DD integer ambiguity resolution in RTK positioning [16,17,18]. In 1999, Teunissen et al. [19] further summarized integer rounding, integer bootstrapping, and integer least squares for a class of integer estimators and demonstrated the best performance of the integer least-squares estimator.

In this study, the LAMBDA method was used for the determination of double-difference (DD) ambiguities. When estimating SD ambiguities from pseudorange measurements, we employed the integer rounding estimator based on the K-means algorithm to determine their integer values and then achieved differential precise time and frequency transfer (DPT). To verify the performance of DPT, common clock time transfer experiments were analyzed and calibrated, as well as three stations’ closure time transfer results. In addition, we evaluated time and frequency transfer performances for four baselines over 30 days. The discussion and conclusion are given in the final section.

2. Differential Precise Time and Frequency Transfer

The key problem to be solved by the differential precise time and frequency transfer of GNSS is to obtain high-precision inter-station clock offset. First, based on the GNSS observations of the time reference station and the user station, the DD mathematical model is established, and the solution of DD carrier-phase ambiguities and the corresponding variance matrix are obtained. Secondly, the reference satellite is selected to calculate the receiver clock offset through the SD pseudorange observation, then the inter-station SD carrier-phase ambiguity of the reference satellite is estimated, and, further, the inter-station SD carrier-phase ambiguities of all visible satellites are estimated. Affected by pseudorange noise and multipath, there may be a bias in the estimation of SD carrier-phase ambiguity by the integer rounding estimator. Then, we determine the final SD integer carrier-phase ambiguities by the K-means clustering method [20,21]. Once the ambiguities are determined, the carrier phase can provide the “absolute” time difference between stations. Finally, the accurate time transfer based on the SD carrier-phase observation is realized. In the following, the detailed DPT method will be introduced.

2.1. Double-Difference Fixed Ambiguity Resolution

In the field of geodesy, the mathematical model for relative positioning [22] based on DD pseudorange and the carrier phase for frequency, $f$, is expressed as:

$$\left\{\begin{array}{l}{P}_{f(m,n)}^{i,j}={\rho }_{f(m,n)}^{i,j}+{eph}_{f(m,n)}^{i,j}+T_{f(m,n)}^{i,j}+I_{f(m,n)}^{i,j}+{{\epsilon }_P}_{f(m,n)}^{i,j}\\ L_{f(m,n)}^{i,j}={\rho }_{f(m,n)}^{i,j}+{eph}_{f(m,n)}^{i,j}+T_{f(m,n)}^{i,j}-I_{f(m,n)}^{i,j}+{\lambda }_f\cdot N_{f(m,n)}^{i,j}+{{\epsilon }_L}_{f(m,n)}^{i,j}\end{array}\right.$$

(1)

where $P_{f(m,n)}^{i,j}$ and $L_{f(m,n)}^{i,j}$ represent DD pseudorange and carrier-phase observations from satellites i and j to receivers m and n, respectively. ${\rho }_{f(m,n)}^{i,j}$ represents the DD geometric distance, ${eph}_{f(m,n)}^{i,j}$ represents the DD ephemeris error residual, $T_{f(m,n)}^{i,j}$ represents the DD tropospheric delay residual, $I_{f(m,n)}^{i,j}$ represents the DD ionospheric delay residual, $N_{f(m,n)}^{i,j}$ represents the DD carrier-phase integer ambiguity, and ${\lambda }_f$ is the corresponding wavelength. ${{\epsilon }_P}_{f(m,n)}^{i,j}$ and ${{\epsilon }_P}_{f(m,n)}^{i,j}$ are the sum of the multipath errors and measurement noise of the DD pseudorange and carrier-phase observations, respectively.

For relative positioning technology, when the distance between two receivers is a short baseline of up to a few kilometers, the DD tropospheric and ionospheric delay residuals can be ignored [23]. The LAMBDA method [16] is used to determine the DD integer carrier-phase ambiguity. In addition, the dual-frequency ionosphere-free (IF) combination can also be used to eliminate the impact of the first-order term of ionospheric delay for longer baselines. DD wide-lane integer ambiguities, $N_{W\left(m,n\right) }^{i,j}$, between two stations are calculated by using the Melbourne–Wübbena combination observation and then obtaining the float DD carrier-phase ambiguity, $N_{IF\left(m,n\right)}^{i,j}$, of the IF combination.

$$\begin{array}{rr}{N}_{IF\left(m,n\right)}^{i,j}& =\frac{1}{{\lambda }_1}·\left(\frac{f_1^2}{f_1^2-f_2^2}·{\lambda }_1·N_{f_1\left(m,n\right)}^{i,j}-\frac{f_2^2}{f_1^2-f_2^2}·{\lambda }_2·N_{f_2\left(m,n\right)}^{i,j}\right)\\ & =\frac{f_1}{f_1+f_2}{N}_{f_1\left(m,n\right)}^{i,j}+\frac{f_1{·f}_2}{f_1^2-f_2^2}{N}_{W\left(m,n\right)}^{i,j}\hfill \\ & \end{array}$$

(2)

where $N_{W\left(m,n\right)}^{i,j}=N_{f_1\left(m,n\right)}^{i,j}-N_{f_2\left(m,n\right)}^{i,j}$. The DD carrier-phase ambiguity of the IF combination can be estimated from IF combination observation. Furthermore, L1 DD integer ambiguity, $N_{f_1\left(m,n\right)}^{i,j}$, can be obtained.

$$N_{f_1\left(m,n\right)}^{i,j}=\left[\frac{f_1+f_2}{f_1}{N}_{IF\left(m,n\right)}^{i,j}-\frac{f_2}{f_1-f_2}{N}_{W\left(m,n\right)}^{i,j}\right]$$

(3)

where $\left[·\right]$ denotes rounding to the nearest integer.

2.2. Single-Difference Ambiguity Estimation

According to Equation (1), it is known that the DD carrier-phase ambiguity solution eliminates the receiver clock offset, which is the most crucial parameter in time and frequency transfer. To retain the clock offset, we perform inter-station SD resolution using the pseudorange and carrier-phase observations. The SD model is shown in Equation (4):

$$\left\{\begin{array}{l}{P}_{f(m,n)}^i={\rho }_{f(m,n)}^i+c{{·t}_r}_{f(m,n)}^j+{eph}_{f(m,n)}^i+T_{f(m,n)}^i+I_{f(m,n)}^i+{{\epsilon }_P}_{f(m,n)}^i\\ L_{f(m,n)}^i={\rho }_{f(m,n)}^i+c{{·t}_r}_{f(m,n)}^j+{eph}_{f(m,n)}^i+T_{f(m,n)}^i-I_{f(m,n)}^i+{\lambda }_f\cdot N_{f(m,n)}^i+{{\epsilon }_L}_{f(m,n)}^i\end{array}\right.$$

(4)

where ${t_r}_{f(m,n)}^j$ represents the inter-station SD receiver clock offset and $c$ represents the speed of light. For the SD carrier-phase observation, receiver clock offset and SD carrier-phase ambiguity are linearly dependent. To separate them, we introduce the initial value of the inter-station receiver clock offset through SD pseudorange observation and then obtain SD carrier-phase ambiguity. When DD carrier-phase ambiguities are fixed, we use the pseudorange and carrier-phase observations to calculate the float SD carrier-phase ambiguity, $N_{f(m,n)}^i(t)$, at epoch $t$. For short baselines of several kilometers, ignoring the SD ionospheric delay residual, the SD tropospheric delay residual, and the SD ephemeris residual, Equation (5) can be obtained:

$$N_{f\left(m,n\right)}^i(t)=\frac{1}{{\lambda }_f}\cdot (L_{f\left(m,n\right)}^i(t)-P_{f\left(m,n\right)}^i(t))$$

(5)

Due to the influence of the noise and multipath of pseudorange, it is generally difficult to accurately determine the SD integer ambiguity using only one epoch of pseudorange measurement. It is known that the ambiguity remains continuous for each visible satellite once the receiver is correctly locked. Therefore, we smoothed Equation (5) for $l$ epochs to reduce the effects of pseudorange noise and multipath and obtained the float SD carrier-phase ambiguity, which is shown in Equation (6). Usually, the pseudorange noise and multipath of the satellite at the highest elevation is the smallest, so the reference satellite can be selected by using the satellite elevations.

$${\overline{N}}_{f(m,n)}^i=\frac{1}{l}{\sum }_{t=1}^{t=l}{N}_{f\left(m,n\right)}^i(t)=\frac{1}{l{\cdot \lambda }_f}\cdot {\sum }_{t=1}^{t=l}(L_{f\left(m,n\right)}^i(t)-P_{f\left(m,n\right)}^i(t))$$

(6)

When using Equation (6) to determine the SD integer ambiguity, the number of epochs required is closely related to the precision of pseudorange measurement. It is generally believed that integer ambiguity can be determined when the accuracy of the average value reaches ±6 cm [24]. Then, we use the integer rounding method to obtain the SD integer carrier-phase ambiguity, $N_{f(m,n)}^i$, of the reference satellite, $i$, as shown in Equation (7).

$$N_{f(m,n)}^i=\left[{\overline{N}}_{f(m,n)}^i\right]=\left[\frac{1}{l{\cdot \lambda }_f}\cdot {\sum }_{t=1}^{t=l}(L_{f\left(m,n\right)}^i(t)-P_{f\left(m,n\right)}^i(t))\right]$$

(7)

where $\left[·\right]$ denotes rounding to the nearest integer. After obtaining the fixed DD integer ambiguity by Equation (1) and the SD integer ambiguity of reference satellite $i$ by Equation (7), it is easy to obtain SD integer ambiguity, $N_{f(m,n)}^j$, of the visible satellite, $j$, by using Equation (8).

$$N_{f(m,n)}^j=N_{f(m,n)}^i+N_{f(m,n)}^{i,j}$$

(8)

After obtaining L1 DD integer ambiguity and DD wide-lane integer ambiguities, $N_{W\left(m,n\right)}^{i,j}$, by using Equations (2) and (3), IF combination SD ambiguity, $N_{IF\left(m,n\right)}^j$, of the visible satellite, $j$, can be obtained.

$$\begin{array}{rr}{N}_{IF\left(m,n\right)}^j& =N_{IF\left(m,n\right)}^i+N_{IF\left(m,n\right)}^{i,j}\hfill \\ & =N_{IF(m,n)}^i+\frac{f_1}{f_1+f_2}{N}_{f_1\left(m,n\right)}^{i,j}+\frac{f_1{·f}_2}{f_1^2-f_2^2}{N}_{W\left(m,n\right)}^{i,j}\\ & \end{array}$$

(9)

From Equation (6), it can be found that estimating the SD carrier-phase ambiguity of the reference satellite, $i$, requires a convergence time. When the precision of pseudorange measurement of P codes is about ±0.3 m, the determination of the ambiguity requires about 25 epochs of pseudorange observations [24]. If the epoch number,$l$, is reduced, the determination of SD carrier-phase ambiguity is inevitably affected by the pseudorange noise and multipath. Next, the pseudorange noise and multipath affect the time transfer performance.

In the double-difference positioning solution, it is a common practice to select the satellite with the maximum elevation as the reference satellite. However, when solving the single-difference ambiguity of the reference satellite using Equation (3), it is inevitable that the influence of the pseudorange noise and multipath will be present. To further reduce the impact of the pseudorange noise and multipath on the SD integer ambiguities, we take each visible satellite as a reference satellite separately and estimate the SD integer ambiguities.

As soon as a set of SD integer ambiguity samplings is available, we can classify them into groups according to their distribution. This group is also called a cluster, a region in which the density of samples is locally higher than in other regions [25]. In this contribution, the K-means clustering method is employed to determine SD integer ambiguities. The K-means algorithm aims to choose the best clusters that minimize the sum square error of the distances between a sample and its cluster centroid [25]. The problem can be expressed mathematically as:

$$\left\{\begin{array}{c}P(S\left\{S_1, \dots ,S_n\right\})=min{{\sum }_{i=1}^n\sum }_{j=1}^m {\left(x_{ij}-{\mu }_i\right)}^2\\ {\mu }_i=\frac{1}{m}\sum _{j=1}^m{x}_{ij} i=1, \dots ,n\end{array}\right.$$

(10)

where $S$ is the set of total samples and $\left\{S_1, \dots ,S_n\right\}$ is the set of n clusters from $S_1$ to $S_n$. For each cluster, $S_i$, there are $m$ samples, $x_{ij}$ is the $j$th sample of the $i$th cluster, and ${\mu }_i$ is the centroid of the $i$th cluster.

In this study, the clustering is applied on SD integer ambiguities, and we set the initial groups according to satellite PRN. For each group, its centroids, ${\mu }_i$, from all members, $x_{ij}$, are calculated using Equation (10). In practice, the solving process of K-means is an optimization problem. In this study, we selected the cluster with the most data sampling points as the final SD integer ambiguity.

2.3. Time Difference Estimation of DPT

In order to obtain the time difference solution between two stations, we introduce the SD integer ambiguities into SD carrier-phase measurements, achieving accurate time transfer. For a short baseline of several kilometers, the inter-station receiver clock offset is shown in Equation (11) after ignoring the SD ephemeris error residuals, tropospheric delay residuals, and ionospheric delay residuals.

$${t_r}_{f(m,n)}=\frac{1}{c}\left(L_{f(m,n)}^j-{\rho }_{f(m,n)}^j-{\lambda }_f\cdot N_{f(m,n)}^j\right)$$

(11)

For the IF combination, the receiver clock offset is shown in Equation (12).

$${t_r}_{IF(m,n)}=\frac{1}{c}\left(L_{IF(m,n)}^j-{\rho }_{IF\left(m,n\right)}^j-{{\lambda }_{IF}\cdot N}_{IF\left(m,n\right)}^j{-T}_{IF\left(m,n\right)}^j\right)$$

(12)

where ${t_r}_{f(m,n)}$ represents the $f$ frequency DPT time transfer result and ${t_r}_{IF(m,n)}$ represents the IF combination DPT time transfer result.

3. Experiment Strategy

To explicitly investigate and analyze the performance of the time and frequency transfer solutions with the proposed DPT method, eight baselines composed of IGS stations were selected. The detailed information of the eight baselines is presented in

Table 1. Details of the eight baselines (B and U represent the base and user stations, respectively).

Baseline	Receiver (B)	Antenna (B)	Receiver (U)	Antenna (U)	Distance (km)
USN7-USN8	SEPT POLARX5TR	TPSCR.G5	SEPT POLARX5TR	TPSCR.G5	0
NYA1-NYA2	TRIMBLE NETR8	ASH701073.1	SEPT POLARX5	JAVRING ANT_G5T	0.17
NYA2-NYAL	SEPT POLARX5	JAVRINGANT_G5T	TRIMBLE NETR9	AOAD/M_B	0.16
NYAL-NYA1	TRIMBLE NETR9	AOAD/M_B	TRIMBLE NETR8	ASH701073.1	0.007
KOKV-KOVB	JAVAD TRE_G3TH DELTA	ASH701945G_M	SEPT POLARX5TR	ASH7019 45G_M	0
CEBR-VILL	SEPT POLARX5TR	SEPCHO KE_B3E6	SEPT POLARX5	SEPCHO KE_B3E6	35.30
KOUR-KOUG	SEPT POLARX5TR	LEIAR25.R3	SEPT POLARX5	SEPCHO KE_B3E6	25.07
ONSA-SPT0	SEPT POLARX5TR	AOAD/M_B	SEPT POLARX5TR	TRM59800.00	67.90

. Using the GPS calibration results provided by BIPM as a reference, the time transfer results for the USN7-USN8 baseline verify the accuracy of DPT time transfer. For the baseline closure consisting of the three short baselines NYA1-NYA2, NYA2-NYAL, and NYAL-NYA1, it verifies the self-consistency of DPT time transfer. In evaluating the long-term time and frequency performance of the baselines, four sets of baselines with external atomic clocks, namely, KOKV-KOVB, CEBR-VILL, KOUR-KOUG, and ONSA-SPT0, are used to verify the DPT time and frequency stability.

Table 2. Details of DPT data processing.

Item	Data Processing
Observation types	GPS L1, L2, C1, P2
Sample rate	30 s
Data period	30 days, from Modified Julian Date (MJD) 59,764.0–59,793.0
Tropospheric delay	Saastamoinen (dual frequency)
Ionospheric delay	IGS global ionosphere maps (single frequency) IF combination (dual frequency)
Cut-off angle	10°

summarizes the information about the data and settings in the DPT data processing.

4. Single-Difference Ambiguity Analysis

The float SD carrier-phase ambiguity of the reference satellite can be obtained by using Equation (3). We took each visible satellite as a reference satellite separately and estimated the L1 and L2 SD ambiguities of the reference satellites. For the KOKV-KOVB baseline, the fractional part of the SD ambiguity of the reference satellite is shown in Figure 1. Note that KOKV-KOVB is the zero baseline. For the baseline of ONSA-SPT0 with a length of 67 km, the fractional part of the reference satellite SD ambiguity is shown in Figure 2. It is believed that when the peak-to-peak value of the SD carrier-phase ambiguity fractional part time series is less than 1/3 cycles and it remains stable in subsequent epochs, the integer ambiguity can be determined.

As shown in the upper panel of Figure 1, the L1 SD carrier-phase ambiguities of eight visible satellites were obtained using the L1 carrier phase and the C1 code in MJD 59,764. It can be clearly seen that within 25 epochs (12.5 min for 30 s simple intervals), the peak-to-peak values of the L1 SD carrier-phase ambiguity fractional parts time series are less than 1/3 cycles, and the converged carrier-phase ambiguities remain stable in subsequent epochs. The bottom panel of Figure 1 plots the SD carrier-phase ambiguities obtained using the L2 carrier phase and the P2 code and shows the convergence of the solutions. It took 8 min to reduce the peak-to-peak values of the L2 SD carrier-phase ambiguity fractional parts time series to less than 1/3 cycles. In addition, for the ONSA-SPT0 baseline in Figure 2, it can be observed that the IF SD carrier-phase ambiguities for each visible satellite were resolved within 40 min.

After determining the SD integer ambiguity of the reference satellite and the DD integer ambiguity, we used each visible satellite as a reference satellite to obtain the SD integer ambiguity of the other visible satellites by Equation (6). Figure 3 plots the L1 and L2 SD ambiguities of G07, G08, G10, G13, G18, G24, G27, and G28 for KOKV-KOVB. It is known that the ambiguity of each visible satellite remains continuous as long as the carrier-phase observation is continuously observed. It can be observed that the estimated SD ambiguities differ depending on the reference satellite used. Taking G18 as an example, when using five different satellites as a reference, the estimated L1 SD ambiguities were 6. However, when using the other two different satellites as a reference respectively, the L1 SD ambiguities were estimated to be 7. Similarly, for the L2 SD ambiguities, when using five different satellites as a reference respectively, the estimated values were 63; when using the other one satellite as a reference, the L1 SD ambiguity was estimated to be 62; and using the rest of the satellites as a reference, the L1 SD ambiguity was estimated to be 61. It can be concluded that there are differences in SD ambiguity resolution using different reference satellites, which affects the accuracy of the time difference solution obtained by Equations (10) and (11). Clustering the results of these SD carrier-phase ambiguities, it can be determined that the L1 SD ambiguity is 6 and the L2 SD ambiguity is 63 for G18.

For the ONSA-SPT0 baseline, Figure 4 illustrates the GPS IF combination SD carrier-phase ambiguities of G12, G13, G24, G25, G27, and G29. Taking G13 as an example, it can be observed that the IF SD carrier-phase ambiguities for each visible satellite are different, and it is not the integer number of cycles because IF combined observations cause the ambiguity to lose integer properties. After clustering the results of these SD carrier-phase ambiguities, the final IF SD ambiguity of G13 was determined to be 5.6.

5. DPT Time and Frequency Transfer Validations

In this section, we first verified the accuracy of time transfer through two methods: common clock time transfer and DPT with closure baseline. Then, we evaluated the precision of long-term DPT time transfer through 30 days of continuous solutions. Finally, we analyzed the time and frequency stability of DPT.

5.1. DPT Time Transfer in a Common Clock

Data from two independent receivers, USN7 and USN8, are being operated at the United States Naval Observatory (USNO) in Washington, DC, using a common reference clock that is regularly calibrated using GPS time transfer. As the time transfer solution of the two stations is expected to constant and known, it is relatively easy to evaluate the expected results of DPT.

To verify the accuracy of the DPT time difference, we conducted time transfer experiments on the USN7 and USN8 stations and performed daily re-initialization. We collected 7 days’ CV time transfer results during MJD 59,764.0–59,770.0, as shown in Figure 5. For the estimated time difference over these 7 days, the standard deviation (STD) for the C1 code was 939 ps and for the P2 code it was 88 ps. As expected, it is well known that the precision of the P code is ten times that of the C1 code due to its narrower code chip width. Since the initial value of DPT time transfer depends on pseudorange time transfer, the accuracy of CV time transfer for different codes will affect DPT time transfer.

The DPT time transfer results are shown in Figure 6. Note that the results have been calibrated using the calibration results published by BIPM, as shown in

Table 3. The calibration values for data from the BIPM Time Department’s published GPS calibration results (ns) [26].

Station	TOT (C1)	TOT (P2)	Calibration Uncertainties
USN7	155.6	149.3	1.2
USN8	202.6	197.4	1.2

, by subtracting the time delays associated with the stations CAB, REF, and INT, along with their respective uncertainties [26]. The calculation is as follows:

$$∆t\left(f\right)={t\left(f\right)-\left\{TOT\left(f\right)\right.}_{USN7}-\left.{TOT\left(f\right)}_{USN8}\right\}\to 0$$

(13)

where f represents the L1 frequency or L2 frequency, $t(f)$ represents the DPT solution, and $TOT\left(f\right)$ represents the total time delay of a geodetic system for frequency f, which is divided into six different parts, including antenna delay, antenna cable delay, receiver internal delay, short cable and splitter delay, internal reference delay, and 1 PPS delay [27].

For the time difference solutions estimated for 7 days in Figure 6, the STD and mean results of the GPS DPT and CV time transfer daily solutions after convergence are shown in

Table 4. STD and mean results of GPS DPT and CV time transfer daily solutions for USN7-USN8 (ps).

MJD	CV Time Transfer				DPT Time Transfer
	Mean_C1	Mean_P2	STD_C1	STD_P2	Mean_L1	Mean_L2	STD_L1	STD_L2
59,764	292	64	1009	99	257	65	12	4
59,765	349	64	952	97	255	65	10	4
59,766	317	64	875	96	254	65	10	3
59,767	322	64	996	87	327	64	14	4
59,768	264	62	976	89	323	64	12	4
59,769	304	61	890	86	324	65	12	5
59,770	321	63	899	91	321	64	12	4

. The STDs of the converged time difference for the GPS L1 and L2 frequencies were 18 ps and 4 ps, respectively. As expected, due to the heavy noise level in the C1 code, time transfer relying on the C1 code requires a long time to converge, making it difficult to determine the absolute time for DPT with high precision. Meanwhile, the estimated time difference solutions for the 7 days had mean convergence values of 291 ps for GPS L1 frequency and 70 ps for the GPS L2 frequency. As a time transfer experiment, the time difference solution between the USN7 and USN8 stations should be zero in theory, and the closure of the two stations shows the accuracy of the DPT time transfer. In Figure 6, there is a constant difference of about 300 ps using the C1 code. This difference may be due to uncertainties of DPT solution or the calibration uncertainties from BIPM, or the combination of DPT and the calibration. Furthermore, we validate the time transfer of DPT through closure experiments in the next sub-section.

5.2. DPT Time Transfer Validations with Closure Baselines

In order to further verify the accuracy of the DPT time difference, we determined the time transfer between each of the NYA1, NYA2, and NYAL stations and performed daily re-initialization for 7 days. We further evaluated the self-consistency of the DPT using the three stations’ closure in Figure 7. Theoretically, the three stations’ baseline closure is given by:

$$\mathit{Closure}=(NYA1-NYA2)+(NYA2-NYAL)+(NYAL-NYA1)=0$$

(14)

As shown in Figure 7, the estimated closure solutions for the 7 days had STD values of 4 ps for the GPS L1 frequency and 3 ps for the GPS L2 frequency after the initialization. Meanwhile, Figure 8 shows the mean value of the baseline closure of the daily solutions. The two dashed lines in Figure 8 represent the average values over 7 days for GPS L1 and L2 frequencies, with the GPS L1 frequency at −14 ps and the GPS L2 frequency at 9 ps after initialization. Notably, the closure did not change over time, indicating that the DPT processing is self-consistent in terms of time transfer. This consistency reinforces the reliability of the DPT technique for accurate time transfer.

5.3. Long-Term DPT Time and Frequency Transfer Validations

We analyzed the common clock experiment between USN7 and USN8 and the three stations’ baseline closure test and verified that DPT can accurately determine the initial value of DPT time transfer. In addition, we also selected four baselines, CEBR-VILL, ONSA-SPT0, KOUR-KOUG, and KOKV-KOVB, for continuous time transfer and further considered the long-term performance of DPT. Table 1 presents the details of the four baselines. To remove the trend effects of the time difference between two station clocks, a common reference time is introduced to subtract the inconsistency between the two clocks. Typically, IGS final station clock solutions are used as the reference times due to their precision of 20 ps (STD) [28]. The data sampling rate of the IGS final station clock solutions is 5 min. In order to better compare with the reference, the time difference solution of DPT is also set with a sampling interval of 300 s. The result of DPT is then:

$$Time transfer=(T_A-T_B)-({IGST}_A-{IGST}_B)$$

(15)

Figure 9 displays the time differences between different L1 and L2 frequencies and the dual-frequency IF combination DPT for four different baselines. It should be noted that none of the four baselines were calibrated and that a constant offset was added to each curve for comparison. It can be seen that the time difference estimated by the five curves shows the same trend. Among them, the KOKV-KOKB baseline was processed using single-frequency L1 or L2 time transfer, with STDs of 18 ps (the L1 frequency for KOKV-KOKB (L1)) and 17 ps (the L2 frequency for KOKV-KOKB (L2)) over a period of 30 days. The reason why the time difference results of the DPT time transfer are so flat is that the initial value of DPT time transfer depends on pseudorange time transfer, and after this initial value is determined, only the carrier phase is used, so the noise in the code is well eliminated. As the baseline distance between the reference station and the user station increases, the DD atmospheric delay cannot be ignored. Therefore, the CEBR-VILL, KOUR-KOUG, and ONSA-SPT0 baselines were processed using the dual-frequency IF combination, with STDs of 69, 75, and 94 ps, respectively. Compared with single-frequency processing, the IF combination used here amplifies the noise to a certain extent, especially for the ONSA-SPT0 baseline shown in Figure 9.

From the above, we analyzed the long-term time transfer performance of DPT of four baselines, CEBR-VILL, ONSA-SPT0, KOUR-KOUG, and KOKV-KOVB. For the evaluation of frequency stability, compared with Allan deviation, modified Allan deviation (MDEV) has the advantage of distinguishing white and flickering phase modulation noise and is often used for frequency stability evaluation [29]. In addition, time deviation (TDEV) is used to characterize the time stability of time links [30]. In this work, we used TDEV and MDEV to analyze the time and frequency stability of DPT, respectively, which were obtained using the Stable32 software [31]. Smaller MDEV and TDEV values indicate better time and frequency stability.

Figure 10 and Figure 11 show the stability results of MDEV and TDEV for the four baselines. As shown in Figure 10, the single frequency DPT has an MDEV of 9 × 10⁻¹⁷ at an averaging time of 1 day and 3 × 10⁻¹⁸ at an averaging time of 7 days. The IF DPT has an MDEV of 3 × 10⁻¹⁶ at an averaging time of 1 day and 2 × 10⁻¹⁷ at an averaging time of 7 days. Figure 11 shows that the time stability of the IF DPT is less than 20 ps for an averaging time of 1 day, especially less than 10 ps within 7 days. As for the single-frequency DPT, the time stability is less than 4 ps for an averaging time of 1 day and approaches 1 ps within 7 days.

6. Discussion

Compared with some existing results of GNSS pseudorange time and frequency transfer techniques, such as CV and AV, DPT indeed improves the time and frequency stability and thus enhances the performance of frequency comparison over averaging times of several hours or longer. In addition, a K-means clustering method was employed to estimate the SD receiver clock offset, which achieves accurate time transfer. The time stability can reach the picosecond level on the averaging time of 7 days, while the frequency stability can reach the sub-10⁻¹⁷ level. Among the existing GNSS time and frequency transfer methods, integer PPP (IPPP) has been proven to have the best clock comparison performance [32,33], and the common clock IPPP approach has achieved a time transfer precision of 61 ps, which is comparatively less precise than DPT time transfer in the KOKB-KOKV baseline. In terms of frequency stability, the MDEV of the common clock IPPP experiment over a 30-day period is no less than 1 × 10⁻¹⁷, while the MDEV of the DPT frequency transfer for the KOKB-KOKV baseline can reach as low as 3 × 10⁻¹⁸ within 30 days. To our knowledge, these results represent the optimal proven time transfer precision and frequency stability performance of GNSS time and frequency transfer.

As expected, for short baselines, single-frequency DPT provides higher precision than dual-frequency IF DPT, and the dual-frequency IF combination amplifies noise to some extent. Therefore, for short baselines, we can prioritize the use of single-frequency DPT for time and frequency comparison. The use of the DPT method is complementary to the PPP time transfer method and can achieve high-precision regional time transfer. It only uses broadcast ephemeris to achieve high-precision time transfer and is easier to implement in real-time dynamic applications. In addition, the excellent frequency stability of DPT frequency transfer helps to evaluate the frequency comparison performance in the region of the latest-generation optical lattice clock frequency standard. It should be noted that this work did not discuss longer baselines. For longer baselines, we will consider incorporating precise ephemeris in future research to improve time transfer accuracy.

7. Conclusions

In this paper, a differential precise time transfer (DPT) with a clustering constraint method was employed to estimate the time difference between two clocks, aiming to achieve accurate and precise time and frequency transfer. For the experiment of common clock time transfer, the DPT time transfer results showed good consistency across different days, with an intra-day accuracy of within 20 ps. Furthermore, we evaluated the self-consistency of the DPT technique using closure among three stations. The baseline DPT closure values of the three stations were always within 10 ps, with a peak-to-peak value of closure of about 25 ps. The closure did not change over time, indicating the self-consistency of the DPT processing in time transfer. Using the station clock solutions provided by IGS as a reference, the STDs of the four baselines were all less than 100 ps within 30 days. It was found that the single-frequency DPT achieved a 10⁻¹⁷ level of MDEV at an averaging time of about 1 day and a 10⁻¹⁸ level at an averaging time of 7 days. The IF DPT achieved a 10⁻¹⁶ level of MDEV at an averaging time of about 1 day and of 10⁻¹⁷ at an averaging time of 7 days. The TDEV of the IF DPT was less than 20 ps for an averaging time of 1 day, especially less than 10 ps within 7 days. As for the single-frequency DPT, the time stability was less than 4 ps for an averaging time of 1 day and approached 1 ps within 7 days. Finally, it can be used for high-accuracy and high-precision time and frequency transfer applications.