Autonomous clock ensemble algorithm for GNSS applications

2.1 Noise spectral density parameters

The white noise processes of the discrete clock model of Equation (1) characterize the stochastic properties of the clock for this model. It is these processes that are predominantly responsible for error in timing systems. We quantify the magnitude of these processes with spectral density values

3

For each component, the subscript indicates the clock state to which the noise process corresponds. Since environmental periodics require two states, we assume the same level of random walk in both noise processes per periodic. The units for the elements of S are S_θ (sec²); S_p,S₁,S₂ (sec²/sec); S_f (sec²/sec³); and S_r sec²/sec⁵.

Equation (4) shows that there is a relationship between each spectral density values of Equation (3) and the Hadamard variance . Note that in the timing literature, the Allan and Hadamard measures are generally applied to average frequency measurements and not phase.

4

Derivations of these relations can be found in Hutsell.^[11]

3.1 Clock state propagation

An expansion of the state transition matrix Φ must also be taken into account in order that it is size-compatible with the expanded state vector, noise vector, and also spectral density matrix S. Since each clock in the collection has the same general model as the others, the structure of the larger transition matrix has the form Φ_E(τ,t_k)=Φ(τ,t_k)⊗I_N where ⊗ is the Kronecker Product. Then, the equations used to propagate the main state vector and its associated state error covariance matrix P are

5

At the propagation stage, the clock states are all advanced forward along the parabolic model as though the drift of the clock is constant. The covariance terms are also propagated and then inflated by the level of process noise in each individual clock state. The process noise inflation is contained in Q(τ,t_k), which consists of accumulated process noise over the interval T between the two epochs. This is defined by the integral

6

The derivation of this integral is involved and can be reviewed further in Brown and Hwang.^[12]

3.2 Measurement update

After the clock state estimates have been propagated, they are compared to actual measurements observed at the epoch t_k. One clock among the member clocks must serve as the reference clock for the data set and that reference can change at each epoch. The ensemble filter is impervious to the reference given. We denote that clock with index r∈C and then establish the vector of measurements:

Here, the notation indicates the measurement of clock i with respect to reference clock r at the epoch t_k. The clock measurements that are ingested by this filter are phase measurements of one clock with respect to another. Hence, the observation model is

where y(t_k) is the vector of measurement errors having covariance matrix R(t_k) = E [y(t_k)y(t_k)^T ]. The observation matrix H takes the form

There is a rank deficiency in the observation model, and therefore, additional constraints on the timescale’s estimates must be introduced. We make the additional requirement that the weighted sum of all member clock corrections must be zero. This ensures that the ensemble does not drift in frequency over time. The equation of these constraints is

Four of these constraints are applied to the observation model, one for each of the independent primary clock states.^[13] The environmental periodic states, which are largely deterministic, are not constrained. The coefficient weighting terms, w_θ, w_p, w_f, and w_r, are independent weight distributions that are assigned to favor clocks in each process noise type. This key multi-weighting technique allows the clock ensemble to exploit the most stable clock(s) at all averaging intervals.

Although dynamic and responsive to process noise parameters, a clock’s weight in any of the ensemble’s state averages is limited to the range [0, 2.5/N] where N is the number of active clocks in the ensemble. An initial distribution uses only the process noise parameters and weights clocks according to the inverse of that state process parameter. A recursive algorithm shifts the distribution several times until each weight is at or below an upper threshold. The use of 2.5/N as an upper limit is an arbitrary choice– larger than 1/N, but well under 1/2 when N gets larger than five clocks. This upper limit of weight prevents any one clock from dominating the ensemble average. Since the Kalman filter itself is a recursive process and its state estimate vector has the ensemble time as reference, an unbounded weight condition will typically yield a solution where one clock becomes more and more stable with respect to the ensemble and eventually attains 100% of the ensemble weight. Also, to protect the ensemble from poorly performing clocks, the weighting conditions allow poor clocks to fall to 0 weight, and thus, there is no lower limit on weights.

The clock weights are dynamic and can change at each epoch. Impulse weighting is a method that allows the new weights to take effect immediately at the epoch of computation. Another method, termed the Gradient Weighting Method, moves the weight distribution towards the newly calculated weights. As described in Coleman,^[14] this weighting approach limits the changes in clock weight from one epoch to the next so that dramatic shifts in the ensemble’s composition do not take place. Protecting the clock state estimates and covariances from sudden impulses, which can happen when the ensemble composition changes suddenly, helps to guard against short term and localized ensemble instabilities.

After determination of clock weights, the Kalman Gain is calculated. This gain matrix dissociates all the noise accumulated in the entire system of clocks and apportions appropriate updates to all clock states based on the covariances thus far estimated. The update equations associated with this gain matrix are

7

Figure 2 shows a diagram of one entire step of the recursive process for this clock ensemble algorithm. Additional details of the Kalman filter approach and stochastic methods applied here can be found in Maybeck.^[15]

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

Recursive process of Kalman filter estimation of clock states

4.1 Clock break and anomaly detection

Throughout operation, a clock is bound to generate measurements that can be classified as outliers or breaks. In the case of real-time operation, a suspect measurement cannot immediately dictate the type of behavior the clock is undergoing. We therefore devise an algorithm to handle clock measurements and assess whether a clock is generating data that are consistent with the clock model of Equation (1), or data that are indicative of a broken or modified clock signal.

To assess whether a clock may have had a phase break, we take note of measurements that deviate substantially from the clock model prediction. These differences are exactly the pre-fit residuals, and they are key indicators of whether a clock may have had a phase break.

From a statistical perspective, pre-fit residuals are normal random variables with mean zero and variance that depends on the noise parameters of the clock and the measurement noise. Pre-fit residuals that are consistently several times larger than the variance are indicative of measurements being affected by non-random influences. The actual cause is of little concern to the filter process only that some event occurred which has shifted the expected clock measurements. The pre-fit residual is

8

where is the variance of the residual and has the form

9

While the magnitude of the pre-fit residual is important, the error covariance matrix P evolves and thus so does the variance here. In order that we can assess all clocks’ residuals using a common scale, we use the normalized pre-fit residual

10

The value of n_i(t_k) is also the number of σ_i(t_k) values by which the measurement deviates from the clock model. We use this value to determine whether the measurement received is an outlier. Typically, this function is set to declare outliers at |n_i(t_k)|>5. Unlike all other variables in the ensemble filter, the sequence of normalized residuals is saved for the past 1000 epochs, which allows a break detection algorithm to search the history of any clock’s behavior.

The following list of clock status types is available for assignment to each clock at each epoch. The assignment depends on the clock’s status at the previous epoch, the recent history of statuses and the present value of n_i(t_k).

• Steady state. Clock measurement is consistent with its model and clock state estimates are stable (determined by state error variances changing over time by less than a small user-defined tolerance);

• Outlier. Deviant clock measurement although a recent history of steady-state measurements is observed;

• Phase break. Suspected after a continuous string of N_p outlier measurements;

• Frequency break. Suspected after a continuous string of N_f phase break alerts;

• Clock reset. Failure of states and covariance to converge on clock measurements after intervention;

• Missing data. Status when no data available for given clock.

In each case of break or reset, the filter responds with modifications to the clock’s state or covariance matrix in an attempt to re-converge the states to the actual observations of the clock. When these mitigation efforts fail, the clock is reset by putting initial state values of zero and large values of variance on each state. The Kalman filter then has a chance to re-converge on new state values.

In the case that a phase break is suspect, the propagation model of Equation (5) is modified to include an impulse term δx, which is the observed discrepancy between the clock state and measurement. Hence,

11

Since phase breaks are easily observed, a correction of this type is typically all that is required for the clock to resume accurate estimation with low pre-fit residuals. If the residuals continue to grow, then the breaking routine proceeds to attempt a frequency break correction.

Frequency breaks are not easily observed by the phase data since the phase series remains continuous across the break. To correct for a frequency break, the filter is allowed to search for and converge on a new frequency value by adding additional process noise at the propagation stage in the clock’s frequency state. This is accomplished by adding δS to the process noise matrix of Equation (6) so that

This one-time addition of process noise inflates the state error covariance matrix so that the filter will accept a wider range of frequency terms and hopefully steer towards a new and correct clock frequency value.

If neither of these two techniques are able to correct a clock’s growing pre-fit residuals, then all eight of this clock’s states are set to 0 and the clock’s covariance values are reset to a large diagonal matrix. This initialization allows the Kalman filter to start a new estimation recursion on this clock in search of a new solution while all other clocks continue to be filtered members of the ensemble. In the case of a failure in the clock hardware, no such autonomous technique will be valuable, and intervention by an operator will be required to handle and correct the solution.

A break in the drift state would be realized in the phase data as an inflection change. Since the filter is taking only phase measurements, a sudden change of this type would be challenging to detect and even more so to correct for. At present, the automatic break detection does not attempt detection or repair to the drift state. Error in the drift state could be more problematic for the ensemble stability than for the clock estimate.

4.2 Clock parameter adaptation

The clock spectral density parameters of Equation (3) can be adaptively updated at each filter step using the clock’s pre-fit residuals. Several duplicate filters, each stepping at nominally different intervals, are implemented for this purpose. Since each clock will typically be dominated by a single noise type at a given averaging interval, a parallel implementation of the ensemble filter can run with a stepping interval that nominally yields a pre-fit residual sequence for the clock that is sensitive to longer term noise types.

The relationship between the variance of the pre-fit residuals and the clock parameters can be written out by utilizing Equations (5) and (9). We obtain

12

where the notation refers to the (i+(m−1)N,i+(n−1) N) entry of the matrix P. Developments and details regarding this approach can be found in Stein.^[16]

If the clock parameters are not correct for a given clock, then adjustments to the densities , and are sought such that Equation (12) holds. Back solving Equation (12) for the delta densities and accounting for the fact that each clock pre-fit residual will be biased low because of its contribution to the ensemble, it can be shown that

13

where α ∈ {θ,p, f, r} and are considered as a point estimate for . The exponent γ depends on the dominant noise type α that is under consideration and is equal to {0,1,3,5} for {θ,p,f,r}, respectively. Also, we define the ensemble average of the spectral densities, seen in Equation (13), by the following weighted sum

For a particular filter stepping interval, τ, only one of the noise-type densities is dominant. For this reason, one may wish to implement several filters running at nominally different time stepping interval so that each density parameter is assessed by some filter. The appropriate form of Equation (13) is chosen in each case. The updates should be made slowly so that the actual density modification is of the form

14

for noise type α and constant β ≪ 1. This constraint is important since severe changes to the clock density parameters cause immediate and large changes to the clock ensemble weights, which can lead to short-term instability.^[14] Severe automatic break alerts can also ensue if the density parameters should become too small. Furthermore, this adaptation should only be employed once the clock has reached steady state in the filter. Note that total phase white noise density may not be distinguished from the white (phase) measurement noise R_i, unless the measurement noise is substantially lower than the white phase density. Also, spectral density adaptation is not performed on those parameters for the environmental periodic states.

5.1 Test I: simulated GPS constellation

Accurate clock estimation and prediction is critical for a system such as GPS since user range errors of approximately 29.9 cm are experienced for each nanosecond of clock error. We test the clock ensemble algorithm with a simulation of the GPS constellation clocks as well as its monitor station clocks. The process noise parameters used to simulate the clock behavior are based on data from the IGS clock products representative of the GPS network at the start of 2018. This simulation is generated for a 50-day period with simulated measurement noise having variance of 1×10⁻²⁰ sec² (equivalent to about 3-cm RMS) that is appropriate for a full GNSS constellation clock solution.^[17] In both the IGST and GPS time applications, the input of this ensembling software has such a clock difference solution and thus that level of noise would be appropriate here.

The clock ensemble must accurately estimate the phase, frequency, and drift of each clock in order for predictions to retain any degree of precision for future epochs. The relative error of the phase estimation for clock i at epoch t_k is calculated by

15

Some reference clock, REF, must be chosen to eliminate a constellation-wide bias. It is best to choose such a clock with the lowest noise parameters; estimates of that clock will give much less contribution to the relative errors owed to clock i. The frequency and drift state relative errors are similarly calculated by simply replacing the estimated and true values of p with estimated and true values of either f or r.

Figure 3 shows distributions of the relative errors computed via Equation (15) for phase, frequency, and drift. For all the distributions shown, the first 5-day period of the solution was discounted in order to dismiss the large errors often encountered during a cold start period. Note the phase converges better than the frequency and that better than the drift. The drift state often takes months to converge in a Kalman filter, especially since only measurements of the clock phase are available for the Kalman update stage.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Distribution of relative error of the estimated clock states as calculated by Equation (15). One series shown for each of phase, frequency, and drift. Note that 1e−2 marks the position of 0.01 (1% error) and 1e−1 the position of 0.1 (10% error)

5.2 Test II: GPS life test clock reference

In earlier literature, NRL presented performance results of the Block IIF life test clocks; see Vannicola et al.^[2,18] In these cases, the reference clock chosen was from among the many hydrogen masers within the NRL PCEF. The older generation of masers was a collection that included MHM2000, SAO, and SigmaTau masers. Stability performance of the reference would therefore be limited by the stability of the best local hydrogen maser. In the past, the best maser was chosen as a reference, but risk is involved with this scenario as the ongoing testing platform becomes reliant on single clock members for a stable and break-free reference.

The reference timescale at NRL’s PCEF that was launched in 2014 was able to maintain sufficiently low drift such that all test clocks could be measured and analyzed in both long- and short-term averaging intervals without regard to the performance of any individual maser.^[4] A more recent acquisition of Microsemi’s MHM2010 masers (identified by H Masers 5, 6, 7, 8) in the laboratory has expanded the collection of clocks that contribute to the ensemble. These new masers have large drift offsets, which translates to the reference ensemble’s drift as well; reducing the timescale drift requires accurate estimates of the masers’ drift.^[19] Further, we must ensure that the ensemble’s drift remains close to a clock member whose longer term stability is confidently known to be low.

After the ensemble’s filter was given several months to converge on the drift estimates for the clock, we selected a period of data to analyze the performance of the ensemble. The data used below are the phase measurements from NRL’s PCEF dual mixer phase measurement system during the period 1 August through 31 December 2018. Table 1 shows the resultant drift estimates, process noise parameters applied, and dynamically determined weights for all contributing member clocks near the start of this data period. The drift of both the GPS IIF cesium clock and an NRL-USNO time link measurement are low and thus can be used to estimate severe drift in the ensemble average. The cesium clock cannot, however, contribute weight to the ensemble less it bias the reference towards its performance and hence also the reporting statistics. For this reason, there are N=8 contributing clocks in the ensemble.

While the weight of drift is skewed towards the USNO master clock link and the best performing long-term masers, short-term ensemble performance is better due to weights that favor the laboratory’s newest masers. Many of the new Microsemi masers exhibit excellent short-term performance with Hadamard deviation well below 1×10⁻¹⁵ at a 1-min averaging interval. Note that this test uses 2.5/N for the upper weight limit. Since N=8 is the number of contributing ensemble members, we see that USNO’s w_r weight achieves the maximum possible value. This is expected since UTC(USNO) is steered to maintain long-term stability.

Data collected during August and September contained a number of phase breaking events that were detected by the ensemble and corrected for immediately. These are catalogued in Table 2. The causes of these breaks are typically the result of hardware issues or necessary lab maintenance that interrupts a clock signal. Generally, the phase correction is negative the pre-fit residual. By the correction model of Equation (11), the corrected clock estimates should reduce the clock’s residuals. If no frequency or drift breaks components exist and if no other non-quadratic behavior takes place, the clock should return as an ensemble member within several epochs of data. Throughout these breaking events, the ensemble remains continuous.

Figure 4 shows the clock residuals over the 3-month period of 1 August through 31 October. Note that the residuals go off the chart in each case as is expected by the value of the breaks determined and listed from Table 2. Event 2 involved uncharacteristic changes that did not fit the clock model of the earlier section nor the breaking model of Equation (11). With the phase correction largely compensating the break, however, the filter was able to re-estimate the clock states accurately within 1 day.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Pre-fit residuals, as computed by Equation (8), for three of the ensemble’s contributing member clocks. Each of these clocks had a breaking event during these 3 months of the test; they are marked in the figure and exhibit an off-chart value for the residual at the time of the clock-phase jump. Event numbers correspond to those listed in Table 2

The breaking event between events 2 and 3 was a complete restart of the measurement system. On rare occasions, a new signal added to the measurement system can cause the system to hang up, requiring a restart. A restart of the system leaves arbitrary discontinuities in the entire set of measurements. Since the entire set of clock difference observations changes, this requires a manual re-alignment of the data outside of the ensemble’s break detection process. Future work could involve using observations from a back-up system to solve for one pair and recover the continuity of the entire clock difference observation set.

Stability of the timescale and associated corrected clock estimates is plotted in Figure 5. The stability curves shown are of time series that are raw measurements (with breaks corrected a posteriori) re-aligned to the timescale reference rather than the measurement system’s physical reference clock. The period of data analyzed in this figure is from 1 August through 31 December 2018. Note that the stability of the laboratory masers versus the reference timescale are well below the stability assessed for the test clocks. The ensemble’s clock weights are used to create an estimate of the ensemble’s stability—also shown in Figure 5. The stability is largely unaffected by the two more sizable breaks listed in Table 2.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

Stability of NRL laboratory clocks. Reference for all clock measurement series is the NRL timescale. The curve marked with black squares shows the estimated ensemble stability based upon the mean weight contribution of each clock

5.3 TEST III: IGS clock product reference

As a significant part of the IGS program for the geophysical community that initially collected GPS tracking data from hundreds of sites around the world, a coherent reference time was needed for these data and subsequent geophysical analyses. The IGS clock products were to provide the time reference for all these worldwide tracking data. The products were initially formed by combining multiple clock solutions from different analysis centers around the globe.^[20] This clock combination resulted in a measure of improved measurement precision for clocks for both stations and GPS satellites. However, the resulting demonstrated performance for this combination, made from a daily linear alignment to broadcast GPST, had unfavorable stability from one day to the next. In 2003, an early version of a time reference, known as IGS Time (IGST), for the clock products was launched to improve IGST for these products and eliminate the severe impulses to the daily solutions.^[21] The result was a greatly improved IGST stability for these clock products and consequently for the analyses dependent upon them.

Since the initial release, a later version of the ensemble algorithm has been built into the IGS processing. That version contains the features that have been presented earlier including the break routine algorithm. The environmental periodic states are an important feature of this algorithm for GNSS applications. Figure 6 shows a typical estimation for a GPS satellite with respect to the IGST ensemble. The upper plot shows both the total phase and phase estimates demonstrating the ability of the ensemble algorithm to isolate the environmental harmonic states. These periodic states have frequencies of once and twice per revolution and are quite common to the GPS constellation clocks. The lower plot is the pre-fit residual sequence for this clock. With the exception of one outlier data point, this sequence shows that the estimation of this clock is consistent with the clock measurements observed.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Example of satellite clock product for G25, which is the IGS’s rapid solution for the operational clock on-board GPS satellite on PRN 25. In the upper plot, the offset from IGRT timescale is shown for both the total phase and phase estimates; a quadratic polynomial was removed from both series. In the lower plot, normalized prefit residual sequence is shown as calculated in Equation (10)

Detection and response to breaks and missing data is a critical feature for the IGS clock ensemble process. There are presently 388 clocks listed in the ensemble, but only 100–150 will have data to filter on any given day. The architecture of the IGS processing and combination software puts the choice of ground network used in the solution to the IGS’ Analysis Centers. Only clocks chosen by sufficiently many centers are admitted to the IGS’s solution.^[17] In addition, there is a boundary jump in the clock solution each day since the orbit and clock combinations are performed in daily batches. For these reasons, clock data in the IGS often contains anomalies, especially for ground stations. The existing break detection algorithm can repair many breaks.

Figure 7 shows an example of the automated break response and correction routine. The IGS rapid timescale (IGRT) detected a phase break in the clock solution at IGS station HOB2 (Hobart, Australia) in mid February 2018. All estimates and measurements plotted in Figure 7 are with respect to the calculated IGRT (IGS Rapid Timescale) reference time. Linear polynomials are removed from both the estimate and measurement series; one from each side of the break so that the details are clearer in the plot. The difference between the P₁ and P₂ coefficients are 5.44×10⁻⁴ sec and 1.89×10⁻¹⁴ sec/sec for the constant and linear, respectively. The actual break value is therefore about 500 microseconds. Such a jump could be caused by a quick swap of the clock signal at the site’s receiver, or a restart of the clock, for example.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

Example of breaking event on a IGS network station clock. Top: phase estimates and measurements around the time of the break; Middle: clock status values assigned as the event unfolds; Bottom: normalized pre-fit residual sequence that drives the breaking procedure

The first two plots of Figure 7 demonstrate that the filter can continue to produce clock estimates in the absence of data, as seen around the end of day 13. At the start of day 14, data resume at a new phase offset yielding large normalized residuals that far exceed the window [−5,5] of the bottom plot. The filter therefore recognizes these data as outliers until after some tolerance when it implements a phase impulse as in Equation (11). This automated response to the station’s clock phase jump allows the timescale to keep that clock’s contribution to the ensemble upon recovery of the correct state values. Without an automated recovery of the state values, this clock would need to drop from the ensemble average thereby reducing the ensemble’s stability in the medium term around the event. Despite the demonstration here, the IGS processing often occurs daily without any human input, so future improvements and increased robustness would be desirable for the IGS’s clock ensemble computation.