Doc: H:\cctwp\gocet03.doc

Testing gridding and filtering of GOCE gradiometer data by Least-Squares Collocation using simulated data.

 

by

C.C.Tscherning

Department of Geophysics,

University of Copenhagen

 

DRAFT 2003-01-23

Abstract: Least-squares 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 noise-free data with 5 s sampling provided by IAG along a 1 month realistic orbit. Data in a 2o x 2o area was used. For these data it was verified that 0.1 degree gridding in an interior 1o x 1o 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.

  1. Introduction.
  2. It will most likely not be possible to use least-squares 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, d2T/dr2=Tzz, 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 time-wise (to 1 s) and spase-wise 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.

  3. Testing the LSC prediction procedure with data without 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 error-degree variances ei2 from degree 2 to 24 scaled with a factor of a=0.5, the depth to the Bjerhammar-sphere of 1061 m and the variance of the gravity anomalies at the mean Earth surface of 408.13 mgal2. A scale factor on the degree-variances from 25 to infinity was used equal to A= 98.79 mgal2/R2, where R is the mean Earth radius. Hence the analytic covariance expression for the anomalous potential is with RB the Bjerhammar sphere radius

Pi 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 mEU2 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

xx-yy

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,x-y

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 along-track, cross-track 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 Lapalce-equation 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.1o grid bounded by -0.5o and 0.5o 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.1o 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 1-derivatives and EU for 2.-derivatives. Note that the instrument x-axis in this area points either North and South, due to the rotation of the Earth in an Earth-Fixed 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 noise-free observations in the calculations were supposed to have an associated noise standard deviation of 5 mEU. Note that the computed error-estimate 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 ai and bi are drawn from a normal distribution with variances .

The auto-covariance 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 ai and bi 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 ai bi

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 auto-covariance 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 pseudo-random noise (s), units EU2.

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 pseudo-random 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 time-difference, hence they will be regarded as uncorrelated. The variance of 0.003027 EU2 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 Tzz 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 EU2 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 measurement-phase 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.

 

  1. Conclusion.

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 long-periodic 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 Steady-State Ocean Circulation Mission, ESA SP-1233, 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/TP-1998-206861, Goddard Space Flight Center, Greenbelt, MD, July, 1998.

Reguzzoni, M.: From the time-wise to the space-wise GOCE observables. Como, 2002.

Sanso', F. and C.C.Tscherning: Fast spherical collocation: A General Implementation. IAG Symposia, Vol. 125, pp. 131-137, 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. 115-123, 1993.