Doc: H:\cctwp\gocet03.doc
Testing gridding and filtering of GOCE gradiometer data by LeastSquares Collocation using simulated data.
by
C.C.Tscherning
Department of Geophysics,
University of Copenhagen
DRAFT 20030123
Abstract: Leastsquares collocation (LSC) has been used to test filtering and gridding of GOCE gradiometer data. Initially the function of the GRAVSOFT program GEOCOL was checked by using noisefree data with 5 s sampling provided by IAG along a 1 month realistic orbit. Data in a 2^{o} x 2^{o} area was used. For these data it was verified that 0.1 degree gridding in an interior 1^{o} x 1^{o} area of one quantity from the same quantity could be done with an error below 5 % of the signal variance. If several components of the gravity gradient were used, the error decreased slightly.
Noise with standard deviations of 0.050 and 0.150 E, respectively, and a 1 s correlation distance was generated and added to the data. For the 0.050 E noise the result of the gridding was as for the data without noise, when 1 s sampling was used. For the 0.150 E data, the use of totally 2 months of simulated data and all 3 diagonal components of the gradient matrix were necessary in order to obtain a result below 10 % of the signal standard deviation. Further improvement is possible if data from the planned second measurement phase of GOCE are used.
It will most likely not be possible to use leastsquares collocation (LSC) in its most general form for the modelling of the gravity field of the Earth using GOCE satellite gradiometer data. In its general form LSC requires the solution of systems of equations with as many unknowns as observations.
Regional solutions may however be constructed, and such solutions may form the basis for the creation of datasets of gridded values at satellite altitude. Certain types of gridded data (e.g. the second order radial derivative, d^{2}T/dr^{2}=T_{zz}, where T is the anomalous potential) may be used as the basis for the application of the method of Fast Spherical Collocation, see Sansò and Tscherning, (2002).
In order to test the use of LSC, the program GEOCOL (Tscherning, 1975) has been used with test data available from the IAG Special Commission 7. These test data have been generated for a period of one month with a 5 s sampling, along a realistic GOCE orbit. It includes all the 6 second order derivatives of the gravity potential, V, given in a local reference frame aligned with the velocity vector of the satellite and the radius vector. No roll of the satellite was included. The data were generated using EGM96 (Lemoine et al., 1997) to degree 300. The data set did not include any noise. However, noise has been simulated and test were carried out using the noisy data. The dataset has furthermore been densified timewise (to 1 s) and spasewise corresponding to 2 months of observations in order to test the influence of the sampling.
In the following we will describe the results of the testing. Not only did the testing reveal some embarrassing programming errors, but it did also give better insight into which gravity gradient components to use for the gridding. In section 2 the test using data without noise are described, and in section 3 is described how noise was generated and the results of the gridding using data with added noise.
As a reference potential EGM96 to degree 24 was used. This has the effect that the fields becomes more statistically homogeneous. This is not very important for the gridding, but very important for the interpretation of the error estimates.
As a covariance function we used a function defined using the EGM96 errordegree variances e_{i}^{2} from degree 2 to 24 scaled with a factor of a=0.5, the depth to the Bjerhammarsphere of 1061 m and the variance of the gravity anomalies at the mean Earth surface of 408.13 mgal^{2}. A scale factor on the degreevariances from 25 to infinity was used equal to A= 98.79 mgal^{2}/R^{2}, where R is the mean Earth radius. Hence the analytic covariance expression for the anomalous potential is with R_{B} the Bjerhammar sphere radius
P_{i} are the Legendre polynomials, r and r’ are the size of the radius vectors of P, Q andis the spherical distance between the points. Covariances of the gravity gradient components may be evaluated as described in Tscherning(1976, 1991). Subroutines COVAX, COVBX and COVCX used for these computations are included in the program GEOCOL.
Two test areas were selected. One bounded at latitudes –1 and 1 degrees and longitude –1 and 1 degrees. This area included 29 points. The other area was bounded by latitude 81 and 83 degrees and longitude 45 to 50 degrees and included 294 points.
The individual gravity gradients were used to predict the other gradients as well as the components of the gravity vector given in an inertial, Cartesian frame (also computed from EGM96 to degree 300).
Below in Table 1, the results from the first test area are presented. The points in which the values are predicted ar the same points which are associated with the observations. Hence, the test is a kind of consistency test, and it was used to detect errors in the subroutine COVCX.
The results are given as the differences and standard deviations of "observed" minus calculated and the mean value of the error of prediction calculated also by LSC. Units are mEU for the gravity gradients and mgal for the gravity vector components. In the calculations an uncorrelated data noise of 5 mEU was used, i.e. 25 mEU^{2} was added to the diagonal of the normal equations.
x 
y 
z 
xx 
2xy 
xz 
yy 
yz 
zz 

x 
0.23 1.20 
003 010 
011 021 
001 006 
001 008 
000 003 
005 014 

0.34 
025 
032 
010 
018 
012 
020 

y 
0.15 2.23 
061 009 
001 021 
001 024 
001 005 
000 009 
062 006 

0.34 
032 
020 
032 
010 
015 
030 

3.37 0.34 
000 006 
003 009 
000 045 
028 009 
003 026 
027 014 

6.54 
012 
020 
019 
065 
065 
072 

xx 
0.09 0.58 
0.09 0.55 
1.34 0.05 
032 041 
006 013 
001 048 
038 005 
001 028 
039 005 
1.18 
1.16 
0.76 
071 
020 
028 
063 
060 
063 

2xy 
0.09 0.47 
0.01 0.53 
0.11 0.48 
008 005 
054 020 
006 078 
010 002 
007 016 
018 006 
1.8 
0.95 
1.8 
045 
085 
030 
043 
035 
087 

yz 
0.02 0.62 
0.01 0.30 
0.43 0.13 
017 009 
009 014 
000 010 
001 008 
021 059 
017 011 
1.10 
0.40 
0.40 
032 
030 
013 
025 
089 
045 

yy 
0.13 0.91 
0.03 0.68 
0.49 0.85 
005 006 
002 009 
001 013 
009 013 
017 019 
005 030 
1.15 
0.80 
1.20 
062 
021 
064 
074 
030 
061 

xz 
0.03 0.52 
0.00 0.38 
0.08 0.16 
020 007 
019 026 
007 164 
009 014 
002 006 
012 008 
1.10 
0.50 
0.28 
022 
030 
089 
032 
018 
045 

zz 
0.07 0.32 
0.02 0.59 
0.20 0.28 
009 013 
043 028 
002 044 
009 012 
006 020 
023 049 
0.26 
0.75 
0.40 
030 
062 
025 
030 
030 
127 

xxyy 
0.07 0.79 
0.07 0.83 
1.64 0.11 
017 003 
004 005 
001 062 
017 002 
011 023 
034 005 
1.80 
1.00 
0.90 
018 
010 
032 
043 
032 
085 

xy,xy 
0.06 0.49 
0.03 0.63 
1.23 0.05 
005 001 
000 049 
005 001 
008 017 
010 002 

1.69 
0.95 
0.80 
040 
023 
042 
032 
082 

xx,yy, zz 
0.05 0.27 
0.01 0.35 
0.17 0.25 
000 001 
003 002 
000 041 
000 000 
006 009 
000 001 
0.41 
0.70 
0.38 
002 
008 
021 
002 
024 
002 
Table 1. Results of predicting gravity gradient components from gravity vector components and from gravity gradients. EGM96 to degree 24 subtracted. The first column gives the kind of data used for prediction. The x, y, z refers to derivatives (1. and 2.) of T alongtrack, crosstrack and up. The column heading of each of the following columns shows the type of quantity predicted. If the quantity in the column and the row is the same, it is the mean value and standard deviation of the quantity in the upper cell and the standard deviation of the quantity as computed from the covariance function. Otherwise it is the mean and standard deviation of the predicted quantity in the upper cell and below the mean value of the prediction error computed using GEOCOL. Units are mEU for gravity gradients and mgal for the other quantities. Note that all quantities refer to the instrument (velocity vector) oriented frame.
The table shows that there is a very good consistency between the different quantities. Furthermore, the analytically computed error of prediction is large when the error is large and small when it is small as it should be. Similar results were obtained in the area close to the Pole, where the velocity vector direction varies much. It was furthermore throughout the computations verified that the Lapalceequation was fulfilled both in the East/North/Up (ENU) frame and in the instrument frame. This shows that the program works correctly at least for this type of quantities. However, before this conclusion could be reached a number of errors were detected and corrected.
Subsequently data in a 0.1^{o} grid bounded by 0.5^{o} and 0.5^{o }in latitude and longitude were generated, now in a East, North, Up frame, radially oriented. The quite remarkable results are given in Table 2.
Data 
e(ast) 
n(north) 
u(p) 
ee 
2*en 
eu 
nn 
nu 
uu 
xx 
0.50 0.15 
1.30 0.01 
0.45 0.15 
032 004 
007 003 
017 010 
001 000 
042 002 
040 002 
1.10 
0.74 
1.13 
060 
011 
057 
003 
027 
059 

2*xy 
0.46 0.04 
1.89 0.03 
0.01 0.24 
011 001 
000 000 
014 003 
009 001 
069 005 
020 002 
0.80 
0.83 
1.76 
043 
003 
030 
043 
030 
017 

xz 
0.20 0.05 
0.10 0.08 
0.56 0.01 
010 007 
016 009 
000 003 
021 003 
000 001 
011 004 
0.35 
0.28 
1.10 
030 
028 
013 
021 
002 
042 

yy 
0.08 0.00 
0.11 0.03 
0.07 0.02 
000 000 
002 001 
003 001 
007 002 
003 006 
007 002 
0.78 
1.12 
1.14 
003 
016 
028 
059 
058 
059 

yz 
0.13 0.02 
0.43 0.07 
0.64 0.01 
000 004 
007 004 
000 000 
019 005 
005 002 
020 006 
0.29 
0.42 
1.13 
024 
026 
003 
032 
017 
042 

zz 
0.53 0.08 
1.21 0.08 
0.13 0.15 
009 005 
041 011 
018 003 
009 005 
040 003 
000 000 
0.76 
0.71 
0.36 
027 
053 
025 
027 
025 
002 

xx,yy,zz 
0.28 0.01 
1.14 0.00 
0.12 0.13 
000 000 
002 001 
009 001 
001 000 
038 001 
000 000 
0.69 
0.66 
0.33 
002 
007 
021 
002 
020 
002 
Table 2. The table gives results of predictions of data in a grid with spacing 0.1^{o} totally 121 points. The first row gives the mean and standard deviation of the differences "observed" minus predicted. The second row gives the LSC computed error estimate. In this table, e, n, u refer to a usual local East, North, (radially) Up system. Units: mgal for 1derivatives and EU for 2.derivatives. Note that the instrument xaxis in this area points either North and South, due to the rotation of the Earth in an EarthFixed frame. This is the reason results of predicting nn from xx or ee from yy are excellent.
Again the results are striking. The mean values are sometimes large, indicating that we are using a too limited area of "observations". Again the noisefree observations in the calculations were supposed to have an associated noise standard deviation of 5 mEU. Note that the computed errorestimate is large when the bias is large as it should be. This again confirms the correctness of the program.
3. Gridding with noisy data.
3.1. Noise generation.
An accepted model for the (total) gradiometer noise is given in the Granada report (ESA, 1999, Fig.8.2 ). The contineous noise spectral density function given in units of mEU/sqrt(Hz) may be discretized at 1 s interval from 1 to N=100000 s and modeled as 3 linear functions,
The noise may then be represented as a (periodic) Fourier series as a function of time, t,
,
where the coefficients a_{i} and b_{i} are drawn from a normal distribution with variances .
The autocovariance function for the corresponding stationary process will have the variances as coefficients in a cosine series with argument equal to the time difference,
.
The resulting noise standard deviation becomes 56 mEU.
The error is much larger than the 5 mEU normally used in simulations. However, it must be reduced using more data and filtering. Averaging may be used , and a 4 s averaging will reduce the standard deviation 50 %, obviously. Secondly the noise may be reduced by filtering using the fact that the low order harmonics are known very well. In the following a dataset with 1 s sampling as well as a dataset corresponding to 2 months of observations will be used.
A FORTRAN program "ecf.f" was written in order to generate (pseudo) random noise for the evaluation of the covariance function for time differences as well as for intervals of time differences (mean values). An example of the generated noise spectral constituents a_{i} and b_{i} as well of the noise for the time difference interval from 0  40 s are given in the tables below.
Table 3. Wavenumber, fourier coefficients and spectral density
i a_{i} b_{i }
0 .1892 .1143 .1200 1 .0577 .0120 .1121 2 .1465 .0152 .1041 3 .0078 .0838 .0962 4 .0288 .0678 .0883 5 .1397 .0211 .0803 6 .0768 .0320 .0724 7 .0588 .0078 .0645 8 .0960 .0821 .0565
9 .0085 .0158 .0486
10 .0195 .0224 .0407
11 .0049 .0172 .0327
12 .0044 .0044 .0248
13 .0083 .0319 .0169
14 .0059 .0013 .0089
15 .0007 .0002 .0010
16 .0007 .0009 .0010
17 .0002 .0010 .0010
18 .0005 .0000 .0010
19 .0016 .0016 .0010
20 .0007 .0014 .0010

Tabel 4. Noise for the first 50 seconds (units of EU).
Time noise time noise time noise time noise time noise
0 .0480 1 .0168 2 .0428 3 .0506 4 .0942 5 .0942 6 .0425 7 .0869 8 .0969 9 .0531 10 .0531 11 .0860 12 .0187 13 .0155 14 .0942 15 .0942 16 .0055 17 .0417 18 .0158 19 .0307
20 .0307 21 .0455 22 .0211 23 .0578 24 .0588
25 .0588 26 .0268 27 .1180 28 .0632 29 .0004
30 .0004 31 .0310 32 .0125 33 .0328 34 .0134
35 .0134 36 .0662 37 .0412 38 .0493 39 .0125
40 .0125 41 .0224 42 .0041 43 .1092 44 .0964
45 .0964 46 .1083 47 .0165 48 .0498 49 .0119

The autocovariance function of the stationary time series and the one computed forming product sums from the series generated as described above are somewhat different, see Table 5. Both functions show a very small correlation already after 1 s.
Table 5. Covariance function computed from power spectral density (p) and estimated from the generated pseudorandom noise (s), units EU^{2}.
Time dif p s
ference
0 .003144 .003027 1 .000546 .000448 2 .000142 .000031 3 .000067 .000074 4 .000041 .000167 5 .000029 .000027
6 .000022 .000067
7 .000018 .000074
8 .000016 .000173
9 .000014 .000045
10 .000013 .000124

The one generated from pseudorandom noise goes faster to zero, so that the noise becomes uncorrelated after 2 s time difference.
3.2. Gridding with noisy data.
The generated noise was added to the "observations" used in section 2. These data have a 5 s timedifference, hence they will be regarded as uncorrelated. The variance of 0.003027 EU^{2} was added to the diagonal of the normal equations.
The results are presented in Table 6.
Data 
x 
y 
z 
xx 
2*xy 
xz 
yy 
yz 
zz 
xx 
0.52 0.07 
2.06 0.06 
0.18 0.26 
012 002 
029 005 
011 005 
003 002 
073 006 
009 004 
1.30 
0.98 
1.31 
068 
061 
071 
016 
046 
071 

2*xy 
1.00 0.15 
3.80 0.26 
0.71 0.69 
041 007 
013 011 
044 008 
071 027 
190 015 
030 033 
1.19 
1.20 
2.02 
059 
018 
062 
063 
062 
111 

xz 
0.33 0.04 
0.58 0.10 
0.83 0.18 
004 007 
015 010 
001 004 
028 004 
004 006 
032 011 
1.18 
0.47 
1.59 
066 
052 
064 
035 
017 
081 

yy 
0.09 0.11 
0.004 0.14 
0.04 0.03 
002 001 
060 005 
000 007 
023 004 
002 014 
021 005 
1.01 
1.31 
1.32 
015 
065 
048 
069 
071 
071 

yz 
0.37 0.07 
0.76 0.24 
1.10 0.03 
015 010 
012 013 
001 005 
065 005 
008 014 
050 013 
0.46 
1.21 
1.61 
037 
050 
017 
065 
068 
082 

zz 
0.65 0.06 
2.78 0.06 
0.25 0.43 
004 004 
029 004 
019 005 
005 013 
123 002 
009 015 
0.95 
0.91 
0.62 
038 
069 
043 
038 
040 
018 
Table 6. Differences between "observed" and predicted data in a grid. Noisy data with a noise standard deviation of 55 mE was used. First row gives the type of predicted quantities and the first column the kind of data used. For 1. derivatives the units are mgal and for 2. derivatives it is mEU. By a mistake the noise added to the 2*xy component was not multiplied by 2.
In the table we note that not only have the predictions become worse, but also the estimated prediction error from LSC  as it should be. There are still some large biases due to the small area used for the "observations". The observation and the prediction grid points are the same as used in Section 2.
It is encouraging to see that rather small errors of prediction can be expected using T_{zz } and predicting the same kind of data in a grid. The estimated error is 18 mEU. Hence for 5 s means the error will be around 9 mEU. Then we are not far from the 5 mEU required.
A dataset with 1 s sampling was generated and noise added. This improved the results considerably, so that when 3 diagonal components were used a "perfect" filtering was obtained, see Table 7.
Data 
ee 
en*2 
eu 
nn 
nu 
uu 
xx, yy, zz 1 s sampling 
007 004 
002 003 
029 006 
001 002 
054 003 
006 005 
360 obs 
006 
031 
029 
006 
029 
006 
Noise stdv 
xx 
xy*2 
xz 
yy 
yz 
zz 
055 mE 
003 005 
001 010 
001 061 
004 007 
012 030 
001 010 
009 
033 
033 
009 
033 
010 
Table 7. Prediction (row 2 and 3) and filtering (row 6 and 7) using 1 s data and xx, yy, zz as observations with noise standard deviation of 55 mE. Note the very good filtering at the data. The noise is reduced to 10 %.
Now what we have used here is an "optimistic" PSD. Regguzzoni (2002) has used a more pessimistic model, with a corresponding noise variance of 0.032 EU^{2} or a standard deviation of 0.181 EU. Here a noise standard deviation slightly smaller of 0.116 EU was used. The results using the 1 s  1 month dataset were unsatisfactory, and a dataset corresponding to a 2 month period was generated. The results are given in Table 8.
Data 
ee 
en*2 
eu 
nn 
nu 
uu 
xx, yy, zz 1 s sampling 
004 002 
028 004 
029 004 
007 001 
058 005 
011 002 
2 months 
008 
041 
031 
008 
032 
009 
720 obs. 

Noise stdf. 
xx 
xy*2 
xz 
yy 
yz 
zz 
116 mE 
009 006 
026 014 
001 068 
003 005 
009 032 
011 004 
012 
043 
036 
012 
035 
013 
Table 8. Gridding and filtering of data with 116 mE standard deviation noise added. Totally 720 observations, corresponding to 2 months of data. Rows 2 and 3 show the results of the gridding and rows 6 and 7 the filtering.
We see that having the double number of data enables us to obtain very satisfactory results when the data have a twice larger noise standard deviation as shown in Table 7. However a second measurementphase of 2 months (as planned for GOCE) should give results of both the gridding and the filtering which are very satisfactory, i.e. like working with data with a 5 mE noise standard deviation.
Here we should keep in mind, that the error of prediction varies from area to area due to the variation in signal variance. We have here (not by purpose) used an area with a rather smooth gravity field. Similar calculations should be repeated in mountainous areas. Using a larger area may also aid in reducing the longperiodic noise, and hence reduce the noise variance. Furthermore LSC may (as demonstrated above) be used to combine different components, which will reduce the error, or equivalently improve the filtering of the noise. This naturally requires that the noise is uncorrelated between the gravity gradient components.
The tests carried out is both a test of the GEOCOL program as well as of the possible results to be expected when gridding noisy data. A larger area will help in reducing the biases. The investigation of these items will be the object of another study.
Acknowledgement: Thanks to IAG SC7 for providing the test data.
References:
ESA: Gravity Field and SteadyState Ocean Circulation Mission, ESA SP1233, 1999.
Lemoine, F.G., S.C. Kenyon, J.K. Factor, R.G. Trimmer, N.K. Pavlis, D.S. Chinn, C.M. Cox, S.M. Klosko, S.B. Luthcke, M.H. Torrence, Y.M. Wang, R.G. Williamson, E.C. Pavlis, R.H. Rapp, and T.R. Olson, The Development of the Joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) Geopotential Model EGM96, NASA/TP1998206861, Goddard Space Flight Center, Greenbelt, MD, July, 1998.
Reguzzoni, M.: From the timewise to the spacewise GOCE observables. Como, 2002.
Sanso', F. and C.C.Tscherning: Fast spherical collocation: A General Implementation. IAG Symposia, Vol. 125, pp. 131137, Springer Verlag, 2002.
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. 115123, 1993.