Error analysis for the Global Positioning System

The concept of geometric dilution of precision was introduced in the section, error sources and analysis. Computations were provided to show how PDOP was used and how it affected the receiver position error standard deviation.

When visible GPS satellites are close together in the sky (i.e., small angular separation), the DOP values are high; when far apart, the DOP values are low. Conceptually, satellites that are close together cannot provide as much information as satellites that are widely separated. Low DOP values represent a better GPS positional accuracy due to the wider angular separation between the satellites used to calculate GPS receiver position. HDOP, VDOP, PDOP and TDOP are respectively Horizontal, Vertical, Position (3-D) and Time Dilution of Precision.

Figure 3.1 Dilution of Precision of Navstar GPS data from the U.S. Coast Guard provide a graphical indication of how geometry affect accuracy.[8]

We now take on the task of how to compute the dilution of precision terms. As a first step in computing DOP, consider the unit vector from the receiver to satellite i with components ${\displaystyle {\frac {(x_{i}-x)}{R_{i}}}}$, ${\displaystyle {\frac {(y_{i}-y)}{R_{i}}}}$, and ${\displaystyle {\frac {(z_{i}-z)}{R_{i}}}}$ where the distance from receiver to the satellite, ${\displaystyle \ R_{i}}$, is given by:

${\displaystyle R_{i}\,=\,{\sqrt {(x_{i}-x)^{2}+(y_{i}-y)^{2}+(z_{i}-z)^{2}}}}$

where ${\displaystyle \ x,y,}$ and ${\displaystyle \ z}$ denote the position of the receiver and ${\displaystyle \ x_{i},y_{i},}$ and ${\displaystyle \ z_{i}}$ denote the position of satellite i. These x, y, and z components may be components in a North, East, Down coordinate system a South, East, Up coordinate system or other convenient system. Formulate the matrix A as:

${\displaystyle A={\begin{bmatrix}{\frac {(x_{1}-x)}{R_{1}}}&{\frac {(y_{1}-y)}{R_{1}}}&{\frac {(z_{1}-z)}{R_{1}}}&1\\{\frac {(x_{2}-x)}{R_{2}}}&{\frac {(y_{2}-y)}{R_{2}}}&{\frac {(z_{2}-z)}{R_{2}}}&1\\{\frac {(x_{3}-x)}{R_{3}}}&{\frac {(y_{3}-y)}{R_{3}}}&{\frac {(z_{3}-z)}{R_{3}}}&1\\{\frac {(x_{4}-x)}{R_{4}}}&{\frac {(y_{4}-y)}{R_{4}}}&{\frac {(z_{4}-z)}{R_{4}}}&1\end{bmatrix}}}$

The first three elements of each row of A are the components of a unit vector from the receiver to the indicated satellite. The elements in the fourth column are c where c denotes the speed of light. Formulate the matrix, Q, as

${\displaystyle Q=\left(A^{T}A\right)^{-1}}$

This computation is in accordance with Chapter 11 of The global positioning system by Parkinson and Spilker where the weighting matrix, P, has been set to the identity matrix. The elements of the Q matrix are designated as:[9]

${\displaystyle Q={\begin{bmatrix}d_{x}^{2}&d_{xy}^{2}&d_{xz}^{2}&d_{xt}^{2}\\d_{xy}^{2}&d_{y}^{2}&d_{yz}^{2}&d_{yt}^{2}\\d_{xz}^{2}&d_{yz}^{2}&d_{z}^{2}&d_{zt}^{2}\\d_{xt}^{2}&d_{yt}^{2}&d_{zt}^{2}&d_{t}^{2}\end{bmatrix}}}$

The Greek letter ${\displaystyle \ \sigma }$ is used quite often where we have used d. However the elements of the Q matrix do not represent variances and covariances as they are defined in probability and statistics. Instead they are strictly geometric terms. Therefore, d as in dilution of precision is used. PDOP, TDOP and GDOP are given by

${\displaystyle {\begin{aligned}PDOP&={\sqrt {d_{x}^{2}+d_{y}^{2}+d_{z}^{2}}}\\TDOP&={\sqrt {d_{t}^{2}}}=|d_{t}|\\GDOP&={\sqrt {PDOP^{2}+TDOP^{2}}}\end{aligned}}}$

in agreement with "Section 1.4.9 of PRINCIPLES OF SATELLITE POSITIONING".

The horizontal dilution of precision, ${\displaystyle HDOP={\sqrt {d_{x}^{2}+d_{y}^{2}}}}$, and the vertical dilution of precision, ${\displaystyle \ VDOP={\sqrt {d_{z}^{2}}}=|d_{z}|}$, are both dependent on the coordinate system used. To correspond to the local horizon plane and the local vertical, x, y, and z should denote positions in either a North, East, Down coordinate system or a South, East, Up coordinate system.

The equations for computing the geometric dilution of precision terms have been described in the previous section. This section describes the derivation of these equations. The method used here is similar to that used in "Global Positioning System (preview) by Parkinson and Spiker"

Consider the position error vector, ${\displaystyle \mathbf {e} }$, defined as the vector from the intersection of the four sphere surfaces corresponding to the pseudoranges to the true position of the receiver.${\displaystyle \mathbf {e} =e_{x}{\hat {x}}+e_{y}{\hat {y}}+e_{z}{\hat {z}}}$ where bold denotes a vector and ${\displaystyle {\hat {x}}}$, ${\displaystyle {\hat {y}}}$, and ${\displaystyle {\hat {z}}}$ denote unit vectors along the x, y, and z axes respectively. Let ${\displaystyle \ e_{t}}$ denote the time error, the true time minus the receiver indicated time. Assume that the mean value of the three components of ${\displaystyle \mathbf {e} }$ and ${\displaystyle \ e_{t}}$ are zero.

${\displaystyle A{\begin{bmatrix}e_{x}\\e_{y}\\e_{z}\\e_{t}\end{bmatrix}}={\begin{bmatrix}{\frac {(x_{1}-x)}{R_{1}}}&{\frac {(y_{1}-y)}{R_{1}}}&{\frac {(z_{1}-z)}{R_{1}}}&1\\{\frac {(x_{2}-x)}{R_{2}}}&{\frac {(y_{2}-y)}{R_{2}}}&{\frac {(z_{2}-z)}{R_{2}}}&1\\{\frac {(x_{3}-x)}{R_{3}}}&{\frac {(y_{3}-y)}{R_{3}}}&{\frac {(z_{3}-z)}{R_{3}}}&1\\{\frac {(x_{4}-x)}{R_{4}}}&{\frac {(y_{4}-y)}{R_{4}}}&{\frac {(z_{4}-z)}{R_{4}}}&1\end{bmatrix}}{\begin{bmatrix}e_{x}\\e_{y}\\e_{z}\\e_{t}\end{bmatrix}}={\begin{bmatrix}e_{1}\\e_{2}\\e_{3}\\e_{4}\end{bmatrix}}\ (1)}$

where ${\displaystyle e_{1}}$, ${\displaystyle e_{2}}$, ${\displaystyle e_{3}}$, and ${\displaystyle e_{4}}$ are the errors in pseudoranges 1 through 4 respectively. This equation comes from linearizing the Newton-Raphson equation relating pseudoranges to receiver position, satellite positions, and receiver clock errors. Multiplying both sides by ${\displaystyle A^{-1}}$ there results

${\displaystyle {\begin{bmatrix}e_{x}\\e_{y}\\e_{z}\\e_{t}\end{bmatrix}}=A^{-1}{\begin{bmatrix}e_{1}\\e_{2}\\e_{3}\\e_{4}\end{bmatrix}}\ (2)}$ .

Transposing both sides:

${\displaystyle {\begin{bmatrix}e_{x}&e_{y}&e_{z}&e_{t}\end{bmatrix}}={\begin{bmatrix}e_{1}&e_{2}&e_{3}&e_{4}\end{bmatrix}}\left(A^{-1}\right)^{T}\ (3)}$ .

Post multiplying the matrices on both sides of equation (2) by the corresponding matrices in equation (3), there results

${\displaystyle {\begin{bmatrix}e_{x}\\e_{y}\\e_{z}\\e_{t}\end{bmatrix}}{\begin{bmatrix}e_{x}&e_{y}&e_{z}&e_{t}\end{bmatrix}}=A^{-1}{\begin{bmatrix}e_{1}\\e_{2}\\e_{3}\\e_{4}\end{bmatrix}}{\begin{bmatrix}e_{1}&e_{2}&e_{3}&e_{4}\end{bmatrix}}\left(A^{-1}\right)^{T}\ (4)}$ .

Taking the expected value of both sides and taking the non-random matrices outside the expectation operator, E, there results:

${\displaystyle E\left({\begin{bmatrix}e_{x}\\e_{y}\\e_{z}\\e_{t}\end{bmatrix}}{\begin{bmatrix}e_{x}&e_{y}&e_{z}&e_{t}\end{bmatrix}}\right)=A^{-1}E\left({\begin{bmatrix}e_{1}\\e_{2}\\e_{3}\\e_{4}\end{bmatrix}}{\begin{bmatrix}e_{1}&e_{2}&e_{3}&e_{4}\end{bmatrix}}\right)\left(A^{-1}\right)^{T}\ (5)}$

Assuming the pseudorange errors are uncorrelated and have the same variance, the covariance matrix on the right side can be expressed as a scalar times the identity matrix. Thus

${\displaystyle {\begin{bmatrix}\sigma _{x}^{2}&\sigma _{xy}^{2}&\sigma _{xz}^{2}&\sigma _{xt}^{2}\\\sigma _{xy}^{2}&\sigma _{y}^{2}&\sigma _{yz}^{2}&\sigma _{yt}^{2}\\\sigma _{xz}^{2}&\sigma _{yz}^{2}&\sigma _{z}^{2}&\sigma _{zt}^{2}\\\sigma _{xt}^{2}&\sigma _{yt}^{2}&\sigma _{zt}^{2}&\sigma _{t}^{2}\end{bmatrix}}=\sigma _{R}^{2}\ A^{-1}\left(A^{-1}\right)^{T}=\sigma _{R}^{2}\ \left(A^{T}A\right)^{-1}\ (6)}$

since ${\displaystyle \ A^{-1}\left(A^{-1}\right)^{T}\left(A^{T}A\right)=I}$

Note: ${\displaystyle \left(A^{-1}\right)^{T}=\left(A^{T}\right)^{-1},}$ since ${\displaystyle I=\left(AA^{-1}\right)^{T}=\left(A^{-1}\right)^{T}A^{T}}$

Substituting for ${\displaystyle \left(A^{T}A\right)^{-1}=Q}$ there follows

${\displaystyle {\begin{bmatrix}\sigma _{x}^{2}&\sigma _{xy}^{2}&\sigma _{xz}^{2}&\sigma _{xt}^{2}\\\sigma _{xy}^{2}&\sigma _{y}^{2}&\sigma _{yz}^{2}&\sigma _{yt}^{2}\\\sigma _{xz}^{2}&\sigma _{yz}^{2}&\sigma _{z}^{2}&\sigma _{zt}^{2}\\\sigma _{xt}^{2}&\sigma _{yt}^{2}&\sigma _{zt}^{2}&\sigma _{t}^{2}\end{bmatrix}}=\sigma _{R}^{2}{\begin{bmatrix}d_{x}^{2}&d_{xy}^{2}&d_{xz}^{2}&d_{xt}^{2}\\d_{xy}^{2}&d_{y}^{2}&d_{yz}^{2}&d_{yt}^{2}\\d_{xz}^{2}&d_{yz}^{2}&d_{z}^{2}&d_{zt}^{2}\\d_{xt}^{2}&d_{yt}^{2}&d_{zt}^{2}&d_{t}^{2}\end{bmatrix}}\ (7)}$

From equation (7), it follows that the variances of indicated receiver position and time are

${\displaystyle \sigma _{rc}^{2}=\sigma _{x}^{2}+\sigma _{y}^{2}+\sigma _{z}^{2}=\sigma _{R}^{2}\left(d_{x}^{2}+d_{y}^{2}+d_{z}^{2}\right)=PDOP^{2}\sigma _{R}^{2}}$ and

${\displaystyle \sigma _{t}^{2}=\sigma _{R}^{2}d_{t}^{2}=TDOP^{2}\sigma _{R}^{2}}$

The remaining position and time error variance terms follow in a straightforward manner.

According to the theory of relativity, due to their constant movement and height relative to the Earth-centered, non-rotating approximately inertial reference frame, the clocks on the satellites are affected by their speed. Special relativity predicts that the frequency of the atomic clocks moving at GPS orbital speeds will tick more slowly than stationary ground clocks by a factor of ${\displaystyle {\frac {v^{2}}{2c^{2}}}\approx 10^{-10}}$, or result in a delay of about 7 μs/day, where the orbital velocity is v = 4 km/s, and c = the speed of light. This time dilation effect has been measured and verified using the GPS.

The effect of gravitational frequency shift on the GPS due to general relativity is that a clock closer to a massive object will be slower than a clock farther away. Applied to the GPS, the receivers are much closer to Earth than the satellites, causing the GPS clocks to be faster by a factor of 5×10⁻¹⁰, or about 45.9 μs/day. This gravitational frequency shift is noticeable.

When combining the time dilation and gravitational frequency shift, the discrepancy is about 38 microseconds per day, a difference of 4.465 parts in 10¹⁰.[16] Without correction, errors of roughly 11.4 km/day would accumulate in the position.[17] This initial pseudorange error is corrected in the process of solving the navigation equations. In addition the elliptical, rather than perfectly circular, satellite orbits cause the time dilation and gravitational frequency shift effects to vary with time. This eccentricity effect causes the clock rate difference between a GPS satellite and a receiver to increase or decrease depending on the altitude of the satellite.

To compensate for the discrepancy, the frequency standard on board each satellite is given a rate offset prior to launch, making it run slightly slower than the desired frequency on Earth; specifically, at 10.22999999543 MHz instead of 10.23 MHz.[18] Since the atomic clocks on board the GPS satellites are precisely tuned, it makes the system a practical engineering application of the scientific theory of relativity in a real-world environment.[19] Placing atomic clocks on artificial satellites to test Einstein's general theory was proposed by Friedwardt Winterberg in 1955.[20]

To calculate the amount of daily time dilation experienced by GPS satellites relative to Earth we need to separately determine the amounts due to special relativity (velocity) and general relativity (gravity) and add them together.

The amount due to velocity will be determined using the Lorentz transformation. This will be:

${\displaystyle {\frac {1}{\gamma }}={\sqrt {1-{\frac {v^{2}}{c^{2}}}}}}$

For small values of v/c, by using binomial expansion this approximates to:

${\displaystyle {\frac {1}{\gamma }}\approx 1-{\frac {v^{2}}{2c^{2}}}}$

The GPS satellites move at 3874 m/s relative to Earth's center.[18] We thus determine:

${\displaystyle {\frac {1}{\gamma }}\approx 1-{\frac {3874^{2}}{2\left(2.998\times 10^{8}\right)^{2}}}\approx 1-8.349\times 10^{-11}}$

This difference below 1 of 8.349×10⁻¹¹ represents the fraction by which the satellites' clocks move slower than Earth's. It is then multiplied by the number of nanoseconds in a day:

${\displaystyle -8.349\times 10^{-11}\times 60\times 60\times 24\times 10^{9}\approx -7214{\text{ ns}}}$

That is, the satellites' clocks lose 7,214 nanoseconds a day due to special relativity effects.

Note that this speed of 3874 m/s is measured relative to Earth's center rather than its surface where the GPS receivers (and users) are. This is because Earth's equipotential makes net time dilation equal across its geodesic surface.[21] That is, the combination of Special and General effects make the net time dilation at the equator equal to that of the poles, which in turn are at rest relative to the center. Hence we use the center as a reference point to represent the entire surface.

The amount of dilation due to gravity will be determined using the gravitational time dilation equation:

${\displaystyle {\frac {1}{\gamma }}={\sqrt {1-{\frac {2GM}{rc^{2}}}}}}$

For small values of M/r, by using binomial expansion this approximates to:

${\displaystyle {\frac {1}{\gamma }}\approx 1-{\frac {GM}{rc^{2}}}}$

We are again only interested in the fraction below 1, and in the difference between Earth and the satellites. To determine this difference we take:

${\displaystyle \Delta \left({\frac {1}{\gamma }}\right)\approx {\frac {GM_{\text{earth}}}{R_{\text{earth}}c^{2}}}-{\frac {GM_{\text{earth}}}{R_{\text{gps}}c^{2}}}}$

Earth has a radius of 6,357 km (at the poles) making R_earth = 6,357,000 m and the satellites have an altitude of 20,184 km[18] making their orbit radius R_gps = 26,541,000 m. Substituting these in the above equation, with M_earth = 5.974×10²⁴, G = 6.674×10⁻¹¹, and c = 2.998×10⁸ (all in SI units), gives:

${\displaystyle \Delta \left({\frac {1}{\gamma }}\right)\approx 5.307\times 10^{-10}}$

This represents the fraction by which the satellites' clocks move faster than Earth's. It is then multiplied by the number of nanoseconds in a day:

${\displaystyle 5.307\times 10^{-10}\times 60\times 60\times 24\times 10^{9}\approx 45850{\text{ ns}}}$

That is, the satellites' clocks gain 45,850 nanoseconds a day due to general relativity effects. These effects are added together to give (rounded to 10 ns):

45850 – 7210 = 38640 ns

Hence the satellites' clocks gain approximately 38,640 nanoseconds a day or 38.6 μs per day due to relativity effects in total.

In order to compensate for this gain, a GPS clock's frequency needs to be slowed by the fraction:

5.307×10⁻¹⁰ –  8.349×10⁻¹¹ = 4.472×10⁻¹⁰

This fraction is subtracted from 1 and multiplied by the pre-adjusted clock frequency of 10.23 MHz:

(1 –  4.472×10⁻¹⁰) × 10.23 = 10.22999999543

That is, we need to slow the clocks down from 10.23 MHz to 10.22999999543 MHz in order to negate the effects of relativity.

GPS observation processing must also compensate for the Sagnac effect. The GPS time scale is defined in an inertial system but observations are processed in an Earth-centered, Earth-fixed (co-rotating) system, a system in which simultaneity is not uniquely defined. A coordinate transformation is thus applied to convert from the inertial system to the ECEF system. The resulting signal run time correction has opposite algebraic signs for satellites in the Eastern and Western celestial hemispheres. Ignoring this effect will produce an east–west error on the order of hundreds of nanoseconds, or tens of meters in position.[22]