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

Computation of spherical harmonic coefficients

and their error estimates

using Least-Squares Collocation.



Department of Geophysics

Juliane Maries Vej 30,

Dk-2100 Copenhagen Ø, Denmark

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, Yij are the surface spherical harmonics, Pi 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 Cij /R and the anomalous potential is obtained by applying the functional Lij on the covariance function, where

Then we have

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


The covariance with an arbitrary functional Lk is then


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

3. Prediction of spherical harmonic coefficients - theory.

An approximation to T determined using LSC will have the form

where bn are the solutions to the normal equations, Ln 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, {Cnm}.

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(Ln,Lij) 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 100 in latitude, and (2) 1600 additional points distributed with 50 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 106.

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, Tzz, 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 Tzz. 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 -45o and -45o in latitude and -90o , 90o in longitude. The LSC estimated error which cf. Table 1 was 0.29*106 became as shown in Table 2.

Order Error Estimate * 106
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.


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.