On the geodetic model and its solution

The full unknown parameter vector

So far in our treatment of geomatics networks we have dealt almost exclusively in two dimensions and with point coordinates, e.g. in \mathbf{F}(\mathbf{x}, \mathbf{l}_{true})=\mathbf{0} we have so far assumed that the observable, \mathbf{x}, has been a coordinate, such as:

    \begin{equation*} \mathbf{x}= \begin{bmatrix} E \\ N \\ \end{bmatrix} \end{equation*}

where E is the easting of the point in question and N is its northing.

But in the bigger picture, i.e. beyond very local networks, each observable needs actually to be represented as a function of a three dimensional point coordinate, \mathbf{r}, and the gravity potential at that point, \mathbf{W}(\mathbf{r}).

In the geodetic sciences we write the desired parameter or observable as follows:

    \begin{equation*} \mathbf{x}= \begin{bmatrix} \text{your position} \\ \text{gravity field at that position} \\ \end{bmatrix}= \begin{bmatrix} \mathbf{r} \\ \mathbf{W}(\mathbf{r}) \\ \end{bmatrix} \end{equation*}

And in that case, \mathbf{F}(\mathbf{x}, \mathbf{l}_{true})=\mathbf{0} should be thought of as follows:

    \begin{equation*} \mathbf{F}(\begin{bmatrix} \mathbf{r} \\ \mathbf{W}(\mathbf{r}) \\ \end{bmatrix}, \mathbf{l}_{true})=\mathbf{0} \end{equation*}

Of course, this is usually non-linear just as it always has been in our discussions. So we need linearize around approximate values which we’ll represent as \mathbf{r}_0 and \mathbf{W}_0, so that:

    \begin{align*} \mathbf{x}&=\mathbf{x}_0+\boldsymbol{\delta} \\ &=\begin{bmatrix} \mathbf{r}_0 \\ \mathbf{W}_0 \end{bmatrix} +\begin{bmatrix} \delta\mathbf{r} \\ \delta\mathbf{W}(\mathbf{r}) \end{bmatrix} \\ &=\begin{bmatrix} \mathbf{r}_0 \\ \mathbf{U}(\mathbf{r}_0) \end{bmatrix} +\begin{bmatrix} \delta\mathbf{r} \\ \mathbf{T}(\mathbf{r}) \end{bmatrix} \end{align*}

i.e. we approximate \mathbf{r} by \mathbf{r}_0 and \mathbf{W} by \mathbf{W}_0 where:

\mathbf{W} is the actual gravity potential

\mathbf{U} is the normal (ellipsoidal) gravity potential

\mathbf{T} is the anomalous gravity potential

And for those not familiar, gravity potential is related to gravity itself as follows:

    \begin{equation*} \mathbf{g}=\nabla(\mathbf{W})=grad(\mathbf{W})=\begin{bmatrix} \dfrac{\partial\mathbf{W}}{\partial x} \\[2ex] \dfrac{\partial\mathbf{W}}{\partial y} \\[2ex] \dfrac{\partial\mathbf{W}}{\partial z} \end{bmatrix} \end{equation*}

where \mathbf{g} is the Newtonian gravitational attraction (with units of m/s^2) and \mathbf{W} is the gravity potential (with units ofm^2/s^2).

So what?!

At this point I’m often asked the question: Why is all of this worth talking about in a networks class?!

Part of the answer to this – as you will see in detail in your geodesy class – is that we need the gravitational potential of the earth in order to define a physically meaningful reference surface for the coordinates that we measure and wish to estimate.

As you will likely already know, the surface of equal potential (equipotential surface) that best approximates mean sea level on earth is known as the geoid and gives us a reference surface of fundamental and very meaningful significance.

One of my students once said in a moment of clarity that: “Where you are depends on gravity!”

And this is exactly the point.

We can (and do) make up references surfaces of mathematical convenience, such as the reference ellipsoids that are commonly used in our navigation systems. And while those do approximate reality in an aggregate sense, they also differ from it in some significant ways. Water might flow uphill if heights are referred to an ellipsoid instead of the geoid, as an example. And a leveling survey will not close if carried out over a long enough distance, as another.

In general, gravity and its potential have significant impacts on our work in geomatics networks. For example:

  • In a one-dimensional (height network), gravitational potential defines the reference surface for the measurements (of height) themselves. This means you couldn’t level for any significant distance without knowing the geoid.
  • In the adjustment of two-dimensional networks, differences between the actual gravitational potential and the normal (approximate) potential – expressed above as the disturbing potential, \mathbf{T} – are needed for reducing measured quantities to the ellipsoid
  • In three-dimensional networks adjustments, the same disturbances are needed to relate geodetic coordinates to their astronomic equivalents (e.g. local geodetic to local astronomic coordinate frames).

And these differences can be significant. Practically speaking, theodolites and leveling instruments are oriented with respect to the true vertical, but that can be different in direction from (“deflect” from) our mathematically convenient models by many tens of arc seconds, e.g. it can reach 100″ in magnitude in the Himalaya region and it’s not unusual for it to be as large as 10-50″ in regions with significant changes in subterranean density, or where mountains are present.

And the undulations in height can vary regionally by as much as 10-50 cm and globally by as much as 100 cm.

We discussed all of this at some length in class and you can find many general descriptions of the disturbing potential online (such as this and this). But I strongly recommend you read the articles1 provided to you on D2L which describe the geoid undulation, N, the deflections of the vertical, \xi and \eta, and their impacts on geomatics networks.

I’m going to ask you about these in the self-assessment questions below.

Being from Canada, I also find pretty fascinating and useful Natural Resources Canada’s page on the modernization of their height reference.

On linearizing the geodetic model

By now you know all about linearizing \mathbf{r} but when linearizing \mathbf{W}(\mathbf{r}) one has to consider that the potentials \mathbf{W} and \mathbf{T} can only be described by an infinite number of parameters – the physics are much more complex than the geometry.

Functionals of \mathbf{T}, such as the geoid undulation, N, and the deflections of the vertical, \xi and \eta, are therefore used to represent the potential field, and for linearization.

And we’re going to stick with convention by representing \mathbf{T}(\mathbf{r}) as \delta\mathbf{s} in the following so that:

    \begin{equation*} \boldsymbol{\delta}= \begin{bmatrix} \delta\mathbf{r} \\ \mathbf{T}(\mathbf{r}) \end{bmatrix}= \begin{bmatrix} \delta\mathbf{r} \\ \delta\mathbf{s} \end{bmatrix} \end{equation*}


\delta\mathbf{r} is the now familiar corrections to the initial coordinate approximations; and

\delta\mathbf{s} are corrections to the normal gravity field.

On solving the geodetic model – getting to the least squares estimation equations

We know that the general functional model \mathbf{F}(\mathbf{x}, \mathbf{l}_{true})=\mathbf{0} has the linear form \mathbf{A}\boldsymbol{\delta} - \mathbf{B}\mathbf{e} + \mathbf{w}=\mathbf{0}.

But if we try to solve this for \boldsymbol{\delta} and \mathbf{e}, we get an infinite number of solutions – which doesn’t help much!

So we put a constraint on the problem. In the most general case, where gravity is considered and where \boldsymbol{\delta}=\begin{bmatrix}\delta\mathbf{r} \delta\mathbf{s}\end{bmatrix}^T, we solve it with the constraint that:

(1)   \begin{equation*} \mathbf{e}^T\mathbf{C}_\mathbf{l}^{-1}\mathbf{e} + \mathbf{s}^T\mathbf{C}_\mathbf{s}^{-1}\mathbf{s} = minimum \end{equation*}

i.e. we try to solve for the corrections to the unknowns \mathbf{x} and \mathbf{W} simultaneously. This is called integrated geodesy and it’s solved by something called least squares collocation.

In our case, we will solve just for \mathbf{x} by a method with which you’re already familiar – the method of least squares. (And then we’ll deal with gravity separately, e.g. through the reductions we will do later in this course).

In this case the solution comes about by constraining that:

(2)   \begin{equation*} \mathbf{e}^T\mathbf{C}_\mathbf{l}^{-1}\mathbf{e} = minimum \end{equation*}

In other words:

  • We derive the functional model \mathbf{F}(\mathbf{x}, \mathbf{l}_{true})=\mathbf{0}
  • Then we linearize it to get it into the form \mathbf{A}\boldsymbol{\delta} - \mathbf{B}\mathbf{e} + \mathbf{w}=\mathbf{0}
  • Then we solve this under the constraint that \mathbf{e}^T\mathbf{C}_\mathbf{l}^{-1}\mathbf{e} = minimum

Luckily we don’t need to solve the latter manually each time we come across a problem. Rather, its general forms have long been solved and can be readily applied as long as we structure things in the right way.

And, also luckily, you already know the solution under the names “estimation using a least squares adjustment” or “least squares estimation”.

There is a set of least squares adjustment equations for each of the combined, parametric, and condition models, as you know.

We won’t be doing the proofs of these three solutions here and we’re only really going to focus on the parametric case, for which you’ll find the equations in the next lesson.

A note about the term “least squares”

A more rigorously correct description of the constraint we applied above might be “minimum quadratic form” or “least quadratic form”, e.g. since we’re actually minimizing the quadratic form \mathbf{e}^T\mathbf{C}_\mathbf{l}\mathbf{e}.

But we never go out into the field and refer to “minimum quadratic form adjustments”! (And I wouldn’t recommend it!)

Instead, we use the term “least squares” universally because of the case where the observations are uncorrelated, e.g. where:

    \begin{equation*} \mathbf{C}_\mathbf{l}^{-1}=\dfrac{1}{\sigma^2}\mathbf{I} \end{equation*}

in which case the “minimum quadratic form” constraint simplifies to:

(3)   \begin{equation*} \dfrac{1}{\sigma^2}\mathbf{e}^T\mathbf{e} = minimum \end{equation*}


(4)   \begin{equation*} \dfrac{1}{\sigma^2}(e_1^2 + e_2^2 + ... + e_n^2)= minimum \end{equation*}

In other words, it becomes a case of minimizing the squares – or “least squares”!

The following schedules are for Alex’s in-class students:

Welcome (back)!

Sign in:

RegisterForgot your password?