GPS Time and Frequency Transfer : PPP and Phase-Only Analysis

To compute precise point positioning (PPP) and precise time transfer using GPS code and phase measurements, a new software named Atomium was developed by the Royal Observatory of Belgium. Atomium was also adapted to perform a phase-only analysis with the goal to obtain a continuous clock solution which is independent of the GPS codes. In this paper, the analysis strategy used in Atomium is described and the clock solutions obtained through the phase-only approach are compared to the results from the PPP mode. It is shown that the phase-only solution improves the stability of the time link for averaging times smaller than 7 days and that the phase-only solution is very sensitive to the station coordinates used. The method is, however, shown to perform better than the IGS clock solution in case of changes in the GPS receiver hardware delays which affects the code measurements.


INTRODUCTION
GPS-based geodetic time and frequency transfer are based on a consistent modeling of both code and carrier phase measurements, and is presently widely recognized for its high precision of about 100 picoseconds at each epoch [1][2][3][4][5].Its major limitation concerns the discontinuities at the day boundaries, which are caused by the noise of the code measurements.Indeed, the carrier phase measurements contain an unknown initial ambiguity (integer number of cycles) and can, therefore, only provide information on the clock evolution or frequency transfer.The absolute synchronization error between two clocks is obtained from the code measurements, so that its accuracy is limited by the code noise (a few nanoseconds).As GPS data are classically analyzed by daily batches, the absolute clock synchronization error is determined for each day from the code measurements gathered during that day.Due to the colored signature of the code noise, implying large and long-term temporal variations in all code measurements, the clock solutions show discontinuities at the day boundaries.The origin of this colored noise in the codes is presently not yet fully understood [6], but it reduces our ability to access the true clock signal and stresses the necessity to develop a rigorous approach for continuous geodetic time transfer.One solution was already proposed in [7]; it consists in performing the PPP analysis on a longer data batch, reducing consequently the number of day-boundary jumps.A second approach was proposed in [8], where two methods were tested: a clock handover, analyzing daily data batches with one common observation epoch between two consecutive days (midnight), or a stacking of the carrier phase ambiguities, through a stacking of normal equations for consecutive days.
In order to have at our disposal a straightforward tool for the rapid computation of clock solutions in a PPP mode, the ROB developed Atomium.Atomium performs a least squares analysis of single station GPS code and carrier phase data to provide for each observation epoch a clock receiver solution consisting of the synchronization error between the receiver clock and the IGS time scale (IGST).Recently, Atomium has been adapted to provide also a single station phase-only solution.Some preliminary results were presented in [9], and a more comprehensive study of the approach and of its sensitivity to the input parameters is presented here.
The basic observable in the phase-only solution is the difference between successive carrier phase observations from the same satellite.Using a least squares adjustment, International Journal of Navigation and Observation this observable is used to obtain the evolution of the clock synchronization error.The advantage of the phase-only approach is that it provides a continuous solution across the day boundaries, but the drawback is that the absolute clock synchronization errors of course cannot be determined due to the presence of the carrier phase ambiguities.This means that the continuity of the clock solution is interrupted after each tracking interruption.
For the computation of TAI [10] or any other time scale, each period of continuous clock solutions can be combined with a calibrated time transfer method such as two way satellite time and frequency transfer (TWSTFT) which allows synchronizing remote clocks with an accuracy better than 100 picoseconds and has a reliable long-term stability.However, TWSTFT measurements are disturbed by a diurnal variation of 1-3 nanoseconds peak to peak and their resolution is poor (one point each 2 hours in the best case).On the other hand, the GPS phase-only analysis has no calibration, but a very high resolution (5 minutes) and good short-term frequency stability.Combining the two methods will take advantage of the strong points of each of them: the accuracy for TWSTFT and the stability for GPS carrier-phase [10].
The first part of this paper presents the analysis procedure used by Atomium to produce station position and clocks using PPP, and shows the clock solutions obtained for several stations in comparison with the IGS clock solution (IGS combined solution [11]) as well as the solutions obtained by other software packages (NRCan and Bernese 5.0).The second part describes the theoretical background of the phase-only approach implemented in Atomium, compares the phase-only with the PPP results, and investigates the sensitivity of the phase-only method to its input parameters.

PPP APPROACH
Precise point positioning (PPP) consists in determining the station position, receiver clock, and tropospheric zenith path delay without using any measurements from other stations.This approach requires the availability of precise satellite orbits and clocks provided by an external source.The PPP procedure is fully described in [12]; we just recall here the main principle and its implementation in Atomium.
The observation equations for carrier phases (L 1 and L 2 ) and pseudoranges (P 1 and P 2 ) are with R the geometric distance receiver-satellite, τ s the satellite clock error, τ r the receiver clock error, τ t the tropospheric delay, τ i the ionospheric delay, λ the carrier wavelength, N the phase ambiguity, τ d the instrumental code delay, and n the noise.
Atomium is based on the ionosphere-free combinations of L 1 and L 2 and of P 1 and P 2 , named L 3 and P 3 , respectively.The observations are used at the 5-minute sampling rate.The satellite positions are obtained using a Neville interpolation on 12 points of the IGS sp3 files in which the satellite positions are given at a 15-minute sampling rate.The satellite clock corrections are extracted from the IGS clock files in which the sampling rate is 5 minutes.The station position is corrected for its time variations due to degree 2 solid Earth tides as given in the IERS conventions [13] and ocean loading effects, using the FES2004 model [14].The respective elevation (no-azimuth) and nadir-dependent absolute corrections for the receiver and satellite antenna phase center variations as made available by the IGS (ftp://igscb.jpl.nasa.gov/pub/station/general/igs05.atx) are applied, and the carrier phase measurements are corrected for phase windup taking into account the satellite attitude.Periods of eclipse events are eliminated.The instrumental code delays are considered as constant and not included in the present version of the software.
The tropospheric delay is expressed as the sum of the hydrostatic and wet delays, both being the product between a given mapping function (mf) and the zenith path delay (zpd).The hydrostatic part is introduced using the Saastamoinen a priori model [15] with the dry Niell mapping function [16] and the wet part uses the Niell wet mapping function.The wet zpd is estimated with a 2hour sampling rate, with linear interpolation between the epochs of consecutively estimated zpd's.Using a least squares scheme, with a weighting for the codes and carrier phases associated with the noise level of each observation type and satellite elevation, Atomium provides an estimation of (i) the receiver clock delay, at each epoch (5-minutes sampling rate), with respect to IGST; (ii) the position for the whole day (optional: the position can be fixed to a given value); (iii) the zpd, with a sampling rate of 2 hours; (iv) the phase ambiguities.
Figure 1 presents the clock solution obtained with Atomium in a PPP mode for three IGS stations over a period of 15 days.Two of them (NRC1 and BRUS) are equipped with an H-Maser, while the station PTBB is equipped with a cesium clock.NRC1 and BRUS have been chosen for their, respectively, characteristic large and small day boundary jumps.Figure 1 also shows the differences between the Atomium clock solution and the combined IGS clock solution, and this for the period January-Febuary 2007.The differences between the Atomium clock results and the IGS clock solution for all the stations analyzed (not all are shown here) show a constant bias ranging between 0 and 200 picoseconds plus some small variations with an rms between 50 and 100 picoseconds.
In order to validate the Atomium software, the clock results obtained from the PPP analysis were also compared with the results obtained from two other software capable of performing PPP: NRCan [12] and Bernese 5.0 [17].The Bernese software does not directly interpolate the IGS orbits, as done by the NRCan and Atomium software, but it transforms the SP3 positions to an inertial frame, performs a numerical integration of the equations of motion and then converts the obtained satellite positions back to the earthfixed system.
Pascale Defraigne et al.   Figure 2 shows the comparison between the results obtained for IENG over one week with Atomium, NRCan, and Bernese.The NRCan solution has been corrected for a constant bias of about 4 nanoseconds with respect to the other solutions.This bias is caused by the fact that the version of the NRCan software used in this work is using relative corrections for the satellite and antenna phase center variations while the Bernese, Atomium, and the IGS use absolute corrections since November 2006.The NRCan bias remains visible in the lower part of Figure 2 which displays the differences between the different clock results and the IGS clock solution.We can see that the differences have similar amplitudes with an rms ranging between 30 and 100 picoseconds over the analyzed week, but note that the rms is affected both by the accuracy of the solution during the different days, and its precision.

Description of the method
As explained above, the use of PPP for time transfer applications is mainly limited by the presence of day-boundary jumps in the clock solution.These day-boundary jumps are station-dependent and can be large for some stations, see Figure 1 for NRC1, or [18].One solution to suppress the day-boundary jumps is to exploit the continuity of the phase measurements across midnight and to provide a continuous solution without using any pseudoranges.
The phase-only approach developed here is based on this principle and has the advantage of being simple and rapid to use.It consists of using the phase differences between two consecutive epochs (using a 5-minute sampling interval) as basic observable, to create the normal equations accordingly, and to solve for the clock differences between the consecutive epochs.
The first derivative of ( 1) is used to generate the phase derivatives of the ionosphere-free carrier phase measurements: so that where dτ s , dτ r , and dτ t are the differences between the satellite clock errors, the receiver clock synchronization errors, and the tropospheric delays, at times (t i ) and (t i−1 ).This new quantity does not contain any ambiguity as long as there is no tracking interruption.In the occurrence of a cycle slip, the corresponding phase difference will be rejected as an outlier.The phase-only normal equations are constructed and solved in a similar way as for PPP.However, using as observations differences between successive carrier phase measurements, the position can only be determined with a limited precision similar to the classic triple difference phase processing.For this reason, we choose in this stage to fix this parameter to a given value.The sensitivity of the results to the station position (fixed or free) will be discussed in Section 3.5.The phase-only version of Atomium therefore provides (i) the receiver clock derivative dτ r at each epoch; (ii) the zpd, with a sampling rate of 2 hours.
The phase-only clock solution is then obtained from integrating the receiver clock derivative at each epoch obtained from the least squares estimation.
Similar to the strategy used in the PPP mode, the phaseonly analysis is also done in daily batches.But, now the successive days are linked by computing the difference between the first observation of the day and the last observations of the previous day exactly as the differences between two successive observations within a day.This guarantees the continuity of the clock solution across the day boundaries.In case of an actual tracking interruption, the phase-only solution is interrupted and the integration must be restarted.

Discontinuities in the IGS products
While the phase-only analysis has the potential to provide a continuous solution, still some jumps with amplitudes up to 120 picoseconds can appear at the day boundaries.These jumps are caused by the nonperfect overlapping of the IGS orbit and clock files at midnight.The jumps are geographically correlated as their amplitudes depend on the visible satellites at midnight GPS time.This is illustrated in Figure 3, where the stations in Europe (BRUS, IENG) show small jumps, while stations in North-America (USN3, USNO, NRC1, NISU) show a large jump (above 100 picoseconds).
The day-boundary jumps in the phase-only solution are reduced in Atomium by eliminating the satellites with a discontinuity in the IGS orbit and clock information using a 2-sigma threshold.This is, however, only possible when the receiver clock is an H-maser as for less stable clocks the dayboundary phase jumps are lost within the clock noise.The impact of the jumps in that case is, however, partly reduced by a proper weighting of the satellite orbits accounting for the orbit accuracy deduced from the 7-day arc evaluations and given in the IGS reports for IGS final orbits.

Comparison between the phase-only, the PPP, and the IGS receiver clock solutions
The continuous clock solution obtained from the phase-only analysis was compared over one month with the IGS and PPP clock solutions.Figure 4 presents the three solutions for the IGS station NRC1 which is, during winter time, characterized by large day-boundary discontinuities.On a daily basis, the phase-only solution and the PPP solution are equivalent, but the day-boundary jumps are absent from this phase-only solution and the improvement brought by this new method is obvious.Thanks to the elimination of the day-boundary jumps, our phase-only solution improves the stability of the time link for averaging times up to 7 days, as shown in Figure 5 for NRC1-IGST.This improvement confirms the results obtained in [8]. the improvement depends on the size of the day-boundary jumps in the PPP (or equivalently IGS) solution.

Comparison with the NRCan multiday batch analysis
We also compared our phase-only analysis with the continuous (or multiday) PPP solution obtained with the NRCan software [7].This software is able to analyze a data batch of several days with sequential least squares, using both code and phase observables.The resulting multiday PPP solution does not suffer from day-boundary discontinuities and one single position can be estimated for the whole period processed.
Figure 6 shows the comparison between the Atomium phase-only results and the NRCan 30-day PPP solution for the stations NRC1, USN3, and PTBB.For NRC1 and USN3, the Atomium solution shows a drift with respect to the NRCan and IGS solution; this drift is due to the position fixed within the Atomium phase-only analysis; this problem will be discussed in more detail in Section 3.5.For PTBB, a large jump occurs in the IGS solution between mjd 54144 and 54146; this jump is caused by a jump in the code data at mjd 54145.4 and was later confirmed as a change of the receiver hardware delay between 2006 and 2007 (G.Petit, personal communication).There is no IGS clock solution for PTBB for the day with the jump as the IGS analysis rejected the solution because of this problem.As our phase-only approach only uses the carrier phases and neglects the code measurements, it is insensitive to the code jump.However, in the multiday NRCan solution, the jump is distributed over about 10 days inducing an artificial drift of the solution during these days.Also the IGS solution is degraded in this case.This is confirmed by the comparison with TWSTFT measurements between PTBB and another station (USN3 in our case, see Figure 7): neither a discontinuity nor a frequency drift is present in the TWSTFT solution and only  the Atomium phase-only approach is in agreement with the results of this independent technique.Note that in that case, the PTBB and USN3 coordinates have been chosen in order to have no drift between the IGS solution and the phase-only solution (see Section 3.5 for the sensitivity of the results to station position).

Sensitivity to the station position
All the results presented up to now have been obtained by fixing the station position in the analysis.As explained before, the station coordinates can be determined using phase differences between successive epochs, but with a International Journal of Navigation and Observation  limited precision due to the strong correlation between the estimated position and the clock solution.To quantify the sensitivity of the phase-only analysis to the station coordinates, a test was performed using the IGS station BRUS which is characterized by small day-boundary discontinuities.The phase-only solution has been computed over one month either by determining the station position for each day, or by fixing the station position to a given value.In the case of the fixed station position, two different positions separated by 3 mm in the Y-component are tested.The results are presented in Figure 8: the three curves deviate with respect to the IGS solution.Furthermore, we can observe a significant drift between the two solutions obtained with fixed positions, leading to a disagreement of 1.5 nanoseconds after only one month.The drift between the IGS solution and the phase-only solution is, therefore, most probably due to the imperfect station coordinates introduced or determined during the analysis and possible small modeling differences with respect to the IGS solution.The choice of the station coordinates with a subcentimeter precision is, therefore, clearly a crucial point in the phase-only analysis.

CONCLUSIONS
A new tool for precise GPS time and frequency transfer has been presented.It is called Atomium and is based on a least-squares analysis of GPS measurements.A PPP solution can be obtained using both codes and carrier phases, and a phase-only solution is proposed in order to reduce the dayboundary jumps or the intraday jumps which are due to the usage of the code measurements.The accuracy of the results obtained in the PPP mode has been determined from the comparison of these results with the IGS combined clock products.The bias between the Atomium clock solutions and the corresponding IGS clock products is station-dependent, while always remaining below 200 picoseconds.In addition, some small variations, with an rms between 50 and 100 picoseconds, are present.They are similar to the variations between two other software performing PPP (Bernese 5.0 and NRCan) and the IGS clock solution.
A continuous frequency transfer method, based on a phase-only analysis, is also presented in this paper.As demonstrated, the advantage of the phase-only analysis is that it can be combined with a calibrated time transfer method as two way satellite time and frequency transfer (TWSTFT) taking advantage of the strong points of each of the methods: the accuracy for TWSTFT and the stability and high resolution for GPS carrier-phase [9].
Our continuous clock solutions are based on the least squares analysis of the differences between the consecutive GPS phase measurements.The analysis is performed in daily data batches, but the successive days are linked by computing the differences between the first observations of the day and the last observations of the previous day.The results show that a satellite-dependent weighting is necessary to assure that the phase-only solution does not show dayboundary discontinuities (up to 120 picoseconds) due to the nonperfect overlapping of the IGS clock files at midnight.We have demonstrated that the phase-only clock solution improves the stability of the time link for averaging times up to 7 days; the amplitude of the improvement depends of course on the amplitude of the day-boundary discontinuities in the PPP solution.The high sensitivity of the phase-only clock solution to the station coordinates fixed in the analysis was also demonstrated, indicating that a subcentimeter precision is needed for the position.Finally, the efficiency of the phase-only solution in case of changes in the GPS receiver hardware delays was shown; contrarily to the results obtained using both GPS code and phase observations, the phase-only solution is not affected by these hardware delay variations.
clock solution IGS clock solution Difference between Atomium and the IGS solution over 2 month (f)

Figure 1 :
Figure 1: Comparison between Atomium PPP results and the IGS combined clock solutions for BRUS (a,b), NRC1 (c,d) and PTBB (e,f).

Figure 2 :
Figure 2: Comparison between the results obtained for IENG over one week with Atomium, NRCan, Bernese 5.0, and the IGS combined clock solution.

Figure 3 :Figure 4 :
Figure 3: Jumps in the phase-only solution due to discontinuities in the daily satellite orbit/clock products.

Figure 5 :
Figure 5: Modified Allan deviation of the phase-only solution and of the daily-independent PPP solution for NRC1-IGST.

Figure 6 :Figure 7 :
Figure 6: Comparison between the phase-only solution and the NRCan solution obtained through the analysis of a 30-day data batch.

Figure 8 :
Figure 8: Phase-only clock solution for BRUS obtained either by determining the station position for each day or by fixing the station position to a given value, using two different positions separated by an offset of 3 millimeters in the Y-component.