Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

GPS Composite Clock Analysis

GPS Composite Clock Analysis Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2008, Article ID 261384, 8 pages doi:10.1155/2008/261384 Research Article James R. Wright Analytical Graphics, Inc., 220 Valley Creek Blvd, Exton, PA 19341, USA Correspondence should be addressed to James R. Wright, jwright@agi.com Received 30 June 2007; Accepted 6 November 2007 Recommended by Demetrios Matsakis 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. Copyright © 2008 James R. Wright. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1. INTRODUCTION The estimation of NAVSTAR orbits would be incomplete without the simultaneous estimation of GPS clock param- GPS time is created by processing GPS pseudorange measure- eters. I use simulated GPS clock phase and frequency devi- ments with the operational GPS Kalman filter. Brown [2] ations, and simulated GPS pseudorange measurements, to refers to the object created by the Kalman filter as the GPS study Kalman filter estimation errors. composite clock, and to GPS time as the implicit ensemble This paper was first prepared for TimeNav’07 [20]. I am mean phase of the GPS composite clock. The fundamental indebted to Charles Greenhall (JPL) for encouragement and goal by the USAF and the USNO is to control GPS time to help in this work. within a specified bound of UTC/TAI. (I refer to TAI/UTC understanding that UTC has an accumulated discontinuity 2. THE COMPLETE ESTIMATION AND (a sum of leap seconds) when compared to TAI. But unique CONTROL PROBLEM two-way transformations between TAI and UTC have been in successful operational use since 1972. I have no need herein The USNO operates two UTC/TAI master clocks, each of to further distinguish between TAI and UTC.) I present here whichprovidesaccess to an estimate of UTC/TAIinrealtime a quantitative analysis of the GPS composite clock, derived (1 pps). One of these clocks is maintained at the USNO, and from detailed simulations and associated graphics. GPS clock the other is maintained at Schriever Air Force Base in Col- diffusion coefficient values used here were derived from Al- orado Springs. This enables the USNO to compare UTC/TAI lan deviation graphs presented by Oaks et al. [12] in 1998. I to the phase of each GPS orbital NAVSTAR clock via GPS refer to them as “realistic,” and in the sequel I claim “realistic” pseudorange measurements, by using a UTC/TAI master results from their use. Figure 1 presents their diffusion coef- clockinaUSNO GPSgroundreceiver. Each GPSclock is ficient values and my derivation of associated Allan deviation a member of (internal to) the GPS ensemble of clocks, but lines. the USNO master clock is external to the GPS ensemble of My interest in the GPS composite clock derives from clocks. Because of this, the difference between UTC/TAI and my interest in performing real-time orbit determination for the phase of each NAVSTAR GPS clock is observable. This GPS NAVSTAR spacecraft from ground receiver pseudor- difference can be (and is) estimated and quantified. The root ange measurements. (James R Wright is the architect of mean square (RMS) on these differences quantifies the differ- ODTK (Orbit Determination Tool Kit), a commercial soft- ence between UTC/TAI and GPS time. Inspection of the dif- ware product offered by Analytical Graphics, Inc. (AGI).) ferences between UTC/TAI and the phase of each NAVSTAR 2 International Journal of Navigation and Observation Case 12 The sequential estimation of GPS clock deviations re- 1e−11 quires the development of a linear TU and nonlinear MU. 8 N 1 2-state clocks N 2 Diffusion coefficient values The nonlinear MU must be linearized locally to enable ap- √ √ ID σ1( s) σ2(1/ s) plication of the linear Kalman MU. Kalman’s MU [9]de- S13.23e−12 S2 S21.12e − 11 5.3e−17 rives from Sherman’s theorem [11, 15, 16], Sherman’s the- N14e − 11 1e−16 1e−12 N23e − 11 8e−17 orem derives from Anderson’s theorem [1], and Anderson’s S1 theorem derives from the Brunn-Minkowki inequality the- 2 orem [5, 17]. The theoretical foundation for my linearized Time durations T = 15552000 s MU derives from these theorems. 1e−13 τ = 32 s τ 0 = 8s 4.1. Initial conditions PRN seed number S11 1e−14 S2 30, 000, 001 Initialization of all sequential estimators requires the use of 90% ensemble last calc. 6 N 1 60, 000, 001 N 2 90, 000, 001 an initial state estimate column matrix X and an intial state 0|0 estimate error covariance matrix P for time t . 0|0 0 1 2 3 4 5 6 10 10 10 10 10 10 τ (s) 4.2. Kalman filter: linear TU and linear MU Figure 1: Allan deviation lines for S1, S2, N1, and N2. Derivation and calculation for the discrete-time Kalman fil- ter, linear in both TU and MU, is best presented by Meditch [11, Chapter 5]. GPS clock enables the USNO to identify GPS clocks that re- quire particular frequency-rate control corrections. Use of this knowledge enables the USAF to adjust frequency rates of 4.3. Linear TU and nonlinear MU selected GPS clocks. Currently, the USAF uses an automated bang-bang controller on frequency-rate. (According to Bill The simultaneous sequential estimation of GPS clock phase and frequency deviation parameters can be studied with the Feess, an improvement in control can be achieved by replac- development of a linear TU and nonlinear MU for the clock ing the existing “bang-bang controller” with a “proportional controller.”) state estimate subset. This is useful to study clock parameter estimation, as demonstrated in Section 6. Let X denote an n × 1 column matrix of state estimate j|i 3. STOCHASTIC CLOCK PHYSICS components, where the left subscript j denotes state epoch t and the right subscript i denotes time-tag t for the last ob- The most significant stochastic clock physics are under- i servation processed, where i, j ∈{0, 1, 2,...}.Let P denote stood in terms of Wiener processes and their integrals. j|i an associated n× n square symmetric state estimate error co- Clock physics are characterized by particular values of clock- variance matrix (positive eigenvalues). dependent diffusion coefficients, and are conveniently stud- ied with aid of a relevant clock model that relates diffusion coefficient values to their underlying Wiener processes. For 4.3.1. Linear TU my presentation here Ihaveselected“Theclock modeland its relationship with the Allan and related variances” pre- For k ∈{0, 1, 2, 3,... , M}, the propagation of the true un- sented as an IEEE paper by Zucca and Tavella [19] in 2005. known n × 1 matrix state X is given by Except for FM flicker noise, this model captures the most significant physics for all GPS clocks. I simulate and vali- X = Φ X + J ,(1) k+1 k+1,k k k+1,k date GPS pseudorange measurements using simulated phase deviations and simulated frequency deviations, according to where J is called the process noise matrix. Propagation of k+1,k Zucca and Tavella. the known n × 1 matrix state estimate X is given by k|k 4. KALMAN FILTERS X = Φ X (2) k+1|k k+1,k k|k I present my approach for the optimal sequential estimation because the conditional mean of J is zero. Propagation k+1,k of clock deviation states and their error covariance functions. of the known n × n matrix state estimate error covariance Sequential state estimates are generated recursively from two matrix P is given by k|k multidimensional stochastic update functions, the time up- date (TU) and the measurement update (MU). The TU moves P = Φ P Φ + Q ,(3) the state estimate and covariance forward with time, accu- k+1|k k+1,k k|k k+1,k k+1,k mulating integrals of random clock deviation process noise in the covariance. The MU is performed at a fixed measure- where the n × n matrix Q is called the process noise co- k+1,k ment time where the state estimate and covariance are cor- variance matrix (see [19] for concrete clock examples of J k+1,k rected with new observation information. and Q ). k+1,k James R. Wright 3 4.3.2. Nonlinear MU 4.5. Kalman filter advantage Calculate the n × 1matrixfilter gain K : Severe computational problems are incurred in any attempt k+1 to estimate unobservable states using iterated batch least −1 T T K = P H H P H + R . (4) k+1 k+1|k k+1 k+1|k k+1 k+1 k+1 squares methods or iterated maximum likelihood methods for navigation, because state-sized inversions of singular ma- The filter measurement update state estimate n × 1matrix trices are required. Here the Kalman filter is distinguished in X , due to the observation y , is calculated with k+1|k+1 k+1 that estimates of unobservable states can be created and used without matrix inversion problems because the Kalman filter X = X + K y − y X ,(5) k+1|k+1 k+1|k k+1 k+1 k+1|k MU is free of state-sized matrix inversions. where R is the scalar variance on the observation resid- k+1 By design, one typically estimates observable states. But ual y − y(X ), and y(X ) is a nonlinear function of k+1 k+1|k k+1|k the Kalman filter enables one to create unobservable states. X . Define the error ΔX in X : k+1|k k+1|k+1 k+1|k+1 The USAF chose to create unobservable GPS clock parameter states for construction of GPS time. ΔX = X − X . (6) k+1|k+1 k+1 k+1|k+1 Define the n × n state estimate error covariance matrix 5. OBSERVABILITY P with k+1|k+1 Ihavedefined observability in terms of a Kalman filter for- P = E ΔX ΔX . (7) k+1|k+1 k+1|k+1 k+1|k+1 mulation, and I have proved simple theorems related thereto. My definition of observability is different than Kalman’s def- Bucy and Joseph [3, page 141] recommend that P k+1|k+1 inition and, unlike Kalman’s definition, is directly applicable should be calculated with to covariance matrices derived from a Kalman filter. P = P − T,(8) k+1|k+1 k+1|k 5.1. Definition where T −1 If the state estimate error variance of a particular state es- T = P H R H P , k+1|k k+1 k+1|k k+1 k+1 (9) timate component is reduced by processing an observation, R = H P H + R . k+1 k+1 k+1|k k+1 k+1 then that state estimate component is observable to that ob- servation. Otherwise, that state estimate component is not Equations (8)and (9) reduce to the form given by Kalman: observable (unobservable) to that observation. P = I − K H P . (10) k+1|k+1 k+1 k+1 k+1|k Theorem 1. If every component of the row matrix H of k+1 Calculation of P by (8)and (9) is numerically stable, k+1|k+1 measurement-state partial derivatives is zero at time t , then k+1 whereas the Kalman form calculation is not. every component of the state estimate X is unobservable at k+1 time t . k+1 4.4. Nonlinear TU and nonlinear MU Proof. H = 0 implies that P = P according to k+1 k+1|k+1 k+1|k Refer to Section 4.3.2 for the nonlinear MU. (10). Thus none of the variances of P are reduced due to k+1|k processing the observation y . Then by definition, X is k+1 k+1 4.4.1. Nonlinear TU unobservable in every component. The nonlinear TU always spans a nonempty time interval Theorem 2. Given values for scalars H , P > 0, R > k+1 k+1|k k+1 and requires the use of a numerical state estimate integra- 0 at time t , and given that H = / 0, then the scalar state k+1 k+1 tor ϕ . Given an initial time t ,afinal time t ,and aforce 0 f estimate X is observable at time t . k+1 k+1 model u(X (t), t), then ϕ propagates the state estimate X (t ) from t to t using forces u(X (τ ), τ)toget X (t ). That is, 0 f f 2 Proof. The obvious inequality P H + R > k+1|k k+1 k+1 P H > 0 implies that k+1|k k+1 X t = ϕ t ; X t , t , u X (τ ), τ , t ≤ τ ≤ t . (11) f f 0 0 0 f This can be shortened to write P H k+1|k k+1 1 > > 0. (14) P H + R k+1|k k+1 k+1 X t = ϕ t ; X t , t , (12) f f 0 0 Multiply through by −1: where the use of forces u(X (τ ), τ ) is tacitly implied. Thus, ϕ is a column matrix with n elements: 2 P H k+1|k k+1 ⎡ ⎤ −1 < − < 0. (15) P H + R x1 k+1|k k+1 k+1 ⎢ ⎥ ⎢ x2⎥ ⎢ ⎥ Add 1: ⎢ ⎥ x3 ϕ = . (13) ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ P H . k+1|k k+1 0 < 1 − < 1. (16) ϕ P H + R xn k+1|k k+1 k+1 4 International Journal of Navigation and Observation Multiply through by P : in the operational algorithms, but state estimate observabil- k+1|k ity is independent of relativity, so observability can be de- Nh P H k+1|k k+1 fined and discussed independent of relativity.) Let δx and 0 < 1 − P <P . (17) k+1|k k+1|k Gi P H + R k+1|k k+1 δx denote Kalman filter estimation errors in clock phase k+1 Nh Gi Dh for t and t . Define time of transmission difference t T R T Now use (4)and (10)towrite Di andtimeofreceipt difference t : P H k+1|k k+1 Dh Nh Nh P = 1 − P . (18) k+1|k+1 k+1|k 2 t = t − δx , T T T P H + R k+1|k k+1 k+1 (20) Di Gi Gi t = t − δx . Insert (18) into the inequality (17) to get the result: R R R 0 <P <P . (19) Thus, k+1|k+1 k+1|k Nh Dh Nh Thus the variance P is reduced due to processing the ob- k+1|k t = t + δx , T T T servation y . Then the scalar state X is observable by def- (21) k+1 k+1 Gi Di Gi t = t + δx . inition. R R R Nh Gi Equation (21)present t and t as additive combinations of T R 5.2. Theoretical foundation Dh Di deterministic times t and t and Kalman filter estimation T R Nh Gi These theorems are referred to expressions given by Kalman errors in clock phase δx and δx . Define the one-way GPS T R pseudorange measurement ρ : for filter gain K and covariance P ,see (4)and (10). k+1 k+1|k+1 NhGi Kalman’s expressions are derived from the rigorous theo- Gi Nh Gi Nh rem chain provided by Sherman, Anderson, and Brunn- ρ = c t − t , t >t . (22) NhGi R T R T Minkowski—the theoretical foundation is deep. Insert (21) into (22): 5.3. Determine observability directly Di Gi Dh Nh ρ = c t + δx − t + δx NhGi R R T T (23) Given an optimal sequential estimator, given a particular col- Di Dh Gi Nh = c t − t + δx − δx , R T R T lection of applicable observations (real or simulated), and given realistic state estimate error covariance matrices P k+1|k where c is speed of light in vacuum. Define and P at each time t , apply the definition of observ- k+1|k+1 k+1 ability directly (note that this is impossible using Kalman’s Di Dh Δt = t − t , (24) R T definition of observability) to distinguish between observable hi Gi Nh and unobservable state elements. An optimal sequential es- δt = δx − δx . (25) R T timator is designed to eliminate significant aliasing between Then, estimated state elements, and thus enables this distinction. hi ρ = c Δt + δt , (26) NhGi 6. UNOBSERVABLE GPS CLOCK STATES where Δt is deterministic and δt is random. 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 6.2. Partition of Kalman filter estimation errors the clock phase of every operational GPS clock is unobserv- Let x denote the phase component of Kalman filter estima- able from GPS pseudorange observations, as demonstrated tion error that is common to every GPS ensemble clock, when below. GPS NAVSTAR orbit parameters are observable from Gi Nh it exists. Define phase differences x and x with OR OT GPS pseudorange observations. The USAF Kalman filter si- multaneously estimates orbit parameters and clock param- Gi Gi x = δx − x , OR R eters from GPS pseudorange observations, so the state esti- (27) mate is partitioned in this manner into a subset of unobserv- Nh Nh x = δx − x OT T able clock parameters and a subset of observable orbit pa- rameters. This partition is performed by application of Sher- for ground station i and NAVSTAR h. Then Kalman filter es- man’s theorem in the MU. Gi timation errors δx , i ∈{1, 2,...}, for ground station clocks Nh and δx , h ∈{1, 2,...}, for NAVSTAR clocks have the ad- 6.1. GPS pseudorange representation ditive partition: Nh Gi Gi Let t denote time of radio wave transmission for the hth δx = x + x , R C OR Gi NAVSTAR clock, and let t denote time of radio wave receipt (28) Nh Nh for the ith ground station clock. (Refer all times to a coordi- δx = x + x . T OT nate time, e.g., to GPS time. Appropriate transformations be- tween proper time and coordinate time must be performed (This partition was introduced by Brown [2].) James R. Wright 5 6.3. The common random phase component (UECC), and to the second component as the observable is unobservable error independent for each clock (OEIC). (Observability is meaningful here only when processing simulated GPS pseu- Gi Insert (28) into (25): dorange data.) See (28). x is the UECC, x is the OEIC OR Nh for ground station clocks and x is the OEIC for NAVS- OT hi Gi Nh δt = δx − δx R T TAR clocks. On initialization of KF1, the variances on the Gi Nh UECC and OEIC are identical. On processing the first GPS = x + x − x + x (29) C C OR OT pseudorange measurements with KF1 the variances on both Gi Nh = x − x . fall quickly. But with continued measurement processing the OR OT variances on the UECC increase without bound while the Insert (29) into (26): variances on the OEIC appoach zero asymptotically. For simulated GPS pseudorange data I create an optimal Gi Nh ρ = c Δt + x − x . (30) OR OT NhGi sequential estimate of the UECC by application of a second Kalman filter KF2 to pseudomeasurements defined by the Thus, the random phase component x that is common to phase components of KF1 estimation errors. the Kalman filter estimation error for every ensemble clock Since there is no physical process noise on the UECC, an has vanished in the range representation ρ .Variations NhGi estimate of the UECC can also be achieved using a batch least Δx in x cannot cause variations Δρ in ρ : C C NhGi NhGi squares estimation algorithm on the phase components of KF1 estimation errors—demonstrated previously by Green- ∂ρ NhGi Δρ = Δx (31) hall [7]. (I apply sufficient process noise covariance for KF2 NhGi ∂x to mask the effects of double-precision computer word trun- cation. Without this, KF2 does diverge.) because the partial derivative H = ∂ρ /∂x is zero: NhGi ∂ρ NhGi 6.6. Unobservable error common to each clock = 0. (32) ∂x There are at least four techniques to estimate the UECC An application of Theorem 1 to (32) demonstrates that x is C when simulating GPS pseudorange data. First, one could take unobservable from ρ . the sample mean of KF1 estimation errors across the clock NhGi But the architect who designs the complete estimator ensemble at each time and form a sample variance about must design an optimal NAVSTAR orbit estimator to prevent the mean; this would yield a sequential sampling procedure, aliasing from NAVSTAR orbit estimation errors into x .It C but where each mean and variance is sequentially uncon- helps to know that there is no coupling between the orbit and nected. Second, one can employ Ken Brown’s implicit ensem- x in the complete state transition function. I have provided C ble mean (IEM) and covariance; this is a batch procedure re- a new method herein to identify this aliasing, and I have pro- quiring an inversion of the KF1 covariance matrix followed vided suggestions on where to look for inadequate modeling by a second matrix inversion of the modified covariance ma- that would be the source of this aliasing. See Section 9. trix 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, 6.4. Independent random phase components and are processed by a batch least squares estimator to ob- are observable tain optimal batch estimates and covariance matrices for the Gi Nh The independent phase deviations x and x are observ- UECC. Fourth, one can treat the KF1 phase estimation er- OR OT able to GPS pseudorange observations because their partial rors as pseudomeasurements, invoke a second Kalman filter derivatives are nonzero: (KF2), and process these phase pseudomeasurements with KF2 to obtain optimal sequential estimates and variances ∂ρ NhGi = +c, for the UECC. I have been successful with this approach. Gi ∂x OR Figure 3 presents an ensemble of “realistic” KF1 phase esti- (33) ∂ρ mation errors, overlaid with “realistic” KF2 sequential esti- NhGi =−c. Nh mates of UECC in phase. (By “realistic” I refer to realistic ∂x OT clock diffusion coefficient values.) Gi Nh Estimation of x and x by the Kalman filter will reduce OR OT their error variances. 6.7. Observable error independent for each clock At each applicable time subtract the estimate of the UECC 6.5. Partition of KF1 estimation errors from the KF1 phase deviation estimate, for each particular Subtract estimated clock deviations from simulated (true) GPS clock, to estimate the OEIC in phase for that clock. Dur- clock deviations to define and quantify Kalman filter (KF1) ing measurement processing, the OEIC is contained within estimation errors. Adopt Brown’s additive partition of KF1 an envelope of a few parts of a nanosecond (see Figure 4). estimation errors into two components. I refer to the first Figure 4 presents a graph of two cases of the OEIC component as the unobservable error common to each clock for ground station clock S1. For the blue line of Figure 4, 6 International Journal of Navigation and Observation KF1 phase errors & KF2 UECC estimate (s) Simulated & estimated phase deviates (s) Case 12 Case 12 2-state clocks 2-state clocks Diffusion coefficient values 2.4e−7 S1: 3.23e-12 0e−17 −2e−8 Diffusion coefficient values S2: 1.12e-11 5.3e−17 S1: 3.23e−12 0e−17 N 1: 4e-11 1e−16 1.9e−7 Estimated N 2: 3e-11 8e−17 N 1 S2: 1.12e−11 5.3e−17 −7e−0 N 1: 4e−11 1e−16 S1 N 2 1.4e−7 S2 N 2: 3e−11 8e−17 KF1 R = (1 cm) 2 for GPS pseudo-range measurements −1.2e−7 time granularity = 30 s, when visible 9e−8 KF1 initialization errors −1.7e−7 N 1 Simulated 4e−8 N 2 S1 KF1 S1, S2, N 1, N 2 phase error ensemble −2.2e−7 −1e−8 S2 overlaid by UECC estimate in light blue 8days 100000 300000 500000 700000 100000 300000 500000 700000 Time (s) Time (s) Figure 3: KF1 phase errors and KF2 UECC estimate. Figure 2: Simulated and estimated phase deviations for four 2-state clocks. 7.3. Case 12 intervals of link visibility and KF1 range measurement pro- cessing are clearly distinguished from propagation intervals The calculation of S = σ /σ , α ∈{1, 2, 3, 4}, according to α 2α 1α with no measurements. During measurement processing, the the diffusion coefficient values presented in Figure 1 shows observable component of KF1 estimation error is contained that PPN is not satisfied for Case 12: within an envelope of a few parts of a nanosecond. 2S1 −6 −1 Calculation of the sequential covariance for the OEIC re- = 0.00 × 10 s , 1S1 quires a matrix value for the cross-covariance between the 2S2 −6 −1 KF1 phase deviation estimation error and the UECC estima- = 4.73 × 10 s , tion errorateachtime. Ihavenot yetbeenabletocalculate 1S2 (38) this cross-covariance. σ 2N 1 −6 −1 = 2.50 × 10 s , 1N 1 7. ALLAN VARIANCE AND PPN RELATIONS 2N 2 −6 −1 = 2.67 × 10 s . 1N 2 7.1. Allan coefficients versus diffusion coefficients 8. KALMAN FILTERS KF1 AND KF2 Denote τ as clock averaging time, σ (τ ) as Allan variance, a as Allan’s FMWN coefficient, a as Allan’s FMRW coef- 0 −2 I have simulated GPS pseudorange measurements for two ficient, σ as the FMWN diffusion coefficient, and σ as the 1 2 GPS ground station clocks S1 and S2, and for two GPS FMRW diffusion coefficient. Then, NAVSTAR clocks N1 and N2. Here I set simulated measure- 2 −1 2 −1 2 ment time granularity to 30 s for the set of all visible link in- σ (τ ) = a τ + a τ = σ τ + σ τ , (34) 0 −2 y 1 2 tervals. Visible and non-visible intervals are clearly evident in the blue line of Figure 4. I set the scalar root-variance R where for both measurement simulations and Kalman filter KF1 to σ = a , √ √ 1 0 R = 1 cm. Typically R∼1 m for GPS pseudorange, but (35) when carrier phase measurements are processed simultane- σ = 3a . 2 −2 ously with pseudorange, the root-variance is reduced by two orders of magnitude. So the use of R = 1cm enables me 7.2. Proportionate process noise (PPN) to quantify lower performance bounds for the simultaneous processing of both measurement types. Let α denote a variable α ∈{1, 2, 3,... , N} to identify each GPS clock in an ensemble of N clocks. For each clock α define 8.1. Create GPS clock ensemble the ratio S between diffusion coefficients σ and σ : α 1α 2α 2α Typically, one processes measurements with a Kalman filter S = , σ > 0. (36) α 1α 1α to derive sequential estimates of a multidimensional observ- able state. Instead, here I imitate the GPS operational pro- Then PPN is defined when, for each GPS clock α and each cedure and process simulated GPS pseudorange measure- associated ratio S ,wehave ments with KF1 to create asequenceofunobservablemul- S = S = S = ··· = S . (37) 1 2 3 N tidimensional clock state estimates. Clock state components James R. Wright 7 S1 observable component phase errors (s) each of the real GPS clocks, then quantitative upper bounds Case 12 can be calculated on OEIC magnitudes. These calculations 2-state clocks require the use of a rigorous simulator. Existence of signifi- 3e−9 KF2 process noise covariance cant cross-correlations between GPS clock phase errors and ∧ ∧ Color Q(phase ( s 2)) Q(freq ( s/ s) 2) other nonclock GPS estimation modeling errors enables sig- 2e−9 BLue 1e−1 1e−2 nificant aliasing into GPS clock phase estimates during op- Red 1e−16 1e−26 1e−9 eration of KF1 on real data. But given rigorous quantita- tive upper bounds on OEIC magnitudes, then significant vi- 0e−0 olation of these bounds when processing real GPS pseudo- range and carrier phase data identifies nonclock modeling −1e−9 errors related to the GPS estimation model. Modeling error (S1 observable component phase) = (S1 phase error) - (UECC estimate) candidates here include NAVSTAR orbit force modeling er- −2e−9 KF1 KF2 rors, ground antenna modeling errors (multipath), and tro- Sim GPS pseudo-range KF1 measurement granularity = 30 s Largest spikes are where KF1 measurements are not visible pospheric modeling errors. NAVSTAR orbit force modeling −3e−9 errors include those of solar photon pressure, albedo, ther- 100000 300000 500000 700000 mal dump, and propellant outgassing. The accuracy of this Time (s) diagnostic tool depends on the use of realistic clock diffusion coefficient values and a rigorous clock model simulation ca- Figure 4: S1 observable component phase error. pability. 10. OBSERVABLE CLOCKS are unobservable from GPS pseudorange measurements. See Figure 2 for an example of an ensemble of estimated unob- In an earlier version of my paper, I reported on KF1 valida- servable clock phase deviation state components created by tion results where clock S1 was specified as a TAI/UTC clock, KF1. external to the GPS clock ensemble consisting of S2, N1, and N2. This brought observablity (see Sections 5 and 6 herein) 8.1.1. Sherman’s theorem to S2, N1, and N2 clock states from GPS pseudorange mea- surements, drove clocks S2, N1, and N2 immediately to the GPS time, the unobservable GPS clock ensemble mean TAI/UTC timescale, and enabled a clean validation of my fil- phase, is created by the use of Sherman’s theorem [11, 18] ter implementation. Also it raised the question: why not the in the USAF Kalman filter measurement update algorithm same thing for the real GPS clock ensemble? Discussions with on GPS range measurements. Satisfaction of Sherman’s The- Ed Powers (USNO) and Bill Feess (Aerospace Corporation) orem guarantees that the mean-squared state estimate error reveal that this approach was tried and discarded after the on each observable state estimate component is minimized. difficulty in recovery from an uplink hardware failure was But the mean-squared state estimate error on each unobserv- blamed on the use of a single TAI/UTC Master Clock. This able state estimate component is not reduced. Thus the un- issue was resolved with Kenneth Brown’s introduction of the observable clock phase deviation state estimate component implicit ensemble mean. The mean phase (GPS time) of the common to every GPS clock is isolated by application of GPS clock ensemble will remain unobservable to GPS pseu- Sherman’s theorem. An ensemble of unobservable state esti- dorange measurements in the USAF Kalman filter for the mate components is thus created by Sherman’s theorem—see foreseeable future. Figure 3 for an example. REFERENCES 8.2. Initial condition errors [1] T. W. Anderson, “The integral of a symmetric unimodal func- A significant result emerges due to the modeling of Kalman tion over a symmetric convex set and some probability in- filter (KF1) initial condition errors in phase and frequency. equalities,” Proceedings of the American Mathematical Society, Initial estimated clock phase deviations are significantly dis- vol. 6, no. 2, pp. 170–176, 1955. placed by the KF1 initial condition errors in phase. As time [2] K. R. Brown, “The theory of the GPS composite clock,” in Pro- evolves estimated clock phase deviation magnitudes diverge ceedings of the 4th International Technical Meeting of the Satel- lite Division of the Institute of Navigation (ION GPS ’91),pp. continuously and increasingly when referred to true (simu- 223–241, Albuquerque, NM, USA, September 1991. lated) phase deviations, and this is due to filter initial condi- [3] R. S. Bucy and P. D. Joseph, Filtering for Stochastic Processes tion errors in frequency. See Figure 2 for an example. with Applications to Guidance, Interscience, New York, NY, USA, 1968. 9. IDENTIFY NONCLOCK MODELING ERRORS [4] W. Feess, “The Aerospace Corporation,” private Communica- tions, 2006. My interest in the GPS NAVSTAR (SV) orbit determination [5] R. J. Gardner, “The Brunn-Minkowski inequality,” Bulletin of problem, combined with that of the clock parameter estima- the American Mathematical Society, vol. 39, no. 3, pp. 355–405, tion problem, has enabled the identification of a useful diag- nostic tool: given realistic values for diffusion coefficients for [6] C. A. Greenhall, private Communications, 2006-2007. 8 International Journal of Navigation and Observation [7] C. A. Greenhall, “A Kalman filter clock ensemble algorithm that admits measurement noise,” Metrologia,vol. 43, no.4,pp. S311–S321, 2006. [8] S. T. Hutsell, “Relating the hadamard variance to MCS Kalman filter clock estimation,” in Proceedings of the 27th Annual Pre- cise Time and Time Interval (PTTI) Applications and Planning Meeting, p. 293, San Diego, Calif, USA, December 1995. [9] R. E. Kalman, “New methods in wiener filtering theory,” in Proceedings of the 1st Symposium on Engineering Applications of Random Function Theory and Probability,J.L.Bogdanoff andF.Kozin, Eds.,JohnWiley &Sons, NewYork, NY,USA, [10] D. Matsakis, private Communication, November 2007. [11] J. S. Meditch, Stochastic Optimal Linear Estimation and Con- trol, McGraw-Hill, New York, NY, USA, 1969. [12] O. J. Oaks, T. B. McCaskill, M. M. Largay, W. G. Reid, and J. A. Buisson, “Performance of GPS on-orbit NAVSTAR frequency standards and monitor station time references,” in Proceedings of the 30th Annual Precise Time and Time Interval (PTTI) Meet- ing, pp. 135–143, Reston, Va, USA, December 1998. [13] E. Powers, private Communications, 2006. [14] W. Riley, private Communications, PTTI meeting, December [15] S. Sherman, “A theorem on convex sets with applications,” The Annals of Mathematical Statistics, vol. 26, no. 4, pp. 763–767, [16] S. Sherman, “Non-mean-square error criteria,” IEEE Transac- tions on Information Theory, vol. 4, no. 3, pp. 125–126, 1958. [17] E. M. Stein and R. Shakarchi, Real Analysis, Princeton Univer- sity Press, Princeton, NJ, USA, 2005. [18] J. R. Wright, “Sherman’s theorem,” in The Malcolm D. Shuster Astronautics Symposium (AAS ’05), Grand Island, NY, USA, June 2005. [19] C. Zucca and P. Tavella, “The clock model and its relationship with the allan and related variances,” IEEE Transactions on Ul- trasonics, Ferroelectrics, and Frequency Control,vol. 52, no.2, pp. 289–295, 2005. [20] J. R. Wright, “GPS composite clock analysis,” in IEEE Interna- tional Frequency Control Symposium, European Frequency and Time Forum, pp. 523–528, Geneva, Switzerland, June 2007. International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Navigation and Observation Hindawi Publishing Corporation

Loading next page...
 
/lp/hindawi-publishing-corporation/gps-composite-clock-analysis-fJ3PGioBTc
Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2008 James R. Wright. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ISSN
1687-5990
DOI
10.1155/2008/261384
Publisher site
See Article on Publisher Site

Abstract

Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2008, Article ID 261384, 8 pages doi:10.1155/2008/261384 Research Article James R. Wright Analytical Graphics, Inc., 220 Valley Creek Blvd, Exton, PA 19341, USA Correspondence should be addressed to James R. Wright, jwright@agi.com Received 30 June 2007; Accepted 6 November 2007 Recommended by Demetrios Matsakis 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. Copyright © 2008 James R. Wright. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1. INTRODUCTION The estimation of NAVSTAR orbits would be incomplete without the simultaneous estimation of GPS clock param- GPS time is created by processing GPS pseudorange measure- eters. I use simulated GPS clock phase and frequency devi- ments with the operational GPS Kalman filter. Brown [2] ations, and simulated GPS pseudorange measurements, to refers to the object created by the Kalman filter as the GPS study Kalman filter estimation errors. composite clock, and to GPS time as the implicit ensemble This paper was first prepared for TimeNav’07 [20]. I am mean phase of the GPS composite clock. The fundamental indebted to Charles Greenhall (JPL) for encouragement and goal by the USAF and the USNO is to control GPS time to help in this work. within a specified bound of UTC/TAI. (I refer to TAI/UTC understanding that UTC has an accumulated discontinuity 2. THE COMPLETE ESTIMATION AND (a sum of leap seconds) when compared to TAI. But unique CONTROL PROBLEM two-way transformations between TAI and UTC have been in successful operational use since 1972. I have no need herein The USNO operates two UTC/TAI master clocks, each of to further distinguish between TAI and UTC.) I present here whichprovidesaccess to an estimate of UTC/TAIinrealtime a quantitative analysis of the GPS composite clock, derived (1 pps). One of these clocks is maintained at the USNO, and from detailed simulations and associated graphics. GPS clock the other is maintained at Schriever Air Force Base in Col- diffusion coefficient values used here were derived from Al- orado Springs. This enables the USNO to compare UTC/TAI lan deviation graphs presented by Oaks et al. [12] in 1998. I to the phase of each GPS orbital NAVSTAR clock via GPS refer to them as “realistic,” and in the sequel I claim “realistic” pseudorange measurements, by using a UTC/TAI master results from their use. Figure 1 presents their diffusion coef- clockinaUSNO GPSgroundreceiver. Each GPSclock is ficient values and my derivation of associated Allan deviation a member of (internal to) the GPS ensemble of clocks, but lines. the USNO master clock is external to the GPS ensemble of My interest in the GPS composite clock derives from clocks. Because of this, the difference between UTC/TAI and my interest in performing real-time orbit determination for the phase of each NAVSTAR GPS clock is observable. This GPS NAVSTAR spacecraft from ground receiver pseudor- difference can be (and is) estimated and quantified. The root ange measurements. (James R Wright is the architect of mean square (RMS) on these differences quantifies the differ- ODTK (Orbit Determination Tool Kit), a commercial soft- ence between UTC/TAI and GPS time. Inspection of the dif- ware product offered by Analytical Graphics, Inc. (AGI).) ferences between UTC/TAI and the phase of each NAVSTAR 2 International Journal of Navigation and Observation Case 12 The sequential estimation of GPS clock deviations re- 1e−11 quires the development of a linear TU and nonlinear MU. 8 N 1 2-state clocks N 2 Diffusion coefficient values The nonlinear MU must be linearized locally to enable ap- √ √ ID σ1( s) σ2(1/ s) plication of the linear Kalman MU. Kalman’s MU [9]de- S13.23e−12 S2 S21.12e − 11 5.3e−17 rives from Sherman’s theorem [11, 15, 16], Sherman’s the- N14e − 11 1e−16 1e−12 N23e − 11 8e−17 orem derives from Anderson’s theorem [1], and Anderson’s S1 theorem derives from the Brunn-Minkowki inequality the- 2 orem [5, 17]. The theoretical foundation for my linearized Time durations T = 15552000 s MU derives from these theorems. 1e−13 τ = 32 s τ 0 = 8s 4.1. Initial conditions PRN seed number S11 1e−14 S2 30, 000, 001 Initialization of all sequential estimators requires the use of 90% ensemble last calc. 6 N 1 60, 000, 001 N 2 90, 000, 001 an initial state estimate column matrix X and an intial state 0|0 estimate error covariance matrix P for time t . 0|0 0 1 2 3 4 5 6 10 10 10 10 10 10 τ (s) 4.2. Kalman filter: linear TU and linear MU Figure 1: Allan deviation lines for S1, S2, N1, and N2. Derivation and calculation for the discrete-time Kalman fil- ter, linear in both TU and MU, is best presented by Meditch [11, Chapter 5]. GPS clock enables the USNO to identify GPS clocks that re- quire particular frequency-rate control corrections. Use of this knowledge enables the USAF to adjust frequency rates of 4.3. Linear TU and nonlinear MU selected GPS clocks. Currently, the USAF uses an automated bang-bang controller on frequency-rate. (According to Bill The simultaneous sequential estimation of GPS clock phase and frequency deviation parameters can be studied with the Feess, an improvement in control can be achieved by replac- development of a linear TU and nonlinear MU for the clock ing the existing “bang-bang controller” with a “proportional controller.”) state estimate subset. This is useful to study clock parameter estimation, as demonstrated in Section 6. Let X denote an n × 1 column matrix of state estimate j|i 3. STOCHASTIC CLOCK PHYSICS components, where the left subscript j denotes state epoch t and the right subscript i denotes time-tag t for the last ob- The most significant stochastic clock physics are under- i servation processed, where i, j ∈{0, 1, 2,...}.Let P denote stood in terms of Wiener processes and their integrals. j|i an associated n× n square symmetric state estimate error co- Clock physics are characterized by particular values of clock- variance matrix (positive eigenvalues). dependent diffusion coefficients, and are conveniently stud- ied with aid of a relevant clock model that relates diffusion coefficient values to their underlying Wiener processes. For 4.3.1. Linear TU my presentation here Ihaveselected“Theclock modeland its relationship with the Allan and related variances” pre- For k ∈{0, 1, 2, 3,... , M}, the propagation of the true un- sented as an IEEE paper by Zucca and Tavella [19] in 2005. known n × 1 matrix state X is given by Except for FM flicker noise, this model captures the most significant physics for all GPS clocks. I simulate and vali- X = Φ X + J ,(1) k+1 k+1,k k k+1,k date GPS pseudorange measurements using simulated phase deviations and simulated frequency deviations, according to where J is called the process noise matrix. Propagation of k+1,k Zucca and Tavella. the known n × 1 matrix state estimate X is given by k|k 4. KALMAN FILTERS X = Φ X (2) k+1|k k+1,k k|k I present my approach for the optimal sequential estimation because the conditional mean of J is zero. Propagation k+1,k of clock deviation states and their error covariance functions. of the known n × n matrix state estimate error covariance Sequential state estimates are generated recursively from two matrix P is given by k|k multidimensional stochastic update functions, the time up- date (TU) and the measurement update (MU). The TU moves P = Φ P Φ + Q ,(3) the state estimate and covariance forward with time, accu- k+1|k k+1,k k|k k+1,k k+1,k mulating integrals of random clock deviation process noise in the covariance. The MU is performed at a fixed measure- where the n × n matrix Q is called the process noise co- k+1,k ment time where the state estimate and covariance are cor- variance matrix (see [19] for concrete clock examples of J k+1,k rected with new observation information. and Q ). k+1,k James R. Wright 3 4.3.2. Nonlinear MU 4.5. Kalman filter advantage Calculate the n × 1matrixfilter gain K : Severe computational problems are incurred in any attempt k+1 to estimate unobservable states using iterated batch least −1 T T K = P H H P H + R . (4) k+1 k+1|k k+1 k+1|k k+1 k+1 k+1 squares methods or iterated maximum likelihood methods for navigation, because state-sized inversions of singular ma- The filter measurement update state estimate n × 1matrix trices are required. Here the Kalman filter is distinguished in X , due to the observation y , is calculated with k+1|k+1 k+1 that estimates of unobservable states can be created and used without matrix inversion problems because the Kalman filter X = X + K y − y X ,(5) k+1|k+1 k+1|k k+1 k+1 k+1|k MU is free of state-sized matrix inversions. where R is the scalar variance on the observation resid- k+1 By design, one typically estimates observable states. But ual y − y(X ), and y(X ) is a nonlinear function of k+1 k+1|k k+1|k the Kalman filter enables one to create unobservable states. X . Define the error ΔX in X : k+1|k k+1|k+1 k+1|k+1 The USAF chose to create unobservable GPS clock parameter states for construction of GPS time. ΔX = X − X . (6) k+1|k+1 k+1 k+1|k+1 Define the n × n state estimate error covariance matrix 5. OBSERVABILITY P with k+1|k+1 Ihavedefined observability in terms of a Kalman filter for- P = E ΔX ΔX . (7) k+1|k+1 k+1|k+1 k+1|k+1 mulation, and I have proved simple theorems related thereto. My definition of observability is different than Kalman’s def- Bucy and Joseph [3, page 141] recommend that P k+1|k+1 inition and, unlike Kalman’s definition, is directly applicable should be calculated with to covariance matrices derived from a Kalman filter. P = P − T,(8) k+1|k+1 k+1|k 5.1. Definition where T −1 If the state estimate error variance of a particular state es- T = P H R H P , k+1|k k+1 k+1|k k+1 k+1 (9) timate component is reduced by processing an observation, R = H P H + R . k+1 k+1 k+1|k k+1 k+1 then that state estimate component is observable to that ob- servation. Otherwise, that state estimate component is not Equations (8)and (9) reduce to the form given by Kalman: observable (unobservable) to that observation. P = I − K H P . (10) k+1|k+1 k+1 k+1 k+1|k Theorem 1. If every component of the row matrix H of k+1 Calculation of P by (8)and (9) is numerically stable, k+1|k+1 measurement-state partial derivatives is zero at time t , then k+1 whereas the Kalman form calculation is not. every component of the state estimate X is unobservable at k+1 time t . k+1 4.4. Nonlinear TU and nonlinear MU Proof. H = 0 implies that P = P according to k+1 k+1|k+1 k+1|k Refer to Section 4.3.2 for the nonlinear MU. (10). Thus none of the variances of P are reduced due to k+1|k processing the observation y . Then by definition, X is k+1 k+1 4.4.1. Nonlinear TU unobservable in every component. The nonlinear TU always spans a nonempty time interval Theorem 2. Given values for scalars H , P > 0, R > k+1 k+1|k k+1 and requires the use of a numerical state estimate integra- 0 at time t , and given that H = / 0, then the scalar state k+1 k+1 tor ϕ . Given an initial time t ,afinal time t ,and aforce 0 f estimate X is observable at time t . k+1 k+1 model u(X (t), t), then ϕ propagates the state estimate X (t ) from t to t using forces u(X (τ ), τ)toget X (t ). That is, 0 f f 2 Proof. The obvious inequality P H + R > k+1|k k+1 k+1 P H > 0 implies that k+1|k k+1 X t = ϕ t ; X t , t , u X (τ ), τ , t ≤ τ ≤ t . (11) f f 0 0 0 f This can be shortened to write P H k+1|k k+1 1 > > 0. (14) P H + R k+1|k k+1 k+1 X t = ϕ t ; X t , t , (12) f f 0 0 Multiply through by −1: where the use of forces u(X (τ ), τ ) is tacitly implied. Thus, ϕ is a column matrix with n elements: 2 P H k+1|k k+1 ⎡ ⎤ −1 < − < 0. (15) P H + R x1 k+1|k k+1 k+1 ⎢ ⎥ ⎢ x2⎥ ⎢ ⎥ Add 1: ⎢ ⎥ x3 ϕ = . (13) ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ P H . k+1|k k+1 0 < 1 − < 1. (16) ϕ P H + R xn k+1|k k+1 k+1 4 International Journal of Navigation and Observation Multiply through by P : in the operational algorithms, but state estimate observabil- k+1|k ity is independent of relativity, so observability can be de- Nh P H k+1|k k+1 fined and discussed independent of relativity.) Let δx and 0 < 1 − P <P . (17) k+1|k k+1|k Gi P H + R k+1|k k+1 δx denote Kalman filter estimation errors in clock phase k+1 Nh Gi Dh for t and t . Define time of transmission difference t T R T Now use (4)and (10)towrite Di andtimeofreceipt difference t : P H k+1|k k+1 Dh Nh Nh P = 1 − P . (18) k+1|k+1 k+1|k 2 t = t − δx , T T T P H + R k+1|k k+1 k+1 (20) Di Gi Gi t = t − δx . Insert (18) into the inequality (17) to get the result: R R R 0 <P <P . (19) Thus, k+1|k+1 k+1|k Nh Dh Nh Thus the variance P is reduced due to processing the ob- k+1|k t = t + δx , T T T servation y . Then the scalar state X is observable by def- (21) k+1 k+1 Gi Di Gi t = t + δx . inition. R R R Nh Gi Equation (21)present t and t as additive combinations of T R 5.2. Theoretical foundation Dh Di deterministic times t and t and Kalman filter estimation T R Nh Gi These theorems are referred to expressions given by Kalman errors in clock phase δx and δx . Define the one-way GPS T R pseudorange measurement ρ : for filter gain K and covariance P ,see (4)and (10). k+1 k+1|k+1 NhGi Kalman’s expressions are derived from the rigorous theo- Gi Nh Gi Nh rem chain provided by Sherman, Anderson, and Brunn- ρ = c t − t , t >t . (22) NhGi R T R T Minkowski—the theoretical foundation is deep. Insert (21) into (22): 5.3. Determine observability directly Di Gi Dh Nh ρ = c t + δx − t + δx NhGi R R T T (23) Given an optimal sequential estimator, given a particular col- Di Dh Gi Nh = c t − t + δx − δx , R T R T lection of applicable observations (real or simulated), and given realistic state estimate error covariance matrices P k+1|k where c is speed of light in vacuum. Define and P at each time t , apply the definition of observ- k+1|k+1 k+1 ability directly (note that this is impossible using Kalman’s Di Dh Δt = t − t , (24) R T definition of observability) to distinguish between observable hi Gi Nh and unobservable state elements. An optimal sequential es- δt = δx − δx . (25) R T timator is designed to eliminate significant aliasing between Then, estimated state elements, and thus enables this distinction. hi ρ = c Δt + δt , (26) NhGi 6. UNOBSERVABLE GPS CLOCK STATES where Δt is deterministic and δt is random. 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 6.2. Partition of Kalman filter estimation errors the clock phase of every operational GPS clock is unobserv- Let x denote the phase component of Kalman filter estima- able from GPS pseudorange observations, as demonstrated tion error that is common to every GPS ensemble clock, when below. GPS NAVSTAR orbit parameters are observable from Gi Nh it exists. Define phase differences x and x with OR OT GPS pseudorange observations. The USAF Kalman filter si- multaneously estimates orbit parameters and clock param- Gi Gi x = δx − x , OR R eters from GPS pseudorange observations, so the state esti- (27) mate is partitioned in this manner into a subset of unobserv- Nh Nh x = δx − x OT T able clock parameters and a subset of observable orbit pa- rameters. This partition is performed by application of Sher- for ground station i and NAVSTAR h. Then Kalman filter es- man’s theorem in the MU. Gi timation errors δx , i ∈{1, 2,...}, for ground station clocks Nh and δx , h ∈{1, 2,...}, for NAVSTAR clocks have the ad- 6.1. GPS pseudorange representation ditive partition: Nh Gi Gi Let t denote time of radio wave transmission for the hth δx = x + x , R C OR Gi NAVSTAR clock, and let t denote time of radio wave receipt (28) Nh Nh for the ith ground station clock. (Refer all times to a coordi- δx = x + x . T OT nate time, e.g., to GPS time. Appropriate transformations be- tween proper time and coordinate time must be performed (This partition was introduced by Brown [2].) James R. Wright 5 6.3. The common random phase component (UECC), and to the second component as the observable is unobservable error independent for each clock (OEIC). (Observability is meaningful here only when processing simulated GPS pseu- Gi Insert (28) into (25): dorange data.) See (28). x is the UECC, x is the OEIC OR Nh for ground station clocks and x is the OEIC for NAVS- OT hi Gi Nh δt = δx − δx R T TAR clocks. On initialization of KF1, the variances on the Gi Nh UECC and OEIC are identical. On processing the first GPS = x + x − x + x (29) C C OR OT pseudorange measurements with KF1 the variances on both Gi Nh = x − x . fall quickly. But with continued measurement processing the OR OT variances on the UECC increase without bound while the Insert (29) into (26): variances on the OEIC appoach zero asymptotically. For simulated GPS pseudorange data I create an optimal Gi Nh ρ = c Δt + x − x . (30) OR OT NhGi sequential estimate of the UECC by application of a second Kalman filter KF2 to pseudomeasurements defined by the Thus, the random phase component x that is common to phase components of KF1 estimation errors. the Kalman filter estimation error for every ensemble clock Since there is no physical process noise on the UECC, an has vanished in the range representation ρ .Variations NhGi estimate of the UECC can also be achieved using a batch least Δx in x cannot cause variations Δρ in ρ : C C NhGi NhGi squares estimation algorithm on the phase components of KF1 estimation errors—demonstrated previously by Green- ∂ρ NhGi Δρ = Δx (31) hall [7]. (I apply sufficient process noise covariance for KF2 NhGi ∂x to mask the effects of double-precision computer word trun- cation. Without this, KF2 does diverge.) because the partial derivative H = ∂ρ /∂x is zero: NhGi ∂ρ NhGi 6.6. Unobservable error common to each clock = 0. (32) ∂x There are at least four techniques to estimate the UECC An application of Theorem 1 to (32) demonstrates that x is C when simulating GPS pseudorange data. First, one could take unobservable from ρ . the sample mean of KF1 estimation errors across the clock NhGi But the architect who designs the complete estimator ensemble at each time and form a sample variance about must design an optimal NAVSTAR orbit estimator to prevent the mean; this would yield a sequential sampling procedure, aliasing from NAVSTAR orbit estimation errors into x .It C but where each mean and variance is sequentially uncon- helps to know that there is no coupling between the orbit and nected. Second, one can employ Ken Brown’s implicit ensem- x in the complete state transition function. I have provided C ble mean (IEM) and covariance; this is a batch procedure re- a new method herein to identify this aliasing, and I have pro- quiring an inversion of the KF1 covariance matrix followed vided suggestions on where to look for inadequate modeling by a second matrix inversion of the modified covariance ma- that would be the source of this aliasing. See Section 9. trix 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, 6.4. Independent random phase components and are processed by a batch least squares estimator to ob- are observable tain optimal batch estimates and covariance matrices for the Gi Nh The independent phase deviations x and x are observ- UECC. Fourth, one can treat the KF1 phase estimation er- OR OT able to GPS pseudorange observations because their partial rors as pseudomeasurements, invoke a second Kalman filter derivatives are nonzero: (KF2), and process these phase pseudomeasurements with KF2 to obtain optimal sequential estimates and variances ∂ρ NhGi = +c, for the UECC. I have been successful with this approach. Gi ∂x OR Figure 3 presents an ensemble of “realistic” KF1 phase esti- (33) ∂ρ mation errors, overlaid with “realistic” KF2 sequential esti- NhGi =−c. Nh mates of UECC in phase. (By “realistic” I refer to realistic ∂x OT clock diffusion coefficient values.) Gi Nh Estimation of x and x by the Kalman filter will reduce OR OT their error variances. 6.7. Observable error independent for each clock At each applicable time subtract the estimate of the UECC 6.5. Partition of KF1 estimation errors from the KF1 phase deviation estimate, for each particular Subtract estimated clock deviations from simulated (true) GPS clock, to estimate the OEIC in phase for that clock. Dur- clock deviations to define and quantify Kalman filter (KF1) ing measurement processing, the OEIC is contained within estimation errors. Adopt Brown’s additive partition of KF1 an envelope of a few parts of a nanosecond (see Figure 4). estimation errors into two components. I refer to the first Figure 4 presents a graph of two cases of the OEIC component as the unobservable error common to each clock for ground station clock S1. For the blue line of Figure 4, 6 International Journal of Navigation and Observation KF1 phase errors & KF2 UECC estimate (s) Simulated & estimated phase deviates (s) Case 12 Case 12 2-state clocks 2-state clocks Diffusion coefficient values 2.4e−7 S1: 3.23e-12 0e−17 −2e−8 Diffusion coefficient values S2: 1.12e-11 5.3e−17 S1: 3.23e−12 0e−17 N 1: 4e-11 1e−16 1.9e−7 Estimated N 2: 3e-11 8e−17 N 1 S2: 1.12e−11 5.3e−17 −7e−0 N 1: 4e−11 1e−16 S1 N 2 1.4e−7 S2 N 2: 3e−11 8e−17 KF1 R = (1 cm) 2 for GPS pseudo-range measurements −1.2e−7 time granularity = 30 s, when visible 9e−8 KF1 initialization errors −1.7e−7 N 1 Simulated 4e−8 N 2 S1 KF1 S1, S2, N 1, N 2 phase error ensemble −2.2e−7 −1e−8 S2 overlaid by UECC estimate in light blue 8days 100000 300000 500000 700000 100000 300000 500000 700000 Time (s) Time (s) Figure 3: KF1 phase errors and KF2 UECC estimate. Figure 2: Simulated and estimated phase deviations for four 2-state clocks. 7.3. Case 12 intervals of link visibility and KF1 range measurement pro- cessing are clearly distinguished from propagation intervals The calculation of S = σ /σ , α ∈{1, 2, 3, 4}, according to α 2α 1α with no measurements. During measurement processing, the the diffusion coefficient values presented in Figure 1 shows observable component of KF1 estimation error is contained that PPN is not satisfied for Case 12: within an envelope of a few parts of a nanosecond. 2S1 −6 −1 Calculation of the sequential covariance for the OEIC re- = 0.00 × 10 s , 1S1 quires a matrix value for the cross-covariance between the 2S2 −6 −1 KF1 phase deviation estimation error and the UECC estima- = 4.73 × 10 s , tion errorateachtime. Ihavenot yetbeenabletocalculate 1S2 (38) this cross-covariance. σ 2N 1 −6 −1 = 2.50 × 10 s , 1N 1 7. ALLAN VARIANCE AND PPN RELATIONS 2N 2 −6 −1 = 2.67 × 10 s . 1N 2 7.1. Allan coefficients versus diffusion coefficients 8. KALMAN FILTERS KF1 AND KF2 Denote τ as clock averaging time, σ (τ ) as Allan variance, a as Allan’s FMWN coefficient, a as Allan’s FMRW coef- 0 −2 I have simulated GPS pseudorange measurements for two ficient, σ as the FMWN diffusion coefficient, and σ as the 1 2 GPS ground station clocks S1 and S2, and for two GPS FMRW diffusion coefficient. Then, NAVSTAR clocks N1 and N2. Here I set simulated measure- 2 −1 2 −1 2 ment time granularity to 30 s for the set of all visible link in- σ (τ ) = a τ + a τ = σ τ + σ τ , (34) 0 −2 y 1 2 tervals. Visible and non-visible intervals are clearly evident in the blue line of Figure 4. I set the scalar root-variance R where for both measurement simulations and Kalman filter KF1 to σ = a , √ √ 1 0 R = 1 cm. Typically R∼1 m for GPS pseudorange, but (35) when carrier phase measurements are processed simultane- σ = 3a . 2 −2 ously with pseudorange, the root-variance is reduced by two orders of magnitude. So the use of R = 1cm enables me 7.2. Proportionate process noise (PPN) to quantify lower performance bounds for the simultaneous processing of both measurement types. Let α denote a variable α ∈{1, 2, 3,... , N} to identify each GPS clock in an ensemble of N clocks. For each clock α define 8.1. Create GPS clock ensemble the ratio S between diffusion coefficients σ and σ : α 1α 2α 2α Typically, one processes measurements with a Kalman filter S = , σ > 0. (36) α 1α 1α to derive sequential estimates of a multidimensional observ- able state. Instead, here I imitate the GPS operational pro- Then PPN is defined when, for each GPS clock α and each cedure and process simulated GPS pseudorange measure- associated ratio S ,wehave ments with KF1 to create asequenceofunobservablemul- S = S = S = ··· = S . (37) 1 2 3 N tidimensional clock state estimates. Clock state components James R. Wright 7 S1 observable component phase errors (s) each of the real GPS clocks, then quantitative upper bounds Case 12 can be calculated on OEIC magnitudes. These calculations 2-state clocks require the use of a rigorous simulator. Existence of signifi- 3e−9 KF2 process noise covariance cant cross-correlations between GPS clock phase errors and ∧ ∧ Color Q(phase ( s 2)) Q(freq ( s/ s) 2) other nonclock GPS estimation modeling errors enables sig- 2e−9 BLue 1e−1 1e−2 nificant aliasing into GPS clock phase estimates during op- Red 1e−16 1e−26 1e−9 eration of KF1 on real data. But given rigorous quantita- tive upper bounds on OEIC magnitudes, then significant vi- 0e−0 olation of these bounds when processing real GPS pseudo- range and carrier phase data identifies nonclock modeling −1e−9 errors related to the GPS estimation model. Modeling error (S1 observable component phase) = (S1 phase error) - (UECC estimate) candidates here include NAVSTAR orbit force modeling er- −2e−9 KF1 KF2 rors, ground antenna modeling errors (multipath), and tro- Sim GPS pseudo-range KF1 measurement granularity = 30 s Largest spikes are where KF1 measurements are not visible pospheric modeling errors. NAVSTAR orbit force modeling −3e−9 errors include those of solar photon pressure, albedo, ther- 100000 300000 500000 700000 mal dump, and propellant outgassing. The accuracy of this Time (s) diagnostic tool depends on the use of realistic clock diffusion coefficient values and a rigorous clock model simulation ca- Figure 4: S1 observable component phase error. pability. 10. OBSERVABLE CLOCKS are unobservable from GPS pseudorange measurements. See Figure 2 for an example of an ensemble of estimated unob- In an earlier version of my paper, I reported on KF1 valida- servable clock phase deviation state components created by tion results where clock S1 was specified as a TAI/UTC clock, KF1. external to the GPS clock ensemble consisting of S2, N1, and N2. This brought observablity (see Sections 5 and 6 herein) 8.1.1. Sherman’s theorem to S2, N1, and N2 clock states from GPS pseudorange mea- surements, drove clocks S2, N1, and N2 immediately to the GPS time, the unobservable GPS clock ensemble mean TAI/UTC timescale, and enabled a clean validation of my fil- phase, is created by the use of Sherman’s theorem [11, 18] ter implementation. Also it raised the question: why not the in the USAF Kalman filter measurement update algorithm same thing for the real GPS clock ensemble? Discussions with on GPS range measurements. Satisfaction of Sherman’s The- Ed Powers (USNO) and Bill Feess (Aerospace Corporation) orem guarantees that the mean-squared state estimate error reveal that this approach was tried and discarded after the on each observable state estimate component is minimized. difficulty in recovery from an uplink hardware failure was But the mean-squared state estimate error on each unobserv- blamed on the use of a single TAI/UTC Master Clock. This able state estimate component is not reduced. Thus the un- issue was resolved with Kenneth Brown’s introduction of the observable clock phase deviation state estimate component implicit ensemble mean. The mean phase (GPS time) of the common to every GPS clock is isolated by application of GPS clock ensemble will remain unobservable to GPS pseu- Sherman’s theorem. An ensemble of unobservable state esti- dorange measurements in the USAF Kalman filter for the mate components is thus created by Sherman’s theorem—see foreseeable future. Figure 3 for an example. REFERENCES 8.2. Initial condition errors [1] T. W. Anderson, “The integral of a symmetric unimodal func- A significant result emerges due to the modeling of Kalman tion over a symmetric convex set and some probability in- filter (KF1) initial condition errors in phase and frequency. equalities,” Proceedings of the American Mathematical Society, Initial estimated clock phase deviations are significantly dis- vol. 6, no. 2, pp. 170–176, 1955. placed by the KF1 initial condition errors in phase. As time [2] K. R. Brown, “The theory of the GPS composite clock,” in Pro- evolves estimated clock phase deviation magnitudes diverge ceedings of the 4th International Technical Meeting of the Satel- lite Division of the Institute of Navigation (ION GPS ’91),pp. continuously and increasingly when referred to true (simu- 223–241, Albuquerque, NM, USA, September 1991. lated) phase deviations, and this is due to filter initial condi- [3] R. S. Bucy and P. D. Joseph, Filtering for Stochastic Processes tion errors in frequency. See Figure 2 for an example. with Applications to Guidance, Interscience, New York, NY, USA, 1968. 9. IDENTIFY NONCLOCK MODELING ERRORS [4] W. Feess, “The Aerospace Corporation,” private Communica- tions, 2006. My interest in the GPS NAVSTAR (SV) orbit determination [5] R. J. Gardner, “The Brunn-Minkowski inequality,” Bulletin of problem, combined with that of the clock parameter estima- the American Mathematical Society, vol. 39, no. 3, pp. 355–405, tion problem, has enabled the identification of a useful diag- nostic tool: given realistic values for diffusion coefficients for [6] C. A. Greenhall, private Communications, 2006-2007. 8 International Journal of Navigation and Observation [7] C. A. Greenhall, “A Kalman filter clock ensemble algorithm that admits measurement noise,” Metrologia,vol. 43, no.4,pp. S311–S321, 2006. [8] S. T. Hutsell, “Relating the hadamard variance to MCS Kalman filter clock estimation,” in Proceedings of the 27th Annual Pre- cise Time and Time Interval (PTTI) Applications and Planning Meeting, p. 293, San Diego, Calif, USA, December 1995. [9] R. E. Kalman, “New methods in wiener filtering theory,” in Proceedings of the 1st Symposium on Engineering Applications of Random Function Theory and Probability,J.L.Bogdanoff andF.Kozin, Eds.,JohnWiley &Sons, NewYork, NY,USA, [10] D. Matsakis, private Communication, November 2007. [11] J. S. Meditch, Stochastic Optimal Linear Estimation and Con- trol, McGraw-Hill, New York, NY, USA, 1969. [12] O. J. Oaks, T. B. McCaskill, M. M. Largay, W. G. Reid, and J. A. Buisson, “Performance of GPS on-orbit NAVSTAR frequency standards and monitor station time references,” in Proceedings of the 30th Annual Precise Time and Time Interval (PTTI) Meet- ing, pp. 135–143, Reston, Va, USA, December 1998. [13] E. Powers, private Communications, 2006. [14] W. Riley, private Communications, PTTI meeting, December [15] S. Sherman, “A theorem on convex sets with applications,” The Annals of Mathematical Statistics, vol. 26, no. 4, pp. 763–767, [16] S. Sherman, “Non-mean-square error criteria,” IEEE Transac- tions on Information Theory, vol. 4, no. 3, pp. 125–126, 1958. [17] E. M. Stein and R. Shakarchi, Real Analysis, Princeton Univer- sity Press, Princeton, NJ, USA, 2005. [18] J. R. Wright, “Sherman’s theorem,” in The Malcolm D. Shuster Astronautics Symposium (AAS ’05), Grand Island, NY, USA, June 2005. [19] C. Zucca and P. Tavella, “The clock model and its relationship with the allan and related variances,” IEEE Transactions on Ul- trasonics, Ferroelectrics, and Frequency Control,vol. 52, no.2, pp. 289–295, 2005. [20] J. R. Wright, “GPS composite clock analysis,” in IEEE Interna- tional Frequency Control Symposium, European Frequency and Time Forum, pp. 523–528, Geneva, Switzerland, June 2007. International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Journal

International Journal of Navigation and ObservationHindawi Publishing Corporation

Published: Jan 3, 2008

References