Methodological framework to transform relative gravity observations to gravity values

Reilly [1] gives the relative gravimeter reading \(x_{i, j}\), made at location \(j\) and time \(t_i\) as

(1)\[ x_{i,j} = g_j − a − b(t_i − t_1) + z_{i, j} − e_i\]

where,

• \(g_j\) is the true gravity value at location \(j\),

• \(−b(t_i − t_1)\) is the value corresponding to the amount meter has drifted by time \(t_i\),

• \(z_{i, j}\) is the tidal gravitational effect at location \(j\) and time \(t_i\) (this is a known quantity that can be computed by the Longman [3] formula),

• \(a\) is the gravimeter zero point (i.e. the absolute gravity value where the relative gravimeters dial reads zero, in the absence of any drift or tidal effects), and

• \(e_i\) is a (small) error term associated to the meter reading at time \(t_i\).

Collecting the known quantities and the error term on the right hand side gives an observation equation of the following form,

(2)\[g_j − a − b(t_i − t_1) = x_{i,j} − z_{i,j} − e_i.\]

Implicit in all surveys using relative gravimeters, there is a subset of the \(N\) observations at \(M\) locations where there are previously recorded gravity values \(\hat g_j\) . With respect to the true gravity value \(g_j\), the values \(\hat g_j\) must also include a (small) error term \(\epsilon_j\). This leads to a second equation of the form,

(3)\[g_j = \hat g_j + \epsilon_j.\]

In order to relate the relative gravity observations to the absolute gravity values, there must be at least one equation of this type. i.e. the relative gravity observations must be made at least one location where gravity is already known. The system of linear equations given by Eqs. (2) and (3) can then be used to determine the best fitting parameters \(\tilde g_j\), \(\tilde a\) and \(\tilde b\) to the observations by least squares.

Method 1: Normal least squares

The system of linear equations given by Eqs. (2) and (3) can be written in the matrix form given by Eq. (4).

(4)\[\begin{split}\begin{vmatrix} P & \boldsymbol{-1} & \boldsymbol{-T}\\ Q & 0 & 0 \end{vmatrix} \quad \begin{vmatrix} \boldsymbol{g}\\ a\\ b \end{vmatrix} \quad = \begin{vmatrix} \boldsymbol{x - z}\\ \boldsymbol{h} \end{vmatrix} \quad + \begin{vmatrix} \boldsymbol{e} \\ \boldsymbol{\epsilon} \end{vmatrix}\end{split}\]

where,

• g is the vector of \(g_j\)’s

• \(P\) is an \(N × M\) matrix with 1’s in the column \(j\) and row \(i\) for observations at locations \(j\) and time \(t_i\), with 0’s everywhere else

• 1 is an \(N× 1\) vector of ones

• \(T\) is a \(N × 1\) vector of times \(t_i − t1\),

• \(Q\) is an \(M × M\) matrix with 1’s on the diagonal for row \(j\) for locations which correspond to Eq. (3) with 0’s everywhere else

• x−z is the vector of meter reading minus the tidal effect for each observation and

• h is a \(M × 1\) vector with entries \(\hat g_i\) for which the rows of matrix \(Q\) have a 1 in them and zero otherwise.

A least squares solution can then be sought which minimizes the sum of the squared error terms in the vectors e and \(\epsilon\). This will give the best fitting parameter estimates of the absolute gravity at the surveyed locations \(\tilde g\), the gravimeter base line \(\tilde a\) and the drift factor \(\tilde b\).

This method does not decouple the reference sites from the system of observation equations and assumes the errors terms \(\epsilon\) are non zero. This allows the least squares solution for the absolute gravity values at the reference sites to differ from the previously recorded values.

Tip

Method 1 is suitable in most cases, and should be applied when there is uncertainty about the reference value.

Method 2: Decoupled least squares

Following Reilly [1]. Reilly’s [1] original method solves Eqs. (2) and (3) by decoupling the reference gravity values from the observation equations.

Decoupling the reference gravity values is performed by modifying Eq. (2). For the observations made at reference gravity sites, the reference gravity values \(\hat g_j\) are subtracted from the right hand side of the observation equations and the value \(g_j\) for that particular location is discarded. This gives a new observation equation of the form,

(5)\[− a − b(t_i − t_1) = x_{i,j} − z_{i,j} − \hat g_j − e_i\]

In terms of Eq. (4), this corresponds to making amendments to \(P\) and x − z. Entries in column \(j\) of \(P\) such that the observation equations are of the form of (5) are set equal to zero, and the values \(\hat g_j\) are subtracted from the vector x − z.

This method fixes the reference gravity values so that they do not change in the least squares solution, forcing the residual error terms \(\epsilon\) to be exactly equal to zero. However the estimated gravity values parameter variances for the reference sites are non-zero, despite the solution fitting exactly to the given previously recorded gravity value.

Tip

Method 2 fixes the reference value by removing it from the solution. It is particularly useful when the primary interest is change in gravity over time. In this case the reference value can be set to an arbitrary value, and changes are determined relative to this value.

Method 3: Constrained least squares

This method solves a set of equations of the form of Eq. (1) with constraint equations that are similar to an error-less form of Eq. (2) using the method of constrained least squares. The least squares matrix equation is formulated as follows.

Consider the equation

(6)\[X\beta = Y − e\]

where, X = [P, −1, −T],

\(\beta\) = \(\begin{vmatrix} \boldsymbol{g} \\ a \\ b\end{vmatrix}\),

Y = [x−z]

with \(P\), \(T\), \(x\), \(z\), \(\begin{vmatrix} \boldsymbol{g} \\ a \\ b\end{vmatrix}\) and \(e\) are as specified in method 1. Constraints are given by the matrix equation,

(7)\[K\beta = c\]

where \(K\) has dimension \(k×(M+2)\) for some \(k\). Here \(k\) is the number of reference gravity sites and so \(K\) has the same form as \(Q\) but with the zero rows removed. \(c\) is a \(k×1\) vector, and contains the specified gravity values in each row for the location corresponding to the column with a 1 in it in matrix \(K\). This equation is similar to an error-less form of Eq. (3). Similar to method 2, the reference gravity values at locations j are fixed. The aim here is to minimize,

(8)\[S = (Y − X\beta)^{\prime}(Y − X\beta) + \lambda ^{\prime}(K\beta − c)^2)\]

where \(\lambda\) is a \(k\) × 1 vector of Lagrange multipliers.

\(S\) has a minimum where (1) \(\frac{\delta S}{\delta \beta} = 0\) and (2) \(\frac{\delta S}{\delta \lambda}= 0\).

This is given by,

(9)\[2X^{\prime}Y − X\lambda) − K^{\prime}\lambda = 0\]

or (for \(\hat \lambda = \lambda /2\))

(10)\[X^{\prime}(Y − X\beta) − K^{\prime} \hat \lambda = 0\]

which implies,

(11)\[X^{\prime}X\beta − K^{\prime} \hat \lambda = X^{\prime}Y\]

and

(12)\[K\beta = c\]

This can be written in matrix form as follows,

(13)\[\begin{split}\begin{vmatrix} X^{\prime}X & K^{\prime} \\ K & 0 \end{vmatrix} \quad \begin{vmatrix} \beta \\ \hat \lambda \end{vmatrix} \quad = \quad \begin{vmatrix} X^{\prime}Y \\ c \end{vmatrix}\end{split}\]

where the best fitting parameters are found by,

(14)\[\begin{split}\begin{vmatrix} \beta \\ \hat \lambda \end{vmatrix} \quad = \begin{vmatrix} X^{\prime}X & K^{\prime} \\ K & 0 \end{vmatrix}^{-1} \quad \begin{vmatrix} X^{\prime}Y \\ c \end{vmatrix}\end{split}\]

The error terms \(\epsilon\) are considered to be zero and the \(\hat g_j\) parameter variances are also zero.

Tip

Method 3 is most appropriate when the reference gravity values have very low uncertainty, since the solution is adjusted to exactly match the reference values. Any errors in the reference values may distort the solution at the other sites, particularly when there is more than one reference station.

Summary

Three methods to transform relative gravity observations to gravity values have been given. Despite their theoretical differences, in general, each method should return almost identical gravity values. However, the derived gravity value variances may differ between the methods due to the theoretical considerations of the precision of the reference gravity values.