27. januar 2000 File: C:\CCTWP\sphcov3.wpd

**Abstract:** The equations expressing the covariances between spherical harmonic coefficients
and linear functionals applied on the anomalous gravity potential, T, are derived. The
functionals are the evaluation functionals, and those associated with first and second order
derivatives of T. These equations forms the basis for the prediction of spherical harmonic
coefficients using Least Squares Collocation (LSC).

This is implemented in the GRAVSOFT program GEOCOL. Initially tests using EGM96
were performed using global and regional sets of geoid heights, gravity anomalies and second
order vertical gravity gradients at ground level and at altitude. The global tests confirm that
coefficients may be estimated consistently using LSC while the error estimates were much
too large for the lower order coefficients.

The validity of an error estimate calculated using LSC with an isotropic covariance function
is based on a hypothesis that the coefficients of a specific degree all belong to the same
normal distribution. However, the coefficients of lower degree do not fulfill this, and this
seems to be the reason for the too pessimistic error estimates. In order to test this the
coefficients of EGM96 were perturbed, so that the pertubations for a specific degree all
belonged to a normal distribution with the variance equal to the mean error variance of the
coefficients. The pertubations were used to generate residual geoid heights, gravity
anomalies and second order vertical gravity gradients. These data were then used to calculate
estimates of the perturbed coefficients as well as error estimates of the quantities, which now
had a very good agreement with the errors computed from the simulated observed minus
calculated coefficients.

Tests with regionally distributed data showed that long-wavelength information is lost, but
also that it seems to be recovered for specific coefficients depending on the location of the
area where the data are located.

**1. Introduction.**

When using least-squares collocation (LSC) available information about the spherical
harmonic coefficients of the anomalous gravity potential, T, can be used directly as
observations, or used in a remove-restore procedure, see Forsberg & Tscherning (1981).
However, if coefficients are to be predicted, we need the explicit covariances (or values of
functionals applied on a reproducing kernel, see Tscherning (1993)). As shown in section 2
the covariances are simply the observation functional applied on the solid spherical harmonic
function of a specific degree and order, multiplied by a constant depending on the degree.

The equations have been implemented in the form of new versions of the subroutines
COVAX, COVBX and COVCX (Tscherning, 1976) in the program COVFIT (Knudsen,
1987).

The equations to be used for the calculation of the estimates of the coefficients and their
error-estimates are derived in Section 3. They are implemented in a new version of the
GRAVSOFT program (Tscherning et al. (1992)) GEOCOL.

In section 4 are described results of coefficient prediction tests initially using the EGM96
(Lemoine et al. (1996)) coefficients from degree 8 to 180. The coefficients were used as
control data and to generate control data sets of geoid heights, gravity anomalies and second
order radial derivatives. The error estimates of the coefficients were compared to the
differences between "observed" and predicted coefficients, and it was found that the errors
were too large for the low degree coefficients. This seemed to be caused by the non-normal
distribution of the low degree coefficients.

The computational experiment was repeated using now as coefficient from degree 2 to 180
pertubations of EGM96 with a normal distribution with a variance which for a given degree
was equal to the error degree variance of the EGM96 coefficients. Using the original EGM96
coefficients as a reference field this resulted in error estimates which were in a very good
agreement with the errors obtained in the calculations.

All derivations and tests have been made in spherical approximation. In order not to use this,
the geodetic latitude must be changed to the geocentric latitude. Furthermore, the value of
the radius-vector (r) which now is calculated as R+h, where R is the Earth mean radius, must
be computed rigorously from the Cartesian coordinates.

**2. Covariances between spherical harmonic coefficients and point related quantities.**

Let P, Q be two points with coordinates (latitude, longitude, r), ,
respectively, and having the spherical distance . R is the mean radius of the Earth, Y_{ij}
are the surface spherical harmonics, P_{i} the Legendre polynomials and the degree-variances. Then the covariance between the values of the anomalous potential T in P, Q is

The covariance between the coefficient GM C_{ij} /R and the anomalous potential is obtained by
applying the functional L_{ij} on the covariance function, where

Then we have

When applying the functional twice we get its norm, squared,

or

The covariance with an arbitrary functional L_{k} is then

Example:

The correlation can also be calculated between for example the value of T in a point P and
the i,j'th coefficient (ev_{P}(T) = T(P) ):

**3. Prediction of spherical harmonic coefficients** - **theory**.

An approximation to T determined using LSC will have the form

where b_{n} are the solutions to the normal equations, L_{n} are the observation functionals and N
is the number of observations. If the data contain errors, the variance-covariance of the noise
is added to the normal-equation matrix, {C_{nm}}.

The value of a predicted quantity is obtained by applying the associated functional on this expression. For the spherical harmonic coefficients we then have

The calculation of the observation functional applied on the solid spherical harmonic
function is generally done by recursion, starting with the (0,0) term. This means that a better
computational strategy would be to calculate all coefficients simultaneously up to and
inclusive degree i. However, we get a problem when we want to calculate the error estimates

because we then must store all the quantities cov(L_{n},L_{ij}) in the expression eq. (10)
simultaneously, and subsequently evaluate the expressions for all degrees up to and inclusive
degree i. This will obviously only be possible when the degree i and the number of
observations are small.

It is worth at this point to recall that since we use an isotropic covariance function (kernel),
cf. eq. 1, the error-estimate is an estimate of the error in a mean square sense, which only
can be interpreted in a standard manner if the coefficients for a given degree all belong to the
same normal distribution, and having a variance given by equation (4).

Finally a remark on the computational procedures. If the coefficients have a normal
distribution with the same variance for each degree, then we may compute error estimates by
generating data from normal distributed perturbed sets of coefficients. We may then
calculate the difference between "observed" coefficients and predicted values using
throughout the same matrix C in equation (10). The mean square difference may then give a
correct estimate of the error.

**4. Numerical tests of coefficient prediction and of error estimation.**

The numerical test data were generated using EGM96 from degree 8 to 180. Here the results
will be illustrated using two point configurations: (1) 400 points distributed in the middle of
equal-area blocks having an extend of 10^{0}** **in latitude, and (2) 1600 additional points
distributed with 5^{0} spacing, totally 2000 points.

Three data types were generated: geoid heights, gravity anomalies and vertical gravity gradients, either located at altitude zero or at varying altitudes, see Table 1 where selected results are presented. The covariance function used, was the one generated from the EGM96 coefficients (the "true" one).

Data | Type | gravity | gravity | gravity | gravity | Tzz |

Height | 0 m | 50 km | 300 km | 300 km | 300 km | |

Order | Number | 400 | 400 | 400 | 2000 | 2000 |

m | EGM96 | |||||

9 | -0.48 | -0.02 | -0.28 | -0.63 | -0.44 | -0.39 |

8 | 1.88 | 0.76 | 1.29 | 1.84 | 1.82 | 1.82 |

7 | -1.18 | -0.24 | -0.37 | -0.84 | -1.18 | -1.26 |

6 | 0.63 | 0.58 | 0.75 | 0.80 | 0.59 | 0.48 |

5 | -0.17 | 0.27 | 0.08 | -0.13 | -0.14 | -0.04 |

4 | -0.09 | -0.08 | -0.19 | -0.23 | -0.06 | 0.04 |

3 | -1.70 | -0.31 | -0.80 | -1.50 | -1.53 | -1.51 |

2 | 0.22 | -0.63 | -0.31 | 0.24 | 0.21 | 0.14 |

1 | 1.43 | 0.81 | 1.00 | 1.12 | 1.34 | 1.22 |

0 | 0.28 | 0.46 | 0.50 | 0.32 | 0.31 | 0.45 |

EGM96 | Mean . | -0.18 | -0.08 | -0.02 | -0.01 | -0.01 |

-Predict | Stand.dev. | 0.64 | 0.40 | 0.18 | 0.07 | 0.14 |

Estim. Err. | 0.79 | 0.60 | 0.34 | 0.29 | 0.32 |

Table 1. Prediction of EGM96 spherical harmonic coefficients of degree 9 and order 9 to 0.

In column 2 are given the "true" values, and in the next column the estimated values. At the
bottom are given the mean of the differences (ERM96 - Predicted), the standard deviation
and the value calculated using equation (10), divided with GM / R, i.e. unitless. All values
have been multiplied by 10^{6}.

From the table we see - as expected - that the result improves with more gravity data and
with a higher altitude. It is remarkable that for this low degree, the second order vertical
gradient, T_{zz}, gives results nearly as good as when gravity was used. However, from Fig. 1 we
see the expected result, namely that geoid heights and gravity give better results for the lower
degrees than T_{zz.} The figure also shows that the main contribution to the coefficients are in
the degrees below 40, corresponding very well to the data spacing of 5 degrees.

The root mean square differences between predicted and "observed" coefficients per degree
are shown in Fig. 2- 4. 2000 values of each of the three data types at an altitude of 300 km
have been used. The figures also show the root mean collocation error estimates (eq. (10))
and the coefficient root variances per degree. These figures show as could already be seen by
inspecting Table 1 that the error estimates are too large. Too pessimistic. The largest
discrepancy occur for the case where geoid heights have been used, and the smallest where
vertical gravity gradients were used.

Since the prediction results were excellent it was believed that the software used in the test
was correct, so the discrepancy would be caused by something else. Recalling the last remark
in the previous section prompted an investigation of whether the EGM96 coefficients were
having a normal distribution for each degree. Histograms per degree were formed, see Fig. 5.
By inspecting the histograms, it immediately becomes clear that the lower order coefficients
do not follow a normal distribution. This also explains why the discrepancy were largest for
geoid heights (which have the largest correlation with the lower order coefficients) and
smallest for the vertical gravity gradient.

In order then to verify that we would be able to obtain valid error estimates for normally
distributed coefficients the following computational experiment was carried out. A generator
of data with a normal distribution with a given variance and zero mean was used to generate
perturbed EGM96 coefficients. The program HARMEXG was used, see the details in
Tscherning et al. (1999). The variance of the normal distribution was put equal to the mean
variance per degree of the EGM96 errors.

The perturbed coefficients were used to generate residual geoid heights, gravity anomalies
and vertical gravity gradients at 300 km altitude as in the first experiment. EGM96 itself was
used as a reference field, and the perturbed coefficients were predicted. The results are given
in Fig. 6 - 9, and we see that the error estimates obtained from comparing observed and
predicted values now agrees well with the collocation error estimates. A small discrepancy is
found for the lowest degrees, which is due to the fact that it is difficult to form a normally
distributed set of numbers from very few values. Furthermore we see that the second order
radial derivative gives very little contribution to the (perturbed) coefficients of the lowest
order.

**5. Prediction of coefficients from a regional data distribution.**

The prediction of spherical harmonic coefficients will only be good, when it is based on a
global data coverage. We know that the estimate of the coefficients will degrade, if we for
example do not have data at the caps at the south and the north pole. But how much will the
estimate degrade if the data only covers a part of the globe ?

Using the error-estimation capability of LSC we are able to study this. A small example is
given here in Table 2, where we have used a subset of 678 values of the 2000 gravity values
used above at 300 km altitude with a distribution limited by -45^{o} and -45^{o} in latitude and ^{ } -90^{o} , 90^{o} in longitude. The LSC estimated error which cf. Table 1 was 0.29*10^{6 } became as
shown in Table 2.

Order | Error Estimate * 10^{6} |

9 | 0.48 |

8 | 0.54 |

7 | 0.56 |

6 | 0.62 |

5 | 0.72 |

4 | 0.71 |

3 | 0.76 |

2 | 0.75 |

1 | 0.79 |

0 | 0.76 |

The interesting fact is that the prediction of some coefficients only becomes slightly worse,

(m = 9). This confirm the conjecture that regional gravity models do contain much correct
information even at the longer wavelengths.

**6. Conclusion.**

It has been demonstrated that LSC may be used successfully used for the prediction of
spherical harmonic coefficients. Meaningful error estimates can be calculated if one works
with pertubations of a set of coefficients and not with the total quantity. This is exactly like
in a network adjustment, where the object not is to find the coordinates, but improvements to
the coordinates.

This will be important in simulation studies, where one tries to understand the influence of
various data types and data distribution on coefficient estimation. The use of the procedures
for much larger data sets will be possible using sparse matrix techniques, a topic where much
progress has been made recently, see Moraux et al. (1999).

**Acknowledgement:** This is a contribution to the ESA supported project Eötvös to Mgal
(E2M), and to the GEOSONAR project funded under the Earth Observation program by the

Danish Research Councils.

**References.**

Forsberg, R. and C.C.Tscherning: The use of Height Data in Gravity Field Approximation by Collocation. J.Geophys.Res., Vol. 86, No. B9, pp. 7843-7854, 1981.

Knudsen, P.: Estimation and Modelling of the Local Empirical covariance Function using gravity and satellite altimeter data. Bulletin Geodesique, Vol. 61, pp. 145-160, 1987.

Lemoine, F.G., D.Smith, R.Smith, L.Kunz, E.Pavlis, N.Pavlis, S.Klosko, M.Torrence, R.Williamson, C.Cox, K.Rachlin, Y.Wang, S.Kenyon, .Salman, R.Trimmer, R.Rapp and S.Nerem: The development of the NASA GSFC and MA joint geopotential model. Proc. Symp. on Gravity, Geoid and Marine Geodesy, Tokyo 1996.

Moreaux, G., C.C.Tscherning & F.Sanso': Approximation of Harmonic Covariance Functions by non Harmonic Locally Supported functions. Journal of Geodesy, Vol. 73, pp. 555-567, 1999.

Tscherning, C.C.: A FORTRAN IV Program for the Determination of the Anomalous Potential Using Stepwise Least Squares Collocation. Reports of the Department of Geodetic Science No. 212, The Ohio State University, Columbus, Ohio, 1974.

Tscherning, C.C.: Covariance Expressions for Second and Lower Order Derivatives of the Anomalous Potential. Reports of the Department of Geodetic Science No. 225, The Ohio State University, Columbus, Ohio, 1976.

Tscherning, C.C.: Computation of covariances of derivatives of the anomalous gravity potential in a rotated reference frame. Manuscripta Geodaetica, Vol. 18, no. 3, pp. 115-123, 1993.

Tscherning, C.C., R.Forsberg and P.Knudsen: The GRAVSOFT package for geoid
determination. Proc. 1. Continental Workshop on the Geoid in Europe, Prague, May 1992,
pp. 327-334, Prague, 1992.

Tscherning, C.C., O. Andersen, D.Arabelos E.Carminati, R.Forsberg, A. Gardi, P.Knudsen, J.N.Larsen, R.Sabadini and G.Strykowski: Refinement of the Current Observation Requirements for GOCE. Final Report ESA - ESTEC Contract 13229/98/NL/GD, Nov. 1999.

Figures: