Abstract
The GPS composite clock defines GPS time, the timescale used today in GPS operations. GPS time is illuminated by examination of its role in the complete estimation and control problem relative to UTC/TAI. The phase of each GPS clock is unobservable from GPS pseudorange measurements, and the mean phase of the GPS clock ensemble (GPS time) is unobservable. A new and useful observability definition is presented, together with new observability theorems, to demonstrate explicitly that GPS time is unobservable. Simulated GPS clock phase and frequency deviations, and simulated GPS pseudorange measurements, are used to understand GPS time in terms of Kalman filter estimation errors.
1. Introduction
GPS time is created by processing GPS pseudorange measurements with the operational GPS Kalman filter. Brown [2] refers to the object created by the Kalman filter as the GPS composite clock, and to GPS time as the implicit ensemble mean phase of the GPS composite clock. The fundamental goal by the USAF and the USNO is to control GPS time to within a specified bound of UTC/TAI. (I refer to TAI/UTC understanding that UTC has an accumulated discontinuity (a sum of leap seconds) when compared to TAI. But unique two-way transformations between TAI and UTC have been in successful operational use since 1972. I have no need herein to further distinguish between TAI and UTC.) I present here a quantitative analysis of the GPS composite clock, derived from detailed simulations and associated graphics. GPS clock diffusion coefficient values used here were derived from Allan deviation graphs presented by Oaks et al. [12] in 1998. I refer to them as “realistic,” and in the sequel I claim “realistic” results from their use. Figure 1 presents their diffusion coefficient values and my derivation of associated Allan deviation lines.
My interest in the GPS composite clock derives from my interest in performing real-time orbit determination for GPS NAVSTAR spacecraft from ground receiver pseudorange measurements. (James R Wright is the architect of ODTK (Orbit Determination Tool Kit), a commercial software product offered by Analytical Graphics, Inc. (AGI).) The estimation of NAVSTAR orbits would be incomplete without the simultaneous estimation of GPS clock parameters. I use simulated GPS clock phase and frequency deviations, and simulated GPS pseudorange measurements, to study Kalman filter estimation errors.
This paper was first prepared for TimeNav'07 [20]. I am indebted to Charles Greenhall (JPL) for encouragement and help in this work.
2. The Complete Estimation and Control Problem
The USNO operates two UTC/TAI master clocks, each of which provides access to an estimate of UTC/TAI in real time (1 pps). One of these clocks is maintained at the USNO, and the other is maintained at Schriever Air Force Base in Colorado Springs. This enables the USNO to compare UTC/TAI to the phase of each GPS orbital NAVSTAR clock via GPS pseudorange measurements, by using a UTC/TAI master clock in a USNO GPS ground receiver. Each GPS clock is a member of (internal to) the GPS ensemble of clocks, but the USNO master clock is external to the GPS ensemble of clocks. Because of this, the difference between UTC/TAI and the phase of each NAVSTAR GPS clock is observable. This difference can be (and is) estimated and quantified. The root mean square (RMS) on these differences quantifies the difference between UTC/TAI and GPS time. Inspection of the differences between UTC/TAI and the phase of each NAVSTAR GPS clock enables the USNO to identify GPS clocks that require particular frequency-rate control corrections. Use of this knowledge enables the USAF to adjust frequency rates of selected GPS clocks. Currently, the USAF uses an automated bang-bang controller on frequency-rate. (According to Bill Feess, an improvement in control can be achieved by replacing the existing “bang-bang controller” with a “proportional controller.”)
3. Stochastic Clock Physics
The most significant stochastic clock physics are understood in terms of Wiener processes and their integrals. Clock physics are characterized by particular values of clock-dependent diffusion coefficients, and are conveniently studied with aid of a relevant clock model that relates diffusion coefficient values to their underlying Wiener processes. For my presentation here I have selected “The clock model and its relationship with the Allan and related variances" presented as an IEEE paper by Zucca and Tavella [19] in 2005. Except for FM flicker noise, this model captures the most significant physics for all GPS clocks. I simulate and validate GPS pseudorange measurements using simulated phase deviations and simulated frequency deviations, according to Zucca and Tavella.
4. Kalman Filters
I present my approach for the optimal sequential estimation of clock deviation states and their error covariance functions. Sequential state estimates are generated recursively from two multidimensional stochastic update functions, the time update (TU) and the measurement update (MU). The TU moves the state estimate and covariance forward with time, accumulating integrals of random clock deviation process noise in the covariance. The MU is performed at a fixed measurement time where the state estimate and covariance are corrected with new observation information.
The sequential estimation of GPS clock deviations requires the development of a linear TU and nonlinear MU. The nonlinear MU must be linearized locally to enable application of the linear Kalman MU. Kalman's MU [9] derives from Sherman's theorem [11, 15, 16], Sherman's theorem derives from Anderson's theorem [1], and Anderson's theorem derives from the Brunn-Minkowki inequality theorem [5, 17]. The theoretical foundation for my linearized MU derives from these theorems.
4.1. Initial Conditions
Initialization of all sequential estimators requires the use of an initial state estimate column matrix and an intial state estimate error covariance matrix for time .
4.2. Kalman Filter: Linear TU and Linear MU
Derivation and calculation for the discrete-time Kalman filter, linear in both TU and MU, is best presented by Meditch [11, Chapter 5].
4.3. Linear TU and Nonlinear MU
The simultaneous sequential estimation of GPS clock phase and frequency deviation parameters can be studied with the development of a linear TU and nonlinear MU for the clock state estimate subset. This is useful to study clock parameter estimation, as demonstrated in Section 6.
Let denote an column matrix of state estimate components, where the left subscript denotes state epoch and the right subscript denotes time-tag for the last observation processed, where . Let denote an associated square symmetric state estimate error covariance matrix (positive eigenvalues).
4.3.1. Linear TU
For , the propagation of the true unknown matrix state is given bywhere is called the process noise matrix. Propagation of the known matrix state estimate is given bybecause the conditional mean of is zero. Propagation of the known matrix state estimate error covariance matrix is given bywhere the matrix is called the process noise covariance matrix (see [19] for concrete clock examples of and ).
4.3.2. Nonlinear MU
Calculate the matrix filter gain : The filter measurement update state estimate matrix , due to the observation , is calculated withwhere is the scalar variance on the observation residual , and is a nonlinear function of . Define the error in :Define the state estimate error covariance matrix with Bucy and Joseph [3, page 141] recommend that should be calculated withwhereEquations (8) and (9) reduce to the form given by Kalman:Calculation of by (8) and (9) is numerically stable, whereas the Kalman form calculation is not.
4.4. Nonlinear TU and Nonlinear MU
Refer to Section 4.3.2 for the nonlinear MU.
4.4.1. Nonlinear TU
The nonlinear TU always spans a nonempty time interval and requires the use of a numerical state estimate integrator . Given an initial time , a final time , and a force model , then propagates the state estimate from to using forces to get . That is,This can be shortened to writewhere the use of forces is tacitly implied. Thus, is a column matrix with elements:
4.5. Kalman Filter Advantage
Severe computational problems are incurred in any attempt to estimate unobservable states using iterated batch least squares methods or iterated maximum likelihood methods for navigation, because state-sized inversions of singular matrices are required. Here the Kalman filter is distinguished in that estimates of unobservable states can be created and used without matrix inversion problems because the Kalman filter MU is free of state-sized matrix inversions.
By design, one typically estimates observable states. But the Kalman filter enables one to create unobservable states. The USAF chose to create unobservable GPS clock parameter states for construction of GPS time.
5. Observability
I have defined observability in terms of a Kalman filter formulation, and I have proved simple theorems related thereto. My definition of observability is different than Kalman's definition and, unlike Kalman's definition, is directly applicable to covariance matrices derived from a Kalman filter.
5.1. Definition
If the state estimate error variance of a particular state estimate component is reduced by processing an observation, then that state estimate component is observable to that observation. Otherwise, that state estimate component is not observable (unobservable) to that observation.
Theorem 1. If every component of the row matrix of measurement-state partial derivatives is zero at time , then every component of the state estimate is unobservable at time .
Proof. implies that according to (10). Thus none of the variances of are reduced due to processing the observation . Then by definition, is unobservable in every component.
Theorem 2. Given values for scalars , , at time , and given that , then the scalar state estimate is observable at time .
Proof. The obvious inequality implies thatMultiply through by : Add 1: Multiply through by : Now use (4) and (10) to write Insert (18) into the inequality (17) to get the result: Thus the variance is reduced due to processing the observation . Then the scalar state is observable by definition.
5.2. Theoretical Foundation
These theorems are referred to expressions given by Kalman for filter gain and covariance , see (4) and (10). Kalman's expressions are derived from the rigorous theorem chain provided by Sherman, Anderson, and Brunn-Minkowski—the theoretical foundation is deep.
5.3. Determine Observability Directly
Given an optimal sequential estimator, given a particular collection of applicable observations (real or simulated), and given realistic state estimate error covariance matrices and at each time , apply the definition of observability directly (note that this is impossible using Kalman's definition of observability) to distinguish between observable and unobservable state elements. An optimal sequential estimator is designed to eliminate significant aliasing between estimated state elements, and thus enables this distinction.
6. Unobservable GPS Clock States
GPS time is created by the operational USAF Kalman filter by processing GPS pseudorange observations. GPS time is the mean phase of an ensemble of many GPS clocks, and yet the clock phase of every operational GPS clock is unobservable from GPS pseudorange observations, as demonstrated below. GPS NAVSTAR orbit parameters are observable from GPS pseudorange observations. The USAF Kalman filter simultaneously estimates orbit parameters and clock parameters from GPS pseudorange observations, so the state estimate is partitioned in this manner into a subset of unobservable clock parameters and a subset of observable orbit parameters. This partition is performed by application of Sherman's theorem in the MU.
6.1. GPS Pseudorange Representation
Let denote time of radio wave transmission for the th NAVSTAR clock, and let denote time of radio wave receipt for the th ground station clock. (Refer all times to a coordinate time, e.g., to GPS time. Appropriate transformations between proper time and coordinate time must be performed in the operational algorithms, but state estimate observability is independent of relativity, so observability can be defined and discussed independent of relativity.) Let and denote Kalman filter estimation errors in clock phase for and . Define time of transmission difference and time of receipt difference :Thus,Equation (21) present and as additive combinations of deterministic times and and Kalman filter estimation errors in clock phase and . Define the one-way GPS pseudorange measurement : Insert (21) into (22):where is speed of light in vacuum. Define Then,where is deterministic and is random.
6.2. Partition of Kalman Filter Estimation Errors
Let denote the phase component of Kalman filter estimation error that is common to every GPS ensemble clock, when it exists. Define phase differences and with for ground station and NAVSTAR . Then Kalman filter estimation errors , , for ground station clocks and , , for NAVSTAR clocks have the additive partition:(This partition was introduced by Brown [2].)
6.3. The Common Random Phase Component Is Unobservable
Insert (28) into (25): Insert (29) into (26):Thus, the random phase component that is common to the Kalman filter estimation error for every ensemble clock has vanished in the range representation . Variations in cannot cause variations in :because the partial derivative is zero: An application of Theorem 1 to (32) demonstrates that is unobservable from .
But the architect who designs the complete estimator must design an optimal NAVSTAR orbit estimator to prevent aliasing from NAVSTAR orbit estimation errors into . It helps to know that there is no coupling between the orbit and in the complete state transition function. I have provided a new method herein to identify this aliasing, and I have provided suggestions on where to look for inadequate modeling that would be the source of this aliasing. See Section 9.
6.4. Independent Random Phase Components Are Observable
The independent phase deviations and are observable to GPS pseudorange observations because their partial derivatives are nonzero:Estimation of and by the Kalman filter will reduce their error variances.
6.5. Partition of KF1 Estimation Errors
Subtract estimated clock deviations from simulated (true) clock deviations to define and quantify Kalman filter (KF1) estimation errors. Adopt Brown's additive partition of KF1 estimation errors into two components. I refer to the first component as the unobservable error common to each clock (UECC), and to the second component as the observable error independent for each clock (OEIC). (Observability is meaningful here only when processing simulated GPS pseudorange data.) See (28). is the UECC, is the OEIC for ground station clocks and is the OEIC for NAVSTAR clocks. On initialization of KF1, the variances on the UECC and OEIC are identical. On processing the first GPS pseudorange measurements with KF1 the variances on both fall quickly. But with continued measurement processing the variances on the UECC increase without bound while the variances on the OEIC appoach zero asymptotically.
For simulated GPS pseudorange data I create an optimal sequential estimate of the UECC by application of a second Kalman filter KF2 to pseudomeasurements defined by the phase components of KF1 estimation errors.
Since there is no physical process noise on the UECC, an estimate of the UECC can also be achieved using a batch least squares estimation algorithm on the phase components of KF1 estimation errors–-demonstrated previously by Greenhall [7]. (I apply sufficient process noise covariance for KF2 to mask the effects of double-precision computer word truncation. Without this, KF2 does diverge.)
6.6. Unobservable Error Common to Each Clock
There are at least four techniques to estimate the UECC when simulating GPS pseudorange data. First, one could take the sample mean of KF1 estimation errors across the clock ensemble at each time and form a sample variance about the mean; this would yield a sequential sampling procedure, but where each mean and variance is sequentially unconnected. Second, one can employ Ken Brown's implicit ensemble mean (IEM) and covariance; this is a batch procedure requiring an inversion of the KF1 covariance matrix followed by a second matrix inversion of the modified covariance matrix inverse; this is not a sequential procedure. Third, one can adopt the new procedure by Greenhall [7] wherein KF1 phase estimation errors are treated as pseudomeasurements, and are processed by a batch least squares estimator to obtain optimal batch estimates and covariance matrices for the UECC. Fourth, one can treat the KF1 phase estimation errors as pseudomeasurements, invoke a second Kalman filter (KF2), and process these phase pseudomeasurements with KF2 to obtain optimal sequential estimates and variances for the UECC. I have been successful with this approach. Figure 3 presents an ensemble of “realistic" KF1 phase estimation errors, overlaid with “realistic" KF2 sequential estimates of UECC in phase. (By “realistic" I refer to realistic clock diffusion coefficient values.)
6.7. Observable Error Independent for Each Clock
At each applicable time subtract the estimate of the UECC from the KF1 phase deviation estimate, for each particular GPS clock, to estimate the OEIC in phase for that clock. During measurement processing, the OEIC is contained within an envelope of a few parts of a nanosecond (see Figure 4).
Figure 4 presents a graph of two cases of the OEIC for ground station clock S1. For the blue line of Figure 4, intervals of link visibility and KF1 range measurement processing are clearly distinguished from propagation intervals with no measurements. During measurement processing, the observable component of KF1 estimation error is contained within an envelope of a few parts of a nanosecond.
Calculation of the sequential covariance for the OEIC requires a matrix value for the cross-covariance between the KF1 phase deviation estimation error and the UECC estimation error at each time. I have not yet been able to calculate this cross-covariance.
7. Allan Variance and PPN Relations
7.1. Allan Coefficients versus Diffusion Coefficients
Denote as clock averaging time, as Allan variance, as Allan's FMWN coefficient, as Allan's FMRW coefficient, as the FMWN diffusion coefficient, and as the FMRW diffusion coefficient. Then,where
7.2. Proportionate Process Noise (PPN)
Let denote a variable to identify each GPS clock in an ensemble of clocks. For each clock define the ratio between diffusion coefficients and :Then PPN is defined when, for each GPS clock and each associated ratio , we have
7.3. Case 12
The calculation of , , according to the diffusion coefficient values presented in Figure 1 shows that PPN is not satisfied for Case 12:
8. Kalman Filters KF1 and KF2
I have simulated GPS pseudorange measurements for two GPS ground station clocks S1 and S2, and for two GPS NAVSTAR clocks N1 and N2. Here I set simulated measurement time granularity to 30 s for the set of all visible link intervals. Visible and non-visible intervals are clearly evident in the blue line of Figure 4. I set the scalar root-variance for both measurement simulations and Kalman filter KF1 to cm. Typically m for GPS pseudorange, but when carrier phase measurements are processed simultaneously with pseudorange, the root-variance is reduced by two orders of magnitude. So the use of cm enables me to quantify lower performance bounds for the simultaneous processing of both measurement types.
8.1. Create GPS Clock Ensemble
Typically, one processes measurements with a Kalman filter to derive sequential estimates of a multidimensional observable state. Instead, here I imitate the GPS operational procedure and process simulated GPS pseudorange measurements with KF1 to create a sequence of unobservable multidimensional clock state estimates. Clock state components are unobservable from GPS pseudorange measurements. See Figure 2 for an example of an ensemble of estimated unobservable clock phase deviation state components created by KF1.
8.1.1. Sherman's Theorem
GPS time, the unobservable GPS clock ensemble mean phase, is created by the use of Sherman's theorem [11, 18] in the USAF Kalman filter measurement update algorithm on GPS range measurements. Satisfaction of Sherman's Theorem guarantees that the mean-squared state estimate error on each observable state estimate component is minimized. But the mean-squared state estimate error on each unobservable state estimate component is not reduced. Thus the unobservable clock phase deviation state estimate component common to every GPS clock is isolated by application of Sherman's theorem. An ensemble of unobservable state estimate components is thus created by Sherman's theorem—see Figure 3 for an example.
8.2. Initial Condition Errors
A significant result emerges due to the modeling of Kalman filter (KF1) initial condition errors in phase and frequency. Initial estimated clock phase deviations are significantly displaced by the KF1 initial condition errors in phase. As time evolves estimated clock phase deviation magnitudes diverge continuously and increasingly when referred to true (simulated) phase deviations, and this is due to filter initial condition errors in frequency. See Figure 2 for an example.
9. Identify Nonclock Modeling Errors
My interest in the GPS NAVSTAR (SV) orbit determination problem, combined with that of the clock parameter estimation problem, has enabled the identification of a useful diagnostic tool: given realistic values for diffusion coefficients for each of the real GPS clocks, then quantitative upper bounds can be calculated on OEIC magnitudes. These calculations require the use of a rigorous simulator. Existence of significant cross-correlations between GPS clock phase errors and other nonclock GPS estimation modeling errors enables significant aliasing into GPS clock phase estimates during operation of KF1 on real data. But given rigorous quantitative upper bounds on OEIC magnitudes, then significant violation of these bounds when processing real GPS pseudorange and carrier phase data identifies nonclock modeling errors related to the GPS estimation model. Modeling error candidates here include NAVSTAR orbit force modeling errors, ground antenna modeling errors (multipath), and tropospheric modeling errors. NAVSTAR orbit force modeling errors include those of solar photon pressure, albedo, thermal dump, and propellant outgassing. The accuracy of this diagnostic tool depends on the use of realistic clock diffusion coefficient values and a rigorous clock model simulation capability.
10. Observable Clocks
In an earlier version of my paper, I reported on KF1 validation results where clock S1 was specified as a TAI/UTC clock, external to the GPS clock ensemble consisting of S2, N1, and N2. This brought observablity (see Sections 5 and 6 herein) to S2, N1, and N2 clock states from GPS pseudorange measurements, drove clocks S2, N1, and N2 immediately to the TAI/UTC timescale, and enabled a clean validation of my filter implementation. Also it raised the question: why not the same thing for the real GPS clock ensemble? Discussions with Ed Powers (USNO) and Bill Feess (Aerospace Corporation) reveal that this approach was tried and discarded after the difficulty in recovery from an uplink hardware failure was blamed on the use of a single TAI/UTC Master Clock. This issue was resolved with Kenneth Brown’s introduction of the implicit ensemble mean. The mean phase (GPS time) of the GPS clock ensemble will remain unobservable to GPS pseudorange measurements in the USAF Kalman filter for the foreseeable future.