## Abstract

The US Naval Research Laboratory (NRL) has developed clock ensembles for a number of vital GNSS projects and services throughout the US and international community. Some efforts include algorithms for local realization of Coordinated Universal Time; participation in the upgrade of the Global Positioning System (GPS) system time; the initial implementation of a coherent geophysical reference time known as IGS Time (IGST); and development of various Department of Defense (DoD) capabilities for the dissemination and transfer of precise time. This paper will cover the major algorithm components that compute and maintain a robust and stable time reference for the most stringent applications mentioned here. These components include automated alert and response to clock anomalies and reduction of measurement weight in response to large clock measurement residuals. Both real and simulated data have been used to validate the ensemble algorithm and its products for reliable and robust usage in these environments.

## 1 INTRODUCTION

Precise and accurate local reference time is an important capability required of any system that expects to achieve and maintain precise synchronization of the system elements or components containing clocks. GNSS is a prime example where the individual system’s reference time enables effective synchronization of its member users onboard clocks enabling the operation of the system. These same techniques to achieve precise and accurate local time are used not only to achieve precise synchronization and stability but also provide the means of maintaining continuity in the reference time by assessing or measuring the individual clock performance in a particular system.

The ensemble algorithm that is used in this analysis stems from that developed and refined at the US Naval Research Laboratory (NRL) over the past decade. NRL has several contributions to the technologies for generating system reference time, contributions to technologies used in international timescale development, and timing applications for precise dissemination of time for ground and space in general. Perhaps the most notable initial effort in this area was known as the Common Time Reference (CTR) project.^{[1]} This project investigated the capabilities needed to generate and maintain a local time reference among a Naval Task Group for maintaining synchronization of their communications and navigation systems. A necessary part of this development was the ability to combine the performance of the varying types of clocks deployed in the Task Group. This need stimulated the development of algorithms to combine the ensembling techniques used in classical timescale generation with system time reference technology for locally realized precise reference time. Coupled with tactical data links to support the intercomparison of the clock outputs, the ability to determine a system reference time across a Task Group that could be compared and steered to the global reference timescale of the US Naval Observatory, UTC(USNO) was realized. This development continued in the development of algorithms for the realization of IGS Time and local time reference systems that could be steered to UTC(USNO) for synchronization in operations and for clock performance evaluation.

Since 1997, NRL has been life testing space-qualified clocks from the NAVSTAR satellites for the GPS Directorate. Initially, two rubidium clocks from the GPS Block IIR satellites were subjected to life test in a specially designed laboratory at NRL for almost 8 years. Subsequent life testing was performed on two cesium beams and two rubidium clocks from the GPS Block IIF satellites beginning in 2008. Currently, one Block IIF rubidium and one Block IIF cesium clock are under life test.^{[2]} GNSS clocks have improved dramatically over the past few decades with Hadamard Deviation at a 1-day averaging interval improving from 8×10^{−14} for Block II and IIA cesium clocks to 2×10^{−14} in Block IIR PerkinElmer Rubidium Atomic Frequency Standards (RAFS).^{[3]} Continued improvements have been seen in GPS’s Block IIF RAFS clocks, which exhibit Hadamard deviations at levels of 5×10^{−15}.^{[2]} For this reason, the ability to measure them both in ground testing and on-orbit has continually required improved time references.^{[4]}

More recent applications of this work have entailed generating a reference time for clock products of Global Navigation Satellite Systems (GNSS) in several efforts. The first was on the platform of the International GNSS Service (IGS) where a reference clock ensemble is generated to serve as a continuous reference with steering and offsets known to UTC(*k*) laboratories. In another instance, the ensemble algorithm was outfitted to serve as the next generation GPS system time as part of the OCX effort.^{[5]} In both instances, outlier data, breaks in the clock solution, or frequency anomalies on a given clock can spoil the continuity and reliability of the generated clock estimates. This has motivated the most recent features of this algorithm, which seek to mitigate the effects of clock anomalies in the most automated way possible; they are presented in the later part of this paper.

In Section 2, we present the clock model that is used to represent a single clock, its states, and stochastic components. This clock model can be used in the simulation of clocks, or a system of clocks in order to evaluate the performance of the overall algorithmic approach. Characteristics of the underlying noise processes in this clock are also discussed as well as their relationships to common clock statistics. This model is sufficiently flexible that it can represent both environmentally controlled laboratory clocks as well as on-orbit GPS clocks presently in service.

In Section 3, the clock model is expanded to a collection of clocks, and the algorithm used to construct a clock ensemble (timescale) is shown in detail. The ensemble algorithm is a Kalman filter, and the description there includes the main cycle of clock state prediction and correction with Kalman Gain, the use of measurements in the ensemble filter to help generate the timescale reference, and the weighing constraints which determine the formation of the ensemble.

Following the main algorithm, in Section 4, we present auxiliary functions utilized for clock anomaly detection and parameter adaptation. Several details that tie these functions to the main Kalman filter loop are shown to clarify how the filter’s output is beneficial to clock error detection. These two functions, as described above, are recent capabilities that have been fielded to several of the scenarios described, but not previously presented. Particulars about the operation of the break detection algorithm and adaptive parameter control are highlighted better here.

Finally, in Section 5, case studies are presented that demonstrate this timescale algorithm’s ability to provide a robust time reference especially in the presence of anomalous data. In the first case, the algorithm is used for a large block of GPS simulated clock data. This case demonstrates convergence to the correct clock difference solutions which can be assessed since the true states were simulated, hence known. In the second case, the algorithm is used to generate a local time reference for precise clocks under special laboratory testing inside NRL’s Precise Clock Evaluation Facility (PCEF).^{[6]} In the third and final case, we test the algorithm with the International GNSS Service (IGS) clock products data, which contains substantially noisier data and larger anomalies than in laboratory data. Although the present IGS rapid and final clock combinations only provide GPS solutions, it is expected that future IGS service will include clocks from other GNSS including Beidou, GLONASS, and Galileo.^{[7]}

## 2 CLOCK MODEL

Both the clock ensemble and clock simulation algorithms utilize a clock model based on standard three-state models developed in the 1980s in several different works including Tryon and Jones^{[8]} and Jones and Tryon.^{[9]} In such models, the phase, frequency, and drift of a clock are included with random walk noise processes associated with each. We expand the model used here to a four-state model to allow for an additional white noise process over the phase state, which yields what we term a *total* phase state. This state also encompasses up to two environmental periodics that affect the clock and/or its signal. The total phase state does not necessarily model or separate lower level electronic noise such as that from a phase lock loop or time keeping system (TKS), but rather allows some stochastic behavior in this state should it be present in raw measurements.

Throughout this report, we define the state vector, **x**, and corresponding noise input vector, **u**, which contain the variables of the states and input noises for a clock. We define these terms below along with the corresponding units.

Each of the noise inputs are white noise sequences, which we will assume to have mean zero. The diagram of Figure 1 shows the dynamics for these eight clock states and how the corresponding white noise sequences integrate to random walks and random runs. Note the relationship of the clock states in Figure 1 that highlights the clock phase (seconds), clock frequency (dimensionless), and clock frequency drift (1/seconds) which are termed phase, frequency, and drift, respectively, throughout this analysis.

Given a discretized set of fixed epochs {*t*_{1},*t*_{2},…,*t*_{k},…}, we can define an interval of continuous time *T*=[ *t*_{k−1},*t*_{k}], which also has the associated epoch-to-epoch time step *τ* = |*T*|. Advancing the clock states from epoch *t*_{k−1} to *t*_{k} is accomplished by the discretized propagation equation

1

where the state transition and process noise pre-multiplier matrices are given by

2

respectively. Here, **0**_{n} is the size *n*×*n* zero matrix, and **I**_{m} is the size *m*×*m* identity matrix. The block components of this transition matrix consist of the basic clock model and the clock environmental periodics given by the matrices

where *ω*_{i} are the frequencies of the environmental periodics. Environmental periodics have been observed in the GPS constellation, and the frequencies of those particular periodics are *ω*_{1}=2.003 and *ω*_{2}=4.006 cycles per day.^{[10]} Since this model will be built into a Kalman filter, flicker noise is not included in the clock model here.

### 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^{2}); *S*_{p},*S*_{1},*S*_{2} (sec^{2}/sec); *S*_{f} (sec^{2}/sec^{3}); and *S*_{r} sec^{2}/sec^{5}.

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 CLOCK ENSEMBLE FORMATION

The NRL clock ensemble algorithm uses the model of the previous section in a Kalman filter in order to form a stable clock ensemble (timescale) reference. In such a formulation, the states of the clocks are predicted forward from one epoch to the next using the clock model and updated using the Kalman Gain. To separate these two steps at each epoch, there is the notion of an *a priori* epoch (before the measurement update) and an *a posteriori* epoch (after the measurement update), which are denoted by , respectively.

For the collection of clocks *C* = {1,2,3,…,*N*}, we build a state vector that contains the states of all *N* clocks rather than just a single clock. For the remainder of this report, we therefore consider the state vector

Here, dependence on the epoch value *t*_{k} is assumed. The clock state estimate vector should not be confused with the true clock states **x**(*t*_{k}). The true states have a reference time of “perfect time,” which can only ever be available to our computations in a simulation. In application, true states are not known just as a “perfect time” reference is never fundamentally tangible. The reference time for the estimates, meanwhile, is exactly the ensemble product that our algorithm realizes.

A white noise vector **u** that has the same size as the expanded vector above is also part of this ensemble model. This noise vector has the same role as in the single clock model. The clock state error covariance matrix for the entire ensemble is similarly large and is defined by the expectation

The covariances of this matrix are estimated in the Kalman filter; these also cannot be known exactly in application since the true states, **x**, are unknown.

### 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]}

## 4 AUXILIARY FUNCTIONS

Some unique features of the clock ensemble algorithm at NRL consist of routines that adapt parameters and react to deviant clock measurements. Each component has the same goal in mind: protect the clock ensemble from bad clock measurements or events. To support this goal, these routines must determine when a clock’s measurements are poor or detrimental to the ensemble and then issue alerts that can cause a user to intervene, or an automated function to downweight measurements in the Kalman update stage.

### 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 ALGORITHM TEST AND EVALUATION

We present three sets of data each occurring on a different platform on which we can test the response of the NRL ensemble algorithm. These data sets differ from one another in complexity, measurement noise, and anomalous data. In the first set, we use a simulation of the GPS constellation and its associated monitor stations, which is free of uncharacteristic clock behavior. This data set can be best employed to assess the timescale’s convergence on correct clock states—even in the presence of external environmental harmonics. In the second test, we show the results of stability for timescale members in a laboratory testing environment where measurement noise in considerably low and the environment is controlled. The final test uses a set of IGS data to form an ensemble and identify clock anomalies. Breaks, outlier, and noisy data are more common in GNSS-type applications, so this test serves as a good test bed for the break control routine.

### 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^{−20} sec^{2} (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.

### 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^{−15} 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.

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.

### 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.

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*_{1} and *P*_{2} coefficients are 5.44×10^{−4} sec and 1.89×10^{−14} 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.

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.

## 6 CONCLUSION

A Kalman filter-based timescale that models clock phase, frequency, and drift is a useful tool for reference in a variety of different applications. In cases where a reference time is essential, the Kalman filter approach offers the advantage of providing real-time service. Accurate estimation of the member clocks and maintenance of the ensemble’s stability are the primary tasks of importance, and they can be disturbed by a host of different events that may introduce discontinuities in the clock observations.

The auxiliary functions of this clock ensemble algorithm improve the reliability of the reference time produced. In the tests presented, we found that the ensemble product is able to continue with near autonomous operation and deliver a stable timescale that is suitable for the application at hand. The laboratory setting at NRL has proven to be the most stable platform thus far where intervention is now needed approximately five to ten times annually. The stability holds at several parts in 10^{16} with contributing masers that have stability several parts higher. In the laboratory setting, future work will consider methods to handle a complete measurement system loss. Such approaches will require other observations, and they would need to be continuously available if the ultimate goal is greater autonomy of the ensemble generation process.

For geodetic and GNSS applications, autonomy is a greater challenge as anomalies can be more common. This ensemble algorithm has been installed in the present IGS architecture and has demonstrated good clock state estimation and some break detection. Nonetheless, future work requires improvements to the frequency break response and management of absent data. Future improvements to these algorithms and perhaps other supporting functions could help to improve the robustness of the IGS ensemble processing.

The test cases and algorithm design that were demonstrated in this paper show that this clock ensemble and resultant time reference provides suitable, break-free, and stable reference for several applications. Serving as an accurate local system reference time for a variety of platforms or timing and navigation technologies in a wider array of environments requires some additional configurations and designs. Continued development and reliance upon GNSS systems and products will only further the need for better and faster synchronization of system components for improved operability in the future.

## HOW TO CITE THIS ARTICLE

Coleman MJ, Beard RL. Autonomous clock ensemble algorithm for GNSS applications. *NAVIGATION-US*. 2020;67:333–346. https://doi.org/10.1002/navi.366

## Footnotes

This article is a U.S. Government work and is in the public domain in the USA.

- Received June 4, 2019.
- Revision received October 29, 2019.
- Accepted January 23, 2020.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.