# GPS Composite Clock Analysis

Authors: James R. Wright
Journal: International Journal of Navigation and Observation
 (2008)

---

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


---

## Body

## 1. Introduction

GPS time iscreated 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.Figure 1
Allan deviation lines for S1, S2, N1, and N2.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 theoptimal 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 matrixX^0|0 and an intial
state estimate error covariance matrix P0|0 for time t0.


### 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 Section6.LetX^j|i denote an n×1 column matrix
of state estimate components, where the left subscript j denotes state
epoch tj and the right
subscript i denotes
time-tag ti for the last
observation processed, where i,j∈{0,1,2,…}. Let Pj|i denote an
associated n×n square
symmetric state estimate error covariance matrix (positive eigenvalues).


#### 4.3.1. Linear TU

Fork∈{0,1,2,3,…,M}, the propagation of the true unknown n×1 matrix state Xk is given by(1)Xk+1=Φk+1,kXk+Jk+1,k,where Jk+1,k is called the
process noise matrix. Propagation of the known n×1 matrix state
estimate X^k∣k is given by(2)X^k+1∣k=Φk+1,kX^k∣kbecause the conditional mean of Jk+1,k is zero.
Propagation of the known n×n matrix state
estimate error covariance matrix Pk|k is given by(3)Pk+1∣k=Φk+1,kPk∣kΦk+1,kT+Qk+1,k,where the n×n matrix Qk+1,k is called the
process noise covariance matrix (see [19] for concrete clock examples of Jk+1,k and Qk+1,k).


#### 4.3.2. Nonlinear MU

Calculate then×1 matrix filter
gain Kk+1:(4)Kk+1=Pk+1∣kHk+1T[Hk+1Pk+1∣kHk+1T+Rk+1]−1.
The filter measurement update state estimate n×1 matrix X^k+1|k+1, due to the observation yk+1, is calculated with(5)X^k+1∣k+1=X^k+1∣k+Kk+1[yk+1−y(X^k+1∣k)],where Rk+1 is the scalar
variance on the observation residual yk+1−y(X^k+1∣k), and y(X^k+1∣k) is a nonlinear
function of X^k+1∣k. Define the error ΔX^k+1∣k+1 in X^k+1∣k+1:(6)ΔX^k+1∣k+1=Xk+1−X^k+1∣k+1.Define the n×n state estimate
error covariance matrix Pk+1∣k+1 with (7)Pk+1∣k+1=E{(ΔX^k+1∣k+1)(ΔX^k+1∣k+1)T}.Bucy and Joseph [3, page 141] recommend
that Pk+1|k+1 should be
calculated with(8)Pk+1∣k+1=Pk+1∣k−T,where(9)T=Pk+1∣kHk+1TR˜k+1−1Hk+1Pk+1∣k,R˜k+1=Hk+1Pk+1∣kHk+1T+Rk+1.Equations (8) and (9) reduce to the form given by
Kalman:(10)Pk+1∣k+1=[I−Kk+1Hk+1]Pk+1∣k.Calculation of Pk+1∣k+1 by (8) and (9) is numerically stable, whereas the Kalman form
calculation is not.


### 4.4. Nonlinear TU and Nonlinear MU

Refer to Section4.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φx. Given an initial time t0, a final time tf, and a force model u(X^(t),t), then φx propagates the
state estimate X^(t0) from t0 to tf using forces u(X^(τ),τ) to get X^(tf). That is,(11)X^(tf)=φx{tf;X^(t0),t0,u(X^(τ),τ),t0≤τ≤tf}.This can be shortened to write(12)X^(tf)=φx{tf;X^(t0),t0},where the use of forces u(X^(τ),τ) is tacitly
implied. Thus, φx is a column
matrix with n elements:(13)φz=[φx1φx2φx3⋮φxn].


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


## 4.1. Initial Conditions

Initialization of all sequential estimators requires
the use of an initial state estimate column matrixX^0|0 and an intial
state estimate error covariance matrix P0|0 for time t0.


## 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 Section6.LetX^j|i denote an n×1 column matrix
of state estimate components, where the left subscript j denotes state
epoch tj and the right
subscript i denotes
time-tag ti for the last
observation processed, where i,j∈{0,1,2,…}. Let Pj|i denote an
associated n×n square
symmetric state estimate error covariance matrix (positive eigenvalues).


### 4.3.1. Linear TU

Fork∈{0,1,2,3,…,M}, the propagation of the true unknown n×1 matrix state Xk is given by(1)Xk+1=Φk+1,kXk+Jk+1,k,where Jk+1,k is called the
process noise matrix. Propagation of the known n×1 matrix state
estimate X^k∣k is given by(2)X^k+1∣k=Φk+1,kX^k∣kbecause the conditional mean of Jk+1,k is zero.
Propagation of the known n×n matrix state
estimate error covariance matrix Pk|k is given by(3)Pk+1∣k=Φk+1,kPk∣kΦk+1,kT+Qk+1,k,where the n×n matrix Qk+1,k is called the
process noise covariance matrix (see [19] for concrete clock examples of Jk+1,k and Qk+1,k).


### 4.3.2. Nonlinear MU

Calculate then×1 matrix filter
gain Kk+1:(4)Kk+1=Pk+1∣kHk+1T[Hk+1Pk+1∣kHk+1T+Rk+1]−1.
The filter measurement update state estimate n×1 matrix X^k+1|k+1, due to the observation yk+1, is calculated with(5)X^k+1∣k+1=X^k+1∣k+Kk+1[yk+1−y(X^k+1∣k)],where Rk+1 is the scalar
variance on the observation residual yk+1−y(X^k+1∣k), and y(X^k+1∣k) is a nonlinear
function of X^k+1∣k. Define the error ΔX^k+1∣k+1 in X^k+1∣k+1:(6)ΔX^k+1∣k+1=Xk+1−X^k+1∣k+1.Define the n×n state estimate
error covariance matrix Pk+1∣k+1 with (7)Pk+1∣k+1=E{(ΔX^k+1∣k+1)(ΔX^k+1∣k+1)T}.Bucy and Joseph [3, page 141] recommend
that Pk+1|k+1 should be
calculated with(8)Pk+1∣k+1=Pk+1∣k−T,where(9)T=Pk+1∣kHk+1TR˜k+1−1Hk+1Pk+1∣k,R˜k+1=Hk+1Pk+1∣kHk+1T+Rk+1.Equations (8) and (9) reduce to the form given by
Kalman:(10)Pk+1∣k+1=[I−Kk+1Hk+1]Pk+1∣k.Calculation of Pk+1∣k+1 by (8) and (9) is numerically stable, whereas the Kalman form
calculation is not.


## 4.3.1. Linear TU

Fork∈{0,1,2,3,…,M}, the propagation of the true unknown n×1 matrix state Xk is given by(1)Xk+1=Φk+1,kXk+Jk+1,k,where Jk+1,k is called the
process noise matrix. Propagation of the known n×1 matrix state
estimate X^k∣k is given by(2)X^k+1∣k=Φk+1,kX^k∣kbecause the conditional mean of Jk+1,k is zero.
Propagation of the known n×n matrix state
estimate error covariance matrix Pk|k is given by(3)Pk+1∣k=Φk+1,kPk∣kΦk+1,kT+Qk+1,k,where the n×n matrix Qk+1,k is called the
process noise covariance matrix (see [19] for concrete clock examples of Jk+1,k and Qk+1,k).


## 4.3.2. Nonlinear MU

Calculate then×1 matrix filter
gain Kk+1:(4)Kk+1=Pk+1∣kHk+1T[Hk+1Pk+1∣kHk+1T+Rk+1]−1.
The filter measurement update state estimate n×1 matrix X^k+1|k+1, due to the observation yk+1, is calculated with(5)X^k+1∣k+1=X^k+1∣k+Kk+1[yk+1−y(X^k+1∣k)],where Rk+1 is the scalar
variance on the observation residual yk+1−y(X^k+1∣k), and y(X^k+1∣k) is a nonlinear
function of X^k+1∣k. Define the error ΔX^k+1∣k+1 in X^k+1∣k+1:(6)ΔX^k+1∣k+1=Xk+1−X^k+1∣k+1.Define the n×n state estimate
error covariance matrix Pk+1∣k+1 with (7)Pk+1∣k+1=E{(ΔX^k+1∣k+1)(ΔX^k+1∣k+1)T}.Bucy and Joseph [3, page 141] recommend
that Pk+1|k+1 should be
calculated with(8)Pk+1∣k+1=Pk+1∣k−T,where(9)T=Pk+1∣kHk+1TR˜k+1−1Hk+1Pk+1∣k,R˜k+1=Hk+1Pk+1∣kHk+1T+Rk+1.Equations (8) and (9) reduce to the form given by
Kalman:(10)Pk+1∣k+1=[I−Kk+1Hk+1]Pk+1∣k.Calculation of Pk+1∣k+1 by (8) and (9) is numerically stable, whereas the Kalman form
calculation is not.


## 4.4. Nonlinear TU and Nonlinear MU

Refer to Section4.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φx. Given an initial time t0, a final time tf, and a force model u(X^(t),t), then φx propagates the
state estimate X^(t0) from t0 to tf using forces u(X^(τ),τ) to get X^(tf). That is,(11)X^(tf)=φx{tf;X^(t0),t0,u(X^(τ),τ),t0≤τ≤tf}.This can be shortened to write(12)X^(tf)=φx{tf;X^(t0),t0},where the use of forces u(X^(τ),τ) is tacitly
implied. Thus, φx is a column
matrix with n elements:(13)φz=[φx1φx2φx3⋮φxn].


## 4.4.1. Nonlinear TU

The nonlinear TU always spans a nonempty time interval
and requires the use of a numerical state estimate integratorφx. Given an initial time t0, a final time tf, and a force model u(X^(t),t), then φx propagates the
state estimate X^(t0) from t0 to tf using forces u(X^(τ),τ) to get X^(tf). That is,(11)X^(tf)=φx{tf;X^(t0),t0,u(X^(τ),τ),t0≤τ≤tf}.This can be shortened to write(12)X^(tf)=φx{tf;X^(t0),t0},where the use of forces u(X^(τ),τ) is tacitly
implied. Thus, φx is a column
matrix with n elements:(13)φz=[φx1φx2φx3⋮φxn].


## 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 definedobservability 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 isobservable to that observation. Otherwise,
that state estimate component is not observable (unobservable) to
that observation.Theorem 1.
If every component of the row matrixHk+1 of
measurement-state partial derivatives is zero at time tk+1, then every component of the state estimate X^k+1 is unobservable at time tk+1.Proof.
Hk+1=0 implies that Pk+1∣k+1=Pk+1∣k according to (10). Thus none of the variances of Pk+1∣k are reduced due
to processing the observation yk+1. Then by definition, X^k+1 is unobservable in every component.Theorem 2.
Given values for scalarsHk+1, Pk+1∣k>0, Rk+1>0 at time tk+1, and given that Hk+1≠0, then the scalar state estimate X^k+1 is observable at time tk+1.Proof.
The obvious inequalityPk+1∣kHk+12+Rk+1>Pk+1∣kHk+12>0 implies that(14)1>Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1>0.Multiply through by −1:(15)−1<−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1<0. Add 1:(16)0<[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]<1. Multiply through by Pk+1∣k:(17)0<[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]Pk+1∣k<Pk+1∣k.
Now use (4) and (10) to write
(18)Pk+1∣k+1=[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]Pk+1∣k.
Insert (18) into the inequality (17) to get the result:
(19)0<Pk+1|k+1<Pk+1|k.
Thus the variance Pk+1∣k is reduced due
to processing the observation yk+1. Then the scalar state X^k+1 is observable
by definition.


### 5.2. Theoretical Foundation

These theorems are referred to expressions given by
Kalman for filter gainKk+1 and covariance Pk+1∣k+1, 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 matricesPk+1∣k and Pk+1∣k+1 at each time tk+1, 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.


## 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 isobservable to that observation. Otherwise,
that state estimate component is not observable (unobservable) to
that observation.Theorem 1.
If every component of the row matrixHk+1 of
measurement-state partial derivatives is zero at time tk+1, then every component of the state estimate X^k+1 is unobservable at time tk+1.Proof.
Hk+1=0 implies that Pk+1∣k+1=Pk+1∣k according to (10). Thus none of the variances of Pk+1∣k are reduced due
to processing the observation yk+1. Then by definition, X^k+1 is unobservable in every component.Theorem 2.
Given values for scalarsHk+1, Pk+1∣k>0, Rk+1>0 at time tk+1, and given that Hk+1≠0, then the scalar state estimate X^k+1 is observable at time tk+1.Proof.
The obvious inequalityPk+1∣kHk+12+Rk+1>Pk+1∣kHk+12>0 implies that(14)1>Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1>0.Multiply through by −1:(15)−1<−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1<0. Add 1:(16)0<[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]<1. Multiply through by Pk+1∣k:(17)0<[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]Pk+1∣k<Pk+1∣k.
Now use (4) and (10) to write
(18)Pk+1∣k+1=[1−Pk+1∣kHk+12Pk+1∣kHk+12+Rk+1]Pk+1∣k.
Insert (18) into the inequality (17) to get the result:
(19)0<Pk+1|k+1<Pk+1|k.
Thus the variance Pk+1∣k is reduced due
to processing the observation yk+1. Then the scalar state X^k+1 is observable
by definition.


## 5.2. Theoretical Foundation

These theorems are referred to expressions given by
Kalman for filter gainKk+1 and covariance Pk+1∣k+1, 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 matricesPk+1∣k and Pk+1∣k+1 at each time tk+1, 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

LettTNh denote time of
radio wave transmission for the hth NAVSTAR
clock, and let tRGi denote time of
radio wave receipt for the ith 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 δx^TNh and δx^RGi denote Kalman
filter estimation errors in clock phase for tTNh and tRGi. Define time of transmission difference tTDh and time of
receipt difference tRDi:(20)tTDh=tTNh−δx^TNh,tRDi=tRGi−δx^RGi.Thus,(21)tTNh=tTDh+δx^TNh,tRGi=tRDi+δx^RGi.Equation (21) present tTNh and tRGi as additive
combinations of deterministic times tTDh and tRDi and Kalman
filter estimation errors in clock phase δx^TNh and δx^RGi. Define the one-way GPS pseudorange measurement ρNhGi:(22)ρNhGi=c(tRGi−tTNh),tRGi>tTNh. Insert (21) into (22):(23)ρNhGi=c([tRDi+δx^RGi]−[tTDh+δx^TNh])=c([tRDi−tTDh]+[δx^RGi−δx^TNh]),where c is speed of
light in vacuum. Define (24)Δt=tRDi−tTDh,(25)δthi=δx^RGi−δx^TNh.Then,(26)ρNhGi=c(Δt+δthi),where Δt is deterministic
and δt is random.


### 6.2. Partition of Kalman Filter Estimation Errors

LetxC denote the
phase component of Kalman filter estimation error that is common to
every GPS ensemble clock, when it exists. Define phase differences xORGi and xOTNh with (27)xORGi=δx^RGi−xC,xOTNh=δx^TNh−xCfor ground station i and NAVSTAR h. Then Kalman filter estimation errors δx^RGi, i∈{1,2,…}, for ground station clocks and δx^TNh, h∈{1,2,…}, for NAVSTAR clocks have the additive partition:(28)δx^RGi=xC+xORGi,δx^TNh=xC+xOTNh.(This partition
was introduced by Brown [2].)


### 6.3. The Common Random Phase Component  Is Unobservable

Insert (28) into (25):(29)δthi=δx^RGi−δx^TNh=[xC+xORGi]−[xC+xOTNh]=xORGi−xOTNh. Insert (29) into (26):(30)ρNhGi=c(Δt+[xORGi−xOTNh]).Thus, the random phase component xC that is common
to the Kalman filter estimation error for every ensemble clock has vanished in
the range representation ρNhGi. Variations ΔxC in xC cannot cause
variations ΔρNhGi in ρNhGi:(31)ΔρNhGi=∂ρNhGi∂xCΔxCbecause the partial derivative H=∂ρNhGi/∂xC is zero:(32)∂ρNhGi∂xC=0.
An application of Theorem 1 to (32) demonstrates
that xC is unobservable
from ρNhGi.But the architect who designs the complete estimator
must design anoptimal NAVSTAR orbit estimator to prevent aliasing from
NAVSTAR orbit estimation errors into xC. It helps to know that there is no coupling between
the orbit and xC 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 deviationsxORGi and xOTNh are observable
to GPS pseudorange observations because their partial derivatives are nonzero:(33)∂ρNhGi∂xORGi=+c,∂ρNhGi∂xOTNh=−c.Estimation of xORGi and xOTNh 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
processingsimulated GPS pseudorange data.) See (28). xC is the UECC, xORGi is the OEIC for
ground station clocks and xOTNh 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 Figure4).Figure4 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.


## 6.1. GPS Pseudorange Representation

LettTNh denote time of
radio wave transmission for the hth NAVSTAR
clock, and let tRGi denote time of
radio wave receipt for the ith 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 δx^TNh and δx^RGi denote Kalman
filter estimation errors in clock phase for tTNh and tRGi. Define time of transmission difference tTDh and time of
receipt difference tRDi:(20)tTDh=tTNh−δx^TNh,tRDi=tRGi−δx^RGi.Thus,(21)tTNh=tTDh+δx^TNh,tRGi=tRDi+δx^RGi.Equation (21) present tTNh and tRGi as additive
combinations of deterministic times tTDh and tRDi and Kalman
filter estimation errors in clock phase δx^TNh and δx^RGi. Define the one-way GPS pseudorange measurement ρNhGi:(22)ρNhGi=c(tRGi−tTNh),tRGi>tTNh. Insert (21) into (22):(23)ρNhGi=c([tRDi+δx^RGi]−[tTDh+δx^TNh])=c([tRDi−tTDh]+[δx^RGi−δx^TNh]),where c is speed of
light in vacuum. Define (24)Δt=tRDi−tTDh,(25)δthi=δx^RGi−δx^TNh.Then,(26)ρNhGi=c(Δt+δthi),where Δt is deterministic
and δt is random.


## 6.2. Partition of Kalman Filter Estimation Errors

LetxC denote the
phase component of Kalman filter estimation error that is common to
every GPS ensemble clock, when it exists. Define phase differences xORGi and xOTNh with (27)xORGi=δx^RGi−xC,xOTNh=δx^TNh−xCfor ground station i and NAVSTAR h. Then Kalman filter estimation errors δx^RGi, i∈{1,2,…}, for ground station clocks and δx^TNh, h∈{1,2,…}, for NAVSTAR clocks have the additive partition:(28)δx^RGi=xC+xORGi,δx^TNh=xC+xOTNh.(This partition
was introduced by Brown [2].)


## 6.3. The Common Random Phase Component  Is Unobservable

Insert (28) into (25):(29)δthi=δx^RGi−δx^TNh=[xC+xORGi]−[xC+xOTNh]=xORGi−xOTNh. Insert (29) into (26):(30)ρNhGi=c(Δt+[xORGi−xOTNh]).Thus, the random phase component xC that is common
to the Kalman filter estimation error for every ensemble clock has vanished in
the range representation ρNhGi. Variations ΔxC in xC cannot cause
variations ΔρNhGi in ρNhGi:(31)ΔρNhGi=∂ρNhGi∂xCΔxCbecause the partial derivative H=∂ρNhGi/∂xC is zero:(32)∂ρNhGi∂xC=0.
An application of Theorem 1 to (32) demonstrates
that xC is unobservable
from ρNhGi.But the architect who designs the complete estimator
must design anoptimal NAVSTAR orbit estimator to prevent aliasing from
NAVSTAR orbit estimation errors into xC. It helps to know that there is no coupling between
the orbit and xC 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 deviationsxORGi and xOTNh are observable
to GPS pseudorange observations because their partial derivatives are nonzero:(33)∂ρNhGi∂xORGi=+c,∂ρNhGi∂xOTNh=−c.Estimation of xORGi and xOTNh 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
processingsimulated GPS pseudorange data.) See (28). xC is the UECC, xORGi is the OEIC for
ground station clocks and xOTNh 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 Figure4).Figure4 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, σy2(τ) as Allan
variance, a0 as Allan's FMWN
coefficient, a−2 as Allan's FMRW
coefficient, σ1 as the FMWN
diffusion coefficient, and σ2 as the FMRW
diffusion coefficient. Then,(34)σy2(τ)=a0τ−1+a−2τ=σ12τ−1+13σ22τ,where(35)σ1=a0,σ2=3a−2.


### 7.2. Proportionate Process Noise (PPN)

Letα denote a
variable α∈{1,2,3,…,N} to identify
each GPS clock in an ensemble of N clocks. For
each clock α define the
ratio Sα between
diffusion coefficients σ1α and σ2α:(36)Sα=σ2ασ1α,σ1α>0.Then PPN is defined when, for each GPS clock α and each
associated ratio Sα, we have(37)S1=S2=S3=⋯=SN.


### 7.3. Case 12

The calculation ofSα=σ2α/σ1α, α∈{1,2,3,4}, according to the diffusion coefficient values
presented in Figure 1 shows that PPN is not satisfied for Case 12:(38)σ2S1σ1S1=0.00×10−6s−1,σ2S2σ1S2=4.73×10−6s−1,σ2N1σ1N1=2.50×10−6s−1,σ2N2σ1N2=2.67×10−6s−1.


## 7.1. Allan Coefficients versus Diffusion Coefficients

Denoteτ as clock
averaging time, σy2(τ) as Allan
variance, a0 as Allan's FMWN
coefficient, a−2 as Allan's FMRW
coefficient, σ1 as the FMWN
diffusion coefficient, and σ2 as the FMRW
diffusion coefficient. Then,(34)σy2(τ)=a0τ−1+a−2τ=σ12τ−1+13σ22τ,where(35)σ1=a0,σ2=3a−2.


## 7.2. Proportionate Process Noise (PPN)

Letα denote a
variable α∈{1,2,3,…,N} to identify
each GPS clock in an ensemble of N clocks. For
each clock α define the
ratio Sα between
diffusion coefficients σ1α and σ2α:(36)Sα=σ2ασ1α,σ1α>0.Then PPN is defined when, for each GPS clock α and each
associated ratio Sα, we have(37)S1=S2=S3=⋯=SN.


## 7.3. Case 12

The calculation ofSα=σ2α/σ1α, α∈{1,2,3,4}, according to the diffusion coefficient values
presented in Figure 1 shows that PPN is not satisfied for Case 12:(38)σ2S1σ1S1=0.00×10−6s−1,σ2S2σ1S2=4.73×10−6s−1,σ2N1σ1N1=2.50×10−6s−1,σ2N2σ1N2=2.67×10−6s−1.


## 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 Figure4. I set the scalar root-variance R for both
measurement simulations and Kalman filter KF1 to R=1 cm. Typically R~1 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 R=1 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 multidimensionalobservable 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.Figure 2
Simulated and estimated phase deviations for four 2-state clocks.Figure 3
KF1 phase errors and KF2 UECC estimate.Figure 4
S1 observable component phase error.


#### 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 Figure2 for an example.


## 8.1. Create GPS Clock Ensemble

Typically, one processes measurements with a Kalman
filter to derive sequential estimates of a multidimensionalobservable 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.Figure 2
Simulated and estimated phase deviations for four 2-state clocks.Figure 3
KF1 phase errors and KF2 UECC estimate.Figure 4
S1 observable component phase error.


### 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.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 Figure2 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 rigoroussimulator. 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 Sections5 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.



---

*Source: 261384-2008-01-03.xml*