High-Rate Monitoring of Satellite Clocks Using Two Methods of Averaging Time

: Knowledge of the global navigation satellite system (GNSS) satellite clock error is crucial in real-time precise point positioning (PPP), seismology, and many other high-rate GNSS applications. In this work, the authors show the characterisation of the atomic GNSS clock’s stability and its dependency on the adopted orbit type using Allan deviation with two methods of averaging time. Four International GNSS Service (IGS) orbit types were used: broadcast, ultra-rapid, rapid and final orbit. The calculations were made using high-rate 1 Hz observations from the IGS stations equipped with external clocks (oscillators). The most stable receiver oscillator was chosen as a reference clock. The results show the advantage of the newest GPS satellite block with respect to the other satellites. Significant differences in the results based on the orbit type used have not been recorded. Many averaging time methods used in Allan deviation (ADEV) show the clock’s fluctuations, usually smoothed in 2 n s averaging times. -10 s for broadcast orbit and between 4.0×10 -12 s to 1.9×10 -10 s in case of ultra-rapid one. Analysis shows a similar behavior for a IIR and IIR-M GPS satellites, while IIF satellites (G26 and G27) are characterized by lower stability for the shortest intervals (1 s and 2 s), clearly higher for longer intervals (4–128 s), while for the highest intervals the stability of TDEV is the same. The same can be observed in the case of uncertainty estimates, error bars of results acquired by using broadcast orbit are much bigger than using ultra rapid ones. Four different orbit types were analysed: three precise with 15 min sampling intervals, and broadcast ephemeris with a 30 min sampling interval. The results from the data indicate GPS satellites frequency stability (ADEV) at a level of 1.7 × 10 − 11 s to 7.2 × 10 − 14 s and time stability (TDEV) at a level of 2.9 × 10 10 s to 3.9 × 10 12 s. The results also show the insignificant advantage of the satellite Block IIF’s frequency and time stabilities in regard to the other two blocks. As used in this article, the many τ method allows for detection of the satellite clock anomaly. This method might also assist in preventing the use of a satellite clock with identified potential instability patterns as a reference clock.


Introduction
GNSS (Global Navigation Satellite System) are widely used for navigation, positioning, scientific [1][2][3][4], and engineering applications [5][6][7]. That includes deformations [8][9][10] and structure monitoring [11][12][13], GIS (Geographic Information System) applications [14][15][16], land use [16][17][18][19], and cadastre [20][21][22][23][24]. Among the widely known positioning techniques, such as static measurements or RTK (Real Time Kinematic), PPP (Precise Point Positioning) is increasingly growing in use; it is a relatively new positioning technique, presented for the first time in [25]. Based on dual-frequency, absolute carrier phase measurements, and high accuracy satellite orbits and clock corrections provided by external services such as CODE (Center for Orbit Determination in Europe) or IGS (International GNSS Service) [26], PPP provides cm-level accuracy. Consequently, high accuracy clock products are the basis for PPP operation. GNSS clocks are also used in a large number of other fields, including mobile phones, computers, and radio transceivers [27]. Simple calculations involving the mean altitude of GNSS satellites as 20,000 km and speed of light as 3 × 10 5 km/s leads to a 30 m range error in the case of 0.1 μs timing error [28]. Due to the construction of GNSS space segments and a motionless receiver, its clock correction at 1 μs level leads to < 1 mm geometrical distance error [29,30].
Each GNSS satellite is equipped with more than one atomic clock [31] which has a direct impact on the navigation positioning accuracy. Any errors in the on-board time result in errors in positioning. Therefore, monitoring of the satellite on-board clocks is very important; in other words the accuracy of the GNSS system depends on the nature of the satellite clocks [32,33]. The stability of the on-board clocks is essential in the GNSS positioning, and the estimation of these clocks requires a large number of tracking ground stations [34]. GNSS on-board clocks have been widely studied to date [35][36][37][38][39], including their influence on the PPP technique [40]. Five minute interval precise clock products are enough for an mm-dm level of PPP precision. For kinematic PPP, a precise satellite clock product with a 30 s or 5 s interval (compared to a 5 min interval) could improve positioning accuracy by up to 30% or more [41].
Calculation of the clock corrections requires at least one reference clock as research shows a more stable oscillation reference clock leads to improvements in positioning accuracy [42]. Satellite clock corrections are available via a few sources/agencies. The first clock corrections are transmitted in real time in navigation (broadcast ephemeris) message. Broadcast ephemeris is forecasted, predicted and interpolated based on the satellite orbits and it does not have enough accuracy for precise applications. The second clock corrections are computed and published in precise SP3 (Standard Product #3 format) ephemeris files together with the satellite's coordinates. Precise orbits are almost two orders of magnitude more precise than broadcast ones in the case of position, and almost two times more accurate in the case of clocks [43]. There are a couple of precise ephemeris models, depending on their accuracy and availability. The fastest available ultra-rapid orbits are accessible in near real-time [44]. The most accurate are rapid and final orbits, but these are only available after 17-41 h and 12-18 days, respectively (Table 1).

Materials and Methods
Each GNSS satellite is equipped with atomic clocks; these clocks are very stable but, due to lack of synchronisation with control segment model time, are subject to drift. In the time frequency domain, satellite clock stability is described by three factors-phase ( ), frequency ( ), and frequency drift ( )-expressed as a quadratic polynomial clock offset model [46]. The GPS broadcast message contains timing elements, and the clock corrections (factors) parameters are clock bias ( ), clock drift ( ), and drift rate ( ). According to [47], the time correction of the particular SV (space vehicle) using [48] naming conventions to the nominal SV time is defined as: Polynomial coefficients ( , , and ) are transmitted in sec, sec/sec, and sec/sec 2 , respectively; -known as a clock data reference time-is also transmitted in a broadcast message in seconds. An external oscillator is mounted only on a minority of IGS (155 out of 504, http://www.igs.org/network, access 24 July 2018) and EPN (European Permanent GNSS Network) stations (40 out of 324, http://epncb.oma.be/_networkdata/stationlist.php, access 24 July 2018). Analysis of the atomic oscillator's characteristics depends on its deviations to nominal frequency. Using two adjacent clock observations ( , ), divided by a nominal time interval ( ), an average fractional deviation ( ) might be computed as [49]: Classical standard deviation does not guarantee convergence for all power-law processes observed in high stability oscillators, therefore in such elaborations Allan variance (AVAR), with square root Allan deviation (ADEV), is used. In atomic clock performance, frequency stability is one of the most important indicators and is usually expressed by the Allan variance [50,51]. The AVAR can be expressed as: where is the average fractional frequency deviation of the interval : Combining Equations (3) and (4) leads to the equation: The time Allan variance (TVAR), with square root (TDEV), is a measure of time stability based on the modified Allan variance and defined as: In simple terms TDEV is MDEV (Modified ADEV), whose slope on a log-log plot is transposed by +1 and normalised by √3 [52], thus TDEV is focused more on time rather than frequency measurements [53].
Plotting the Allan standard deviation is usually observations/data versus time interval τ in a log-log graph. It is also referred to as a two-sample variance, which describes the analysed clock behaviour as a function of time with τ-size samples [38].
The satellite clock stability analysis was based on all of the available types of orbits produced by IGS's CODE center: broadcast, ultra-rapid, rapid, and final (Table 1). Modelling was performed using Bernese GNSS Software 5.2 (BSW). For the calculations 1 Hz observation data for the period 1 January 2018 00:00-02:00 UTC (Universal Time Coordinated, 7200 epochs), registered by the stations as shown in Error! Reference source not found. Figure 1, were used. Earth's orientation parameters, precise orbits and clocks, global ionosphere maps, and DCB (Differential Code Biases) files were CODE's products. Observation data with significant gaps were automatically removed, the reference clock was chosen manually to minimalize its influence (Section 3 and Figure 2). Coordinates (based on the final weekly solutions) and troposphere estimates were adopted as fixed. The final solution was combination of phase and code measurements, the code observations for individual station was subject of station specific weighting. In the final solution phase-only observations were adopted on the basis of the back-substitution step of clock parameters [30]. IGS stations located on European territory operating at a 1 second (1 Hz) sampling rate were included for the experiment's tracking network. Only four of them were equipped with the external hydrogen maser frequency standard: ONS1, MATE, WTZR, and YEBE (Table 3).

Results
The most stable station, equipped with a hydrogen maser oscillator, was chosen as a reference clock. Based on 30 s IGS clock products, the three remaining stations' ADEV are presented inFigure 2. Due to the lack of IGS products for YEBE, this station was removed from the reference station processing. The WTZR graph is clearly outstanding for comparing two of the IGS stations (MATE, ONS1); the stability of this oscillator varies between 9.3 × 10 -14 s and 1.7 × 10 -11 s. The other two graphs are very similar for short τ (averaging time) between 30 s and 240 s. However, for other averaging times, the ONS1 oscillator is much more stable due to the fact that it was chosen as the reference clock for further calculations.
For a clock estimation, code and phase data were used together. Observations from all of introduced stations were taken into consideration to extend the number of visible satellites. Firstly, code measurements were smoothed to check the consistency between code and phase observations. Next, for further analysis smoothed code and original phase observations were taken into consideration. Then clock parameters were estimated using the epoch by epoch method described in more detail by [38,54]. During the calculations, 30 s IGS clocks corrections were introduced as the reference for the consecutive solutions in combination with the manually adopted reference clock selected in Section 4. As the result, 1 s clock corrections were calculated on the basis of the four different orbit types. Different orbit types impact directly on the range and result in the differences in clock error estimation [55]. ADEV was calculated based on the computed 1 s clock corrections.   Figure 3 shows ADEV based on calculations according to two different types of orbit: broadcast and ultra-rapid. The results between three precise orbits were negligibly small, therefore the authors in the following part of this paper focused only on the navigation orbit (broadcast) and ultra-rapid (precise orbit available in real time). The ONS1 station clock was applied as a reference one for each solution. In each part of Figure 3 the most stable are the pair of satellites 26 and 27 (SVN 66 and 71), which highlights the advantage of Block IIF's (red lines) stability in comparison to Blocks IIR (blue) and IIR-M (green). For a G26 satellite, the error bars are calculated that represent standard deviation with a 95% confidence interval. Depending on the averaging time and orbit type adopted for modelling error, the magnitudes are different. In general, a longer averaging time is affected by a larger error due to the smaller number of samples. However, a comparison of broadcast and ultrarapid orbit dependent results shows a smaller error for a precise one orbit for each analysed averaging time. Moreover, there are no error bars on the figures presenting many τ method due to relatively much greater number of results using the 2 n method (Figures 3 and 4 Figure 4. ADEV depending on orbit type using the many τ method.

Allan Deviation: τ = 2 n
The stabilities of Blocks IIR (blue line) and IIR-M (green line) are very similar. A comparison of the orbit types between them does not show any significant differences. Table 5 presents the ADEV of the selected intervals for two analysed orbit types. For τ = 1 and 4 s, the frequency stabilities are close to 1 × 10 −12 s; for τ = 16 and 64 s, every stability is close to 1 × 10 −13 s; and for bigger averaging intervals stabilities-close to even 7-8 × 10 −14 s.

Allan Deviation: All τ Method
An all τ analysis, demanding superior computational power than standard intervals, is usually applied as a consecutive power of 2, or multiplies of selected numbers. However, this method shows a wider spectrum of oscillator deviations in other than successive powers of the number 2, or other rigidly determined numbers and their multiples.
The many τ (averaging times) method shows similar results to the τ = 2 n method, however a greater noise appears, especially for the higher averaging time values. The clock corrections based on precise orbits are slightly more accurate than those based on broadcast ones. Frequency stabilities generally agree with each other for both the 2 n and many τ method, except for PRNs 26 and 27 (Block IIF). Analysis of the 2 n and τ method graphs ( Figure 3) show a very high consistency for all the averaging times between 1 and 256 s included, except for PRN18, which is less stable between 32 s to 128 s averaging times. The reason of the PRN18 stability degradation is not clear at present. For the two highest averaging times, the results are very random; in addition, they are not correlated with the satellite block type. A comparison of the analysed blocks shows almost identical stabilities for Blocks IIR and IIR-M (except for slightly worse results for PRN18). The Block IIF oscillators are slightly more stable than the other two blocks. In the case of the satellite orbit used for the calculations, there are no significant differences, excepting PRN 26 and 27 for the longest 1024 s τ on the basis of BRDC (GPS broadcast ephemeris file) calculations (Figures 3 and 4, top left). Figure 4 presents ADEV calculated using consecutive (with 1 s step) averaging times beginning from 1. This approach is aimed at detecting a possible clock anomaly that was not detected by the 2 n averaging times. If detected, it might have been responsible for smoothing the results. The most visible is the PRN26 oscillator stability graph using the broadcast orbit ( Figure 4, top left), which is jagged compared to the other graphs, but in Figure 3 the same phenomena is also visible for this satellite only. Figure 5 presents TDEV of clock products based on the broadcast and ultra-rapid orbits together with the error bars for a G26 satellites. TDEV varies between 3.9×10 -12 s to 2.9×10 -10 s for broadcast orbit and between 4.0×10 -12 s to 1.9×10 -10 s in case of ultra-rapid one. Analysis shows a similar behavior for a IIR and IIR-M GPS satellites, while IIF satellites (G26 and G27) are characterized by lower stability for the shortest intervals (1 s and 2 s), clearly higher for longer intervals (4-128 s), while for the highest intervals the stability of TDEV is the same. The same can be observed in the case of uncertainty estimates, error bars of results acquired by using broadcast orbit are much bigger than using ultra rapid ones.

Time Deviation: All τ Method
TDEV values of broadcast and ultra-rapid orbit using the many τ method are visible in Figure  6. The use of the many τ method shows the same results and conclusions as described in Section 4.3. Compared to Figure 4, there are no significant fluctuations in stability between successive intervals; this is due to the fact that TDEV is based on modification of the ADEV. Otherwise, TDEV presents time stability of phase versus observation interval.
G31 G26 G27 Figure 6. ADEV depending on orbit type using the many τ method.

Conclusions
The novelty of this work compared with other studies and the current state of knowledge results from three elements: an analysis of the orbit type impact on satellite clock stability, the usage of the Allan's variance all τ method for this type of elaboration, and the high sampling frequency on the basis of 1 s observation data on the IGS stations. The main purpose of the study presented was the determination of the effect of the orbit type used on the stability of GPS clocks. Calculations in BSW on the basis of three precise orbits (ultra-rapid, rapid, and final) show very similar results, thus precise orbits available in less time (ultra-rapid in real time, rapid after 17-41 h) than the final one (after 10-18 days) can be used to analyse clock stability. Unlike in precise positioning, where final orbit is used, ultra-rapid or rapid orbit might be used without loss of results quality for this kind of elaboration. Moreover, in this paper so-called all τ method was use, unlike in all the other similar publications. In this paper, apart from the most widely used averaging times as a consecutive powers of 2 (2 0 s = 1 s, 2 1 s = 2 s, 2 2 = 4 s, etc.) or 10 (10 0 s = 1 s, 10 1 s = 10 s, 10 2 s = 100 s, etc.) or 30/300 s (based on existing clock products), the authors adopted 1 s data. The third element supporting the novelty of this paper compared to existing research is the fact of using clock corrections based on modelling using BSW software with 1 Hz observation frequency. In a future work, the authors plan to analyse clock corrections with higher frequencies than 1 Hz, up to 20 Hz.
An analysis of high-rate 1 Hz satellite clock correction is presented in this paper. IGS's analysis centres provide rapid and final clock products with 5 s and 30 s correction intervals. For precise elaborations such as PPP, seismology or high-rate GNSS monitoring, 1 s interval clock corrections are required. In this paper, the authors show 1 s elaboration of satellite and stations clock corrections depending on the adopted orbit type and the method of ADEV and TDEV averaging time used. The emphasis has been put on the GPS satellites clock's stability depending on the type of orbit used. The use of GPS-only signals results from the fact that the adopted software does not allow for the use of BeiDou, Galileo, and QZSS signals in this type of analysis. Currently, Bernese GNSS Software gives the opportunity to compute only GPS and GLONASS (Russian: ГЛОНАСС -Глобальная навигационная спутниковая система, transliteration: Globalnaya navigatsionnaya sputnikovaya sistema) satellite clock corrections. Moreover, this work should be treated as preliminary research in the field of determining the impact of the orbit type in the satellite clock stability analysis. The next part of the research will include analysis using GNSS signals (GPS GLONASS) and, in the future, using all currently available satellite signals (GPS, GLONASS, Galileo, BeiDou, and QZSS). Four different orbit types were analysed: three precise with 15 min sampling intervals, and broadcast ephemeris with a 30 min sampling interval. The results from the data indicate GPS satellites frequency stability (ADEV) at a level of 1.7 × 10 −11 s to 7.2 × 10 −14 s and time stability (TDEV) at a level of 2.9 × 10 10 s to 3.9 × 10 12 s. The results also show the insignificant advantage of the satellite Block IIF's frequency and time stabilities in regard to the other two blocks. As used in this article, the many τ method allows for detection of the satellite clock anomaly. This method might also assist in preventing the use of a satellite clock with identified potential instability patterns as a reference clock.