The Design of Kalman Smoothers for Global Positioning System Data

Measurements of clocks aboard global positioning system (GPS) satellites as well as GPS system time are used for positioning, navigation, and timing. Different biases in a sequence of time difference measurements repeated each sidereal day produce apparent diurnal effects in the data. A composite time and frequency Kalman estimator is used here. A one day average of frequency at each measurement allows updates of clocks at time intervals of less than one day while aliasing diurnal variations. This technique is applied to data taken many times per day at time standards laboratories around the world according to a tracking schedule issued by the Bureau International des Poids et Mesures (BIPM). In this paper, we compute Kalman smoothed estimates of the time of clocks in the standards labs to define a global time base, then compute Kalman smoothed estimates of the time of GPS clocks against this time base.


I. INTRODUCTION HE MEASUREMENTS of international time stan-
T dards against the global positioning system (GPS) clocks provide many opportunities for interesting studies.First, we can study the relative behavior of the reference clocks themselves.Second, we can study the GPS as a whole.Since the reference clocks are extremely stable, their differences can be known to better than 10 ns with measurements only a few times per day.Since all of the international timing laboratories use the same tracking schedule, satellites are tracked simultaneously at several locations, and these tracks are repeated each day when the satellites are nominally in the same location relative to the rotating earth.It was shown previously that these data can be used to separate variances of various noise terms [l].
Here we use these variances to set noise parameters in Kalman smoothers to obtain linear estimates of the time of clocks.We make two different kinds of estimates: 1) estimates of the time of clocks at the standards labs against the time scale generated by the algorithm, and 2) estimates of the time of clocks aboard GPS satellites and of GPS system time against the scale.The first type of estimate is used to define a global time base to measure the time of the satellite clocks against as they circle the globe.The second type of estimate is smoothed from these measurements of the satellite clocks against the global time standard.
We use two types of Kalman smoothers: a Kalman time smoother ' and a mixed Kalman time-and-frequency smoother.The time smoother is used initially as a standard technique for obtaining optimal estimates of clock states in the presence of white phase noise.From this we discover significant biases in the data.The mixed timeand-frequency smoother allows us to alias the biases while still using measurements at time intervals of less than one day.
11.A GLOBAL TIME BASE Measurements of GPS clocks against time standards laboratories can be used for time transfer between the labs with cancellation of common view errors.In this way we can transfer time with stabilities of a few nanoseconds, accuracies of the order of 10 ns, and frequency stabilities of the order of one part in 1014 for measurement times of about four days and longer [2].This technique is done between pairs of locations with data averaged over 1 day (d).There are biases between measurements during the day [2], that make it difficult to compare clocks at intervals less than once per day.
We design a Kalman smoother using common view time measurements, hence, a Kalman time smoother, filtering forward and smoothing backward, in a manner similar to others [2]- [4].The smoother defines a time scale with estimates of three states against the scale for each clock: time offset, frequency offset, and frequency drift.We employ a model of white frequency modulation (FM), and random walk FM for each clock, and white phase noise (PM) as measurement noise.The level of measurement noise is determined for each track using the multistation separation of variance technique [ 11.The problem here is that biases in measurements repeated once per day produce an apparent diurnal variation [5].
We use nine stations in the smoother: the reference standards at the GPS receivers at each of the National Bureau of Standards (NBS), Boulder, CO; the Jet Propulsion Laboratory (JPL) Deep Space Network station at Goldstone, CA; the Applied Physics Laboratory (APL), Laurel, MD; the U.S. Naval Observatory (USNO), Washington DC; the Paris Observatory (OP), Paris, France; the Physikalisch Technische Bundesanstalt (PTB) in Braunschweig, Federal Republic of Germany; the Tokyo Astronomical Observatory (TAO) and the Communications Research Lab (CRL) both in Tokyo, Japan; and the WWVH radio station of NBS in Kauai, Hawaii.The data covers the 50-d period from February 19 to April 9, 1988, modified Julian days (MJD's) 47210 to 47260.Measurements U.S. Government work not protected by U.S. copyright are made regularly at each of these locations of GPS satellite clocks against the local clock averaged for 13 min, using the space vehicle's (SV) transmitted ephemeris and ionospheric model.In addition a transmitted SV clock correction is applied to this measurement to obtain GPS system time against the local clock.These measurements are repeated every sidereal day.Since the SV's are in 12 h sidereal orbits, we thus maintain essentially the same geometry of each SV measurement.Biases between measurements using different SV's or different geometries appear as diurnal variations, and are aliased away with 1-d averaging times.
Figs. 1-3 show a few of the results of the time smoother, estimating JPL's hydrogen maser against UTC (NBS).We show, respectively, a plot of the time of JPL's clock minus UTC (NBS), the square root of the Allan variance ( ~~( 7) ) of this data, and the magnitude of the Fourier spectrum of the JPL clock minus OTC (NBS).We see how the diurnal variations in the data corrupt our attempts to characterize the clocks.Our estimate of time offsets in Fig. 1 have diurnal ripples, and the variance plots, Figs. 2 and 3 have a sin (x)/x deviation with the high point of the deviation at 1 / 2 d, the low points at multiples of 1 d [ 5 ] .Thus we are able to see oY( 7 ) at 7 values of 1 d and longer, but not at less than 1 d.The magnitude of the Fourier transform plot of the time data, Fig. 3, shows the "bright line" at 1 d, a Fourier component with an amplitude of about 1 ns.

A COMPOSITE TIME AND FREQUENCY SMOOTHER
The measurements of the reference clock minus the SV clock, using a fixed geometry and SV, which are repeated once per sidereal day, have a high stability.Any consistent bias in the measurement is cancelled if we subtract yesterday's measurement from today's; this gives us an average frequency measurement over 1 d.If we average time measurements over 1 d we may update our smoothed estimates at 1 d and longer with minimal effects from biases.We may use our first difference of time, today's minus yesterday's, to update our estimates of frequency for periods less than 1 d.We use these frequency estimates to interpolate our time estimates between 1-d measurements.Fig. 4 illustrates the use of l-d frequency averages to update time over intervals less than 1 d.If we simply filter forward with this method, we would still have a diurnal Fourier component, due to our sudden change of time estimate with our daily time update.Since this is not a real time filter, we may use the benefit of hindsight and smooth backwards as well, allowing our time updates to influence the future and the past in a smooth way.
There are two complications with the composite time and frequency ( T / F ) smoother design.First, we are using some of the same information twice.This must be compensated for by breaking the covariance matrix into two parts: the frequency part, and the time and cross timefrequency part.We update the former part only with frequency measurements, and update the latter part with once-per-day time updates.If this is done, from simula- tion we have seen that the covariance matrix gives a realistic estimate of the smoother's knowledge of the states.
A second complication is in estimating the measurement noise for the smoother.Since we are updating with frequency averaged over 1 d, we are prefiltering the data, removing some of the white PM measurement noise, as well as some of the high frequency white FM noise of the clocks themselves.Of course, without this averaging, the biases prevent our seeing the clocks altogether.However, some information about the clocks is lost.The measurement noise needs to be lowered to account for the prefilter.The prefilter of 1-d frequency averaging allows noise to pass through proportional to the length of time since the previous update, up to 1 d.The clocks themselves are filtered also by using this 1-d frequency average.To the extent that the dominant noise at integration times of 1 /4 to 1 d is white noise FM, the square root of the Allan variance, ay ( T ) , will need to be increased by a slope of

T -~/ ~
extending back from 1 d.This design of smoother was applied to both simulated and real data.Fig. 5 shows typical U,( 7 ) results on simulated data, for the difference between a pair of clocks.The figure has the "true" value based on the generated data in hexagons, and the ( T / F ) smoothed estimate in X's.We see that the smoother can see the behavior of the clocks down to 1 d.For higher frequencies, some of the clock behavior has been filtered by the 1-d frequency average.The variance estimates are too small.The plus symbols indicate the values obtained by increasing ay ( 7 ) by a slope of 7-'12 extending back from 1 d, multiplying the 1 / 2 value by A, and the 1/4-d value by 2. We see that we have indeed recovered the clock differences, in a variance sense, down to 1 / 4 day.Fig. 6 shows the stability of the scale itself, the hexagons, and some of the simulated ensemble of clocks.We see that the scale is better than any individual clock.

IV. A STUDY OF THE GPS
The SV's are denoted by two different numbers: the Navstar number for the sequential order in which the SV was launched, called the SVN number and the pseudorandom code number (PRN), corresponding to the code the SV transmits.For this work we use Navstar (PRN) numbers 3 (6), 6 (9), 8 ( l l ) , 9 (13), 10 (12), and 11 (3), and the GPS system clock from MJD 47210 to 47260, February 19-April 9, 1988.
First we look at the results on simulated data in Fig. 10.Again we find that the prefilter of 1-d frequency av- Fig. 8.The square root of the Allan variance of the data in Fig. 7.We have an estimate of the relative behavior of these clocks for integration times less than 1 d.The estimates are too low for integration times less than one day, due to the one day frequency averages, and weie adjusted up by the reciprocal of the square root of the integration time in days.These values are represented by the plus symbols.
eraging has decreased the white FM level.The theoretical level of the variance has been reconstructed at 1 / 2 d and 1/4 d as before, indicated in the plus symbols.The data are really not dense enough to see the clock at 1 /4 d.We use Navstar 10, PRN 12, as an example of results based on actual clocks.Fig. 11 shows the plot of the estimated time offset of Navstar 10 from the Kalman time smoother against the global time scale.We see a significant diurnal systematic error.Fig. 12 shows the magnitude of the Fourier transform of this data, giving an estimate of this systematic at about 10 ns.Figs. 13 and 14 show the parallel plots using the T / F smoother.The systematic is gone.corrected values at 1 / 2 and 1 /4 d, though there is not dense enough data to see the clock at 1 /4 d.Fig. 16 shows the behavior of the GPS system time from the T / F smoother.This gives a lower limit to the possible stability of time available from GPS.As before, the plus symbols indicate the reconstructed ay ( 7 ) values for integration times less than 1 d.In this case the data is dense enough to see the GPS time at 1 / 4 d and perhaps even 1/8 d.Fig. 17 gives the ~~( 7 ) for a T / F forward filter of GPS time without a backward smoother.This shows the potential for a real-time clock based on GPS.The steering algorithm applied to the local clock based on such filtered measurements could improve short term stability, with the results in Fig. 16 giving a lower limit at about the level.-13 -14 4 5 6

LCG TAU (Seconds)
Fig. 17.The fractional frequency stability of GPS time after applying our time and frequency Kalman filter in the forward direction only.This shows the behavior of the GPS plus T / F Kalman filter as a real time clock.An appropriate steering algorithm designed around the local oscillator could help the short term behavior, with the limit being the performance illustrated in Fig. 16.

V. CONCLUSIONS
The GPS can be used to link reference clocks around the world producing a global time scale at a few parts in at integration times of 1 d.Biases in the satellites produce a diurnal variation with amplitude of the order of 1 ns.This can be circumvented using 1-d frequency averages to update time for intervals less than 1 d.The results gives us estimates of the fractional frequency stability of clock differences down to integration times of 1 /2 d.
The global time scale we generate by linking reference clocks can be used to look at the behavior of satellite clocks.Here we see diurnal effects in clocks of the order of 10 ns.Again we can use a composite time and frequency Kalman smoother to alias the diurnal systematic error.We obtain estimates of the behavior of clocks at integration times less than 1 d, but these are biased low.It is probably possible to estimate the error in our variance estimates by increasing O,,(T) by a slope of T -' / ~ extending back from 1 d.
Finally we look at the behavior of GPS system time.The time and frequency smoothed estimates of this "clock" has its ay ( 7 ) values at about for integration times from l /2 d to 4 d.This is the stability available to users of GPS with post-processing of time data.Using a forward filter only we get stabilities of the order of lo-'* at integration times of 1/2 d.The quality of a real-time clock based on GPS could be better than this, depending strongly on the local oscillator and the steering algorithm.

Figs. 7 SIKMTEDFig. 5 .Fig. 6 .
Fig.5.The time and frequency Kalman smoothed estimate of reference clock differences for simulated data.The X ' s are the estimates, the hexagons the generated data.The estimates are too low for integration times less than one day, due to the one day frequency averages, and were adjusted up by the reciprocal of the square root of the integration time in days.These values are represented by the plus symbols.

Fig. 7 .
Fig. 7.The time offset of JPL against UTC(NBS) as measured by GPS and estimated by a composite time and frequency Kalman smoother.Diurnal ripples apparent in Fig. 1 are now gone.

Fig. 9 .
Fig.9.The magnitude of the Fourier transform of the data in Fig.7.There is no apparent diurnal component.

Fig. 11 .Fig. 12 .Fig. 14 .
Fig. 11.The time estimates of Navstar 10 against the global time scale from the time Kalman smoother with a mean first difference and mean time removed.Diurnal fluctuations are clearly visible.

Fig 16
Fig 16The fractional frequency stability of GPS time after applying our time and frequency Kalman smoother This shows the limits of performance available to a user Again, at integration times less than 1 day the values have been reconstructed as represented by the plus symbols.