Ronni Grapenthin - Notes email: |
New Mexico Tech Dept. of Earth & Environmental Science 801 Leroy Place Socorro, NM-87801 |

Published: 2019/03/11

This is the weby-fied version of a handout for my course on geodetic methods and modeling; might as well share it more broadly. How can we get a position estimate from signals a GPS / GNSS satellite sends?

A **range** describes the geometric distance between two points, in our case they are
a satellite and a receiver. This could be inferred by (1) measuring the transit time,
$\tau $,
of a signal that travels from
satellite to receiver at the speed of light,
$c$,
or (2) by calculating the Euclidean distance between known satellite and known receiver positions.

The problem with (2) is that we don't have the receiver's position and (1) depends
on the receiver clock that is biased (non-atomic clock!) and the signal is affected by
path delays due to ionosphere and toposphere impacts and other error sources. That's why
a travel-time derived observable of distance between satellite and receiver is called **pseudorange**.

In the derivation below we will use pseudorange observables in a simple method to estimate receiver positions. We will not deal with all the error terms (e.g., ionosphere or troposphere) and instead assume that they can be easily integrated.

The pseudorange from receiver $u$ to satellite $s$, ${\rho}^{(s)}$, can be expressed as:

$${\rho}^{(s)}={r}^{(s)}+c(\delta {t}_{u}-\delta {t}^{(s)})+I+T+\u03f5$$ |

where ${r}^{(s)}$ is the true range to satellite $s$, $c$ remains the speed of light, $\delta {t}_{u}$ is the receiver clock bias, $\delta {t}^{(s)}$ is clock bias of satellite $s$ and $I,T$ are ionospheric and tropospheric delays. The last term, $\u03f5$, captures unmodeled effects, such as multipath, measurement errors, etc. Note that subscripts (e.g., $u$) reflect receiver specific values, while superscripts identify individual satellites; these are not powers of $(s)$!

As mentioned in the introduction, we have two ways to express the geometric range: the transit time, $\tau $ multiplied by the speed of light, $c$, and the Euclidean distance. The former does not include the receiver position, which we would like to estimate. Thus, we use Euclidean distance here:

$${r}^{(s)}=\sqrt{{({x}^{(s)}-x)}^{2}+{({y}^{(s)}-y)}^{2}+{({z}^{(s)}-z)}^{2}}$$ |

giving the Euclidean distance between a receiver at position $(x,y,z)$ and the satellite, $s$, at position ${p}^{(s)}=({x}^{(s)},{y}^{(s)},{z}^{(s)})$. Note that here, we assume that the satellite position is given. Several levels of quality of these positions exist from broadcast (sent by satellite; about 100 cm precision) to high-precision-two-weeks-after-the-fact-satellite-positions (about 2.5 cm precision), see IGS website for more.

Note that range ${r}^{(s)}$ is also time dependent, which is not explicitly stated above. Obviously, the satellite moves and its position changes. But our receiver may also move (which is actually the interesting application for geophysics). Furthermore, both positions have to be given in the same coordinate system, which we will require to be earth centered earth-fixed, ECEF (details on coordinate systems may come later, for now assume that we can get the satellite position in ECEF).

Welcome to the meat of this post! Our goal is some sort of solution that gives us the receiver position. Let's have at it!

Inserting the Euclidean distance expression for ${r}^{(s)}$ into the pseudorange measurement model we introduced above, we get:

$${\rho}^{(s)}=\sqrt{{({x}^{(s)}-x)}^{2}+{({y}^{(s)}-y)}^{2}+{({z}^{(s)}-z)}^{2}}+c\delta {t}_{u}-c\delta {t}^{(s)}+\u03f5$$ |

where $x,y,z,\delta {t}_{u}$ are unknown and $\u03f5$ now captures all delays including ionosphere and troposphere delays given separately before. Unfortunately, this equation is non-linear in $x,y,z$. This prevents us from directly setting up a linear system of equations that could be easily solved with, for instance, least-squares approximations. What to do? Instead of throwing up our hands and walking away from the problem, we can try to find a linear approximation of the problem and solve that for receiver position and receiver clock error.

To achieve this, your Calc III knowledge will come in quite handy! We'll use the linear parts of a **multivariate Taylor Series expansion** of
${\rho}^{(s)}$, which
assumes that we can approximate ${\rho}^{(s)}$ with a linear function in the vicinity of a point.
For any function
$f(x,y)$
that is at least differentiable once, a linear approximation about the point
$(a,b)$
is given by the sum of the function at this point and its partial derivatives at that point:

$$f(x,y)=f(a,b)+\frac{\partial f}{\partial x}(a,b)(x-a)+\frac{\partial f}{\partial y}(a,b)(y-b)$$ |

So, if we linearize ${\rho}^{(s)}$ about an approximate position and expected receiver clock bias $({x}_{0},{y}_{0},{z}_{0},{t}_{{e}_{0}})$ using the linear part of the multivariate Taylor Series expansion, we get:

$${\rho}^{(s)}(x,y,z,{t}_{e})={\rho}^{(s)}({x}_{0},{y}_{0},{z}_{0},{t}_{{e}_{0}})+\frac{\partial {\rho}^{(s)}}{\partial x}(x-{x}_{0})+\frac{\partial {\rho}^{(s)}}{\partial y}(y-{y}_{0})+\frac{\partial {\rho}^{(s)}}{\partial z}(z-{z}_{0})+\frac{\partial {\rho}^{(s)}}{\partial {t}_{e}}({t}_{e}-{t}_{{e}_{0}})+\u03f5$$ |

Note that we substituted $\delta {t}_{u}$ with ${t}_{e}$ to avoid double deltas here. We can simplify this a bit:

${\rho}^{(s)}(x,y,z,{t}_{e})-{\rho}^{(s)}({x}_{0},{y}_{0},{z}_{0},{t}_{{e}_{0}})$ | $=$ | $\frac{\partial {\rho}^{(s)}}{\partial x}}\mathrm{\Delta}x+{\displaystyle \frac{\partial {\rho}^{(s)}}{\partial y}}\mathrm{\Delta}y+{\displaystyle \frac{\partial {\rho}^{(s)}}{\partial z}}\mathrm{\Delta}z+{\displaystyle \frac{\partial {\rho}^{(s)}}{\partial {t}_{e}}}\mathrm{\Delta}{t}_{e}+\u03f5$ | ||

$\mathrm{\Delta}{\rho}^{(s)}$ | $=$ | $\left[\begin{array}{cccc}\hfill {\displaystyle \frac{\partial {\rho}^{(s)}}{\partial x}}\hfill & \hfill {\displaystyle \frac{\partial {\rho}^{(s)}}{\partial y}}\hfill & \hfill {\displaystyle \frac{\partial {\rho}^{(s)}}{\partial z}}\hfill & \hfill {\displaystyle \frac{\partial {\rho}^{(s)}}{\partial {t}_{e}}}\hfill \end{array}\right]\left[\begin{array}{c}\hfill \mathrm{\Delta}x\hfill \\ \hfill \mathrm{\Delta}y\hfill \\ \hfill \mathrm{\Delta}z\hfill \\ \hfill \mathrm{\Delta}{t}_{e}\hfill \end{array}\right]+\u03f5$ |

Remember that $\mathrm{\Delta}{\rho}^{(s)}$ is the difference between the measured pseudorange and the expected geometric range between a satellite position and the apriori position. We can calculate an updated absolute position and clock bias by adding $[\mathrm{\Delta}x,\mathrm{\Delta}y,\mathrm{\Delta}z,\mathrm{\Delta}{t}_{e}]$ to the apriori values $[{x}_{0},{y}_{0},{z}_{0},{t}_{{e}_{0}}]$. We're pretty close to a solution here, but we first need to calculate the partial derivatives:

$\left[\begin{array}{cccc}\hfill \frac{\partial {\rho}^{(s)}}{\partial x}\hfill & \hfill \frac{\partial {\rho}^{(s)}}{\partial y}\hfill & \hfill \frac{\partial {\rho}^{(s)}}{\partial z}\hfill & \hfill \frac{\partial {\rho}^{(s)}}{\partial {t}_{e}}\hfill \end{array}\right]$ |

Let's work on this for the term $\frac{\partial {\rho}^{(s)}}{\partial x}$. We will need the chain rule (yay, Calc I!):

$$\frac{\partial {u}^{n}}{\partial x}=n{u}^{n-1}\frac{\partial u}{\partial x}$$ |

and we set $u$ to be the term under the square-root in the range expression for the Euclidean distance given above:

$$u={({x}^{(s)}-x)}^{2}+{({y}^{(s)}-y)}^{2}+{({z}^{(s)}-z)}^{2}$$ |

we can write:

$\frac{\partial {\rho}^{(s)}}{\partial x}$ | $=$ | $\frac{\partial \sqrt{{({x}^{(s)}-x)}^{2}+{({y}^{(s)}-y)}^{2}+{({z}^{(s)}-z)}^{2}}}{\partial x}$ | ||

$=$ | $\frac{\partial \sqrt{u}}{\partial x}$ | |||

$=$ | $\frac{1}{2}}{u}^{-\frac{1}{2}}{\displaystyle \frac{\partial u}{\partial x}$ | |||

$=$ | $\frac{1}{2{u}^{\frac{1}{2}}}}{\displaystyle \frac{\partial [{({x}^{(s)}-x)}^{2}+{({y}^{(s)}-y)}^{2}+{({z}^{(s)}-z)}^{2}]}{\partial x}$ | |||

$=$ | $\frac{1}{2{u}^{\frac{1}{2}}}}{\displaystyle \frac{\partial [{({x}^{(s)}-x)}^{2}]}{\partial x}$ | |||

$=$ | $\frac{2({x}^{(s)}-x)}{2{u}^{\frac{1}{2}}}}(-1)$ | |||

$=$ | $\frac{x-{x}^{(s)}}{{\rho}^{(s)}}$ |

Doing this for all the partial derivatives at the apriori position gives us:

$\frac{\partial {\rho}^{(s)}}{\partial x}$ | $=$ | $\frac{{x}_{0}-{x}^{(s)}}{{\rho}_{0}^{(s)}}$ | ||

$\frac{\partial {\rho}^{(s)}}{\partial y}$ | $=$ | $\frac{{y}_{0}-{y}^{(s)}}{{\rho}_{0}^{(s)}}$ | ||

$\frac{\partial {\rho}^{(s)}}{\partial z}$ | $=$ | $\frac{{z}_{0}-{z}^{(s)}}{{\rho}_{0}^{(s)}}$ | ||

$\frac{\partial {\rho}^{(s)}}{\partial {t}_{e}}$ | $=$ | $c$ |

The result for the partial derivative with respect to the clock error follows from earlier expressions of $\delta {t}_{u}$. Note that ${\rho}_{0}^{(s)}$ is the geometric range from the apriori position to satellite $s$, which can be calculated without needing the precise position. Any model corrections that could be applied (e.g., troposphere, ...) could go in there, too. With these expressions for the partial derivatives, we can rewrite the equation for $\mathrm{\Delta}{\rho}^{(s)}$ as:

$\mathrm{\Delta}{\rho}^{(s)}$ | $=$ | $\left[\begin{array}{cccc}\hfill {\displaystyle \frac{{x}_{0}-{x}^{(s)}}{{\rho}_{0}^{(s)}}}\hfill & \hfill {\displaystyle \frac{{y}_{0}-{y}^{(s)}}{{\rho}_{0}^{(s)}}}\hfill & \hfill {\displaystyle \frac{{z}_{0}-{z}^{(s)}}{{\rho}_{0}^{(s)}}}\hfill & \hfill c\hfill \end{array}\right]\left[\begin{array}{c}\hfill \mathrm{\Delta}x\hfill \\ \hfill \mathrm{\Delta}y\hfill \\ \hfill \mathrm{\Delta}z\hfill \\ \hfill \mathrm{\Delta}{t}_{e}\hfill \end{array}\right]+\u03f5$ |

Assuming that we have $n$ satellites in view, each of which giving us a pseudorange measurement ${\rho}^{(1)},\mathrm{\dots},{\rho}^{(n)}$, we can now set up a linear system of equations:

$$\left[\begin{array}{c}\hfill \mathrm{\Delta}{\rho}^{(1)}\hfill \\ \hfill \mathrm{\Delta}{\rho}^{(2)}\hfill \\ \hfill \mathrm{\vdots}\hfill \\ \hfill \mathrm{\Delta}{\rho}^{(n)}\hfill \end{array}\right]=\left[\begin{array}{cccc}\hfill \frac{{x}_{0}-{x}^{(1)}}{{\rho}_{0}^{(1)}}\hfill & \hfill \frac{{y}_{0}-{y}^{(1)}}{{\rho}_{0}^{(1)}}\hfill & \hfill \frac{{z}_{0}-{z}^{(1)}}{{\rho}_{0}^{(1)}}\hfill & \hfill c\hfill \\ \hfill \frac{{x}_{0}-{x}^{(2)}}{{\rho}_{0}^{(2)}}\hfill & \hfill \frac{{y}_{0}-{y}^{(2)}}{{\rho}_{0}^{(2)}}\hfill & \hfill \frac{{z}_{0}-{z}^{(2)}}{{\rho}_{0}^{(2)}}\hfill & \hfill c\hfill \\ \hfill \mathrm{\vdots}\hfill & \hfill \mathrm{\vdots}\hfill & \hfill \mathrm{\vdots}\hfill & \hfill \mathrm{\vdots}\hfill \\ \hfill \frac{{x}_{0}-{x}^{(n)}}{{\rho}_{0}^{(n)}}\hfill & \hfill \frac{{y}_{0}-{y}^{(n)}}{{\rho}_{0}^{(n)}}\hfill & \hfill \frac{{z}_{0}-{z}^{(n)}}{{\rho}_{0}^{(n)}}\hfill & \hfill c\hfill \end{array}\right]\left[\begin{array}{c}\hfill \mathrm{\Delta}x\hfill \\ \hfill \mathrm{\Delta}y\hfill \\ \hfill \mathrm{\Delta}z\hfill \\ \hfill \mathrm{\Delta}{t}_{e}\hfill \end{array}\right]+\u03f5$$ |

This matrix equation is of the form $Gm=d$ where $G$ is a matrix containing the partial derivatives, $d$ is a vector with the pseudorange differences, and $m$ is the vector with the unknowns; specifically the receiver position and clock error. We can solve this with least squares techniques to minimize the sum of squared residuals, for instance, using the normal equations:

$$m={({G}^{T}G)}^{-1}{G}^{T}d$$ |

We can also introduce a weight matrix $W$ to, for instance, put less emphasis on satellites at low elevation angles (I leave it up to you to define reasonable relationships here):

$$m={({G}^{T}WG)}^{-1}{G}^{T}Wd$$ |

Once we have a solution $m=[\mathrm{\Delta}x,\mathrm{\Delta}y,\mathrm{\Delta}z,\mathrm{\Delta}{t}_{e}]$, we can add update the apriori values:

$$\left[\begin{array}{c}\hfill {x}_{new}\hfill \\ \hfill {y}_{new}\hfill \\ \hfill {z}_{new}\hfill \\ \hfill {t}_{{e}_{new}}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {x}_{0}\hfill \\ \hfill {y}_{0}\hfill \\ \hfill {z}_{0}\hfill \\ \hfill {t}_{{e}_{0}}\hfill \end{array}\right]+\left[\begin{array}{c}\hfill \mathrm{\Delta}x\hfill \\ \hfill \mathrm{\Delta}y\hfill \\ \hfill \mathrm{\Delta}z\hfill \\ \hfill \mathrm{\Delta}{t}_{e}\hfill \end{array}\right]$$ |

and iterate until improvements are small.

You can now do what your phone does. If you wish, practice it here in a lab that's part of my class. Use any programming language you like.

The position quality you're getting here isn't fantastic, but good enough to keep you on the road. In science we're usually using dual frequency GPS and track the phase of the signal, which allows for much higher precision. More on that another time, I hope.

rg <at> nmt <dot> edu | Created: 2019/03/11 | Last modified: March 12 2019 15:14.