Transcript Geoid determination by 3D least-squares collocation by C.C.Tscherning

Geoid determination by
3D least-squares collocation
by
C.C.Tscherning
Niels Bohr Institute, University of Copenhagen
C.C.Tscherning, University of Copenhagen, 2008-09-10
1
Purpose:
Guide to gravity field modeling, and especially to geoid
determination, using least-squares collocation (LSC).
DATA

GRAVITY FIELD MODEL

EVERYTHING
=
Height anomalies, gravity anomalies, gravity
disturbances, deflections of the vertical, gravity
gradients, spherical harmonic coeffients
C.C.Tscherning, University of Copenhagen, 2008-09-10
2
Quasi-geoid:
Important:
the term geoid = the quasi-geoid,
i.e. the height anomaly evaluated at the surface of
the Earth.
C.C.Tscherning, University of Copenhagen, 2008-09-10
3
Gravsoft
The use of the GRAVSOFT package of FORTRAN
programs will be explained.
A general description of the GRAVSOFT programs
are given in:
Forsberg, R. and C.C.Tscherning: An overview
manual for the GRAVSOFT Geodetic Gravity
Field Modelling Programs. 2.edition. Contract
report for JUPEM, 2008. http
C.C.Tscherning, University of Copenhagen, 2008-09-10
4
Gravsoft
• We will only consider 3D applications.
C.C.Tscherning, University of Copenhagen, 2008-09-10
5
FORTRAN 77
• All programs in FORTRAN77.
• Have been run on many different computers under
many different operating systems.
• Available commercially, but free of charge if used
for scientific or educational purposes.
C.C.Tscherning, University of Copenhagen, 2008-09-10
6
General methodology
• General methodology for (regional or local)
gravity field modelling :
• Transform all data to a global geocentric geodetic
datum (ITRF05/GRS80/WGS84), (but NO tides,
NO atmosphere) GEOCOL
• “geoid-heights” must be converted to height
anomalies N2ZETA
• Use the remove-restore method.
C.C.Tscherning, University of Copenhagen, 2008-09-10
7
Remove-restore method
• The effect of a spherical harmonic expansion and
of the topography is removed from the data and
• subsequently added to the result. GEOCOL, TC,
• TCGRID
• This is used for all gravity modelling methods
including LSC.
• This will produce what we will call residual data.
C.C.Tscherning, University of Copenhagen, 2008-09-10
8
Covariance
Determine a Reproducing Kernel so that it agrees
with a covariance function for the residual data in
the region in question.
EMPCOV, COVFIT
Thereby we have an analytic representation of the
covariance function.
C.C.Tscherning, University of Copenhagen, 2008-09-10
9
Select
Make a homogeneous selection of the data to be
used for geoid determination using rule-ofthumbs for the required data density, SELECT
X
o
o
X
x o
xo
o
o
x
x
If many data select those with the smallest error X
Selection of points O closest to the middle. 6 points
selected
C.C.Tscherning, University of Copenhagen, 2008-09-10
10
Errors
• check for gross-errors (make histograms and
contour map of data), GEOCOL
• verify error estimates of data, GEOCOL.
C.C.Tscherning, University of Copenhagen, 2008-09-10
11
Gravity field approximation and datum
– Determine using LSC a gravity field approximation,
including contingent systematic parameters such as
height system bias N0. GEOCOL
– Compute estimates of the height-anomalies and their
errors. GEOCOL
– If the error is too large, and more data is available add
new data and repeat.
–
C.C.Tscherning, University of Copenhagen, 2008-09-10
12
Restoring and checking.
• Check model, by comparison with data not used to
obtain the model. GEOCOL.
• Restore contribution from Spherical Harmonic
model and topography. GEOCOL, TC.
• Convert height anomalies to geoid heights if
needed N2ZETA.
• The whole process can be carried through using
the GRAVSOFT programs
• Compare with results using other methods !
C.C.Tscherning, University of Copenhagen, 2008-09-10
13
Test Data and programs
• GRAVSOFT includes data from New Mexico,
USA, which can be used to test the programs and
procedures. (Files: nmdtm, nmfa, nmzeta etc.)
• They have here been used to illustrate the use of
the programs.
• If pyGravsoft has been loaded on your computer,
programs are found in directory src, executables in
bin and data in data. Documentation in doc.
C.C.Tscherning, University of Copenhagen, 2008-09-10
14
Anomalous potential.
• The anomalous gravity potential, T, is equal to the
difference between the gravity potential W and the socalled normal potential U,
T = W-U.
• T is a harmonic function, and may as such be
approximated arbitrarily well by a series in solid
spherical harmonics, Snm
N
~
T ( ,  , r )  GM 
n

n  2 m  n
•
Cnm S nm ( ,  , r )
GM is the product of the gravitational constant and the mass of the Earth and C nm the fully
normalized spherical harmonic coefficients.
C.C.Tscherning, University of Copenhagen, 2008-09-10
15
Coordinates used.
GEOCOL accepts geocentric, geodetic and Cartesian
(X,Y,Z) coordinates but output only in geodetic.
C.C.Tscherning, University of Copenhagen, 2008-09-10
16
Solid spherical harmonics.
S nm ( P)  S nm ( ,  , r ) 
 cos m ,0  m  n
a
n  1 Pnm (sin  ) 
r
 sin| m|  , n  m  0
• where a is a scale factor (generally the semi-major
axis) and Pnm the Legendre functions.
• We have orthogonality:


1
S ( ,  , R) S ( ,  , R) cos( )d d


4  
n
/2
nm
ij
 /2
 1, n  i , m  j
 
 0, n  i  m  j
C.C.Tscherning, University of Copenhagen, 2008-09-10
17
Bjerhammar-sphere
The functions Ynm(P) are orthogonal basefunctions in a Hilbert
space with an isotropic innerproduct, harmonic down to a
so-called Bjerhammar-sphere
totally enclosed in the Earth.
T will not necessarily be an
element of such a space, but may be approximated arbitrarily
well with such functions. In spherical approximation the
ellipsoid is replaced by a sphere with radius 6371 km.
C.C.Tscherning, University of Copenhagen, 2008-09-10
18
Spherical approximation
• NOT used when evaluating an EGM.
• When used: r=Mean radius+h.
• Geodetic latitude = geocentric latitude.
C.C.Tscherning, University of Copenhagen, 2008-09-10
19
Reproducing Kernel

K ( P, Q) 

n 2



n 2
a 
n  
 rr ' 
2
(2n  1)a 2 n
n

m  n
S nm ( P) S nm (Q)
n 1
Pn (cos )
r
P
ψ
r’
Q
where ψ is the spherical distance between P and Q,
Pn the Legendre polynomials and σn are positive
constants, the (potential) degree-variances.
C.C.Tscherning, University of Copenhagen, 2008-09-10
20
Inner product, Reproducing property
1

2
(2n  1)a  n

S nn ( P), S nn ( P)

Ynn
2
T ( P)  (T (Q), K ( P, Q)) 

n
n 2
m  n
 




i
i2
ji
GM Cnm S nm (Q),  (2i  1)a 2 i  Sij ( P)Yij (Q)
n
 
n 2
m  n

n
n 2
m  n
 

i
GM Cnm S nm ( P) (2i  1)a 2 i 
i2
j i


Sij (Q) Sij (Q)

GM Cnm S nm ( P)  T ( P)
C.C.Tscherning, University of Copenhagen, 2008-09-10
21
Closed expression – no summation to

• the degree-variances are selected equal to simple
polynomial functions in the degree n multiplied by
exponential expressions like qn, where q < 1, then
K(P,Q) given by a closed expression. Example:
2 n 1
 RB 
 n   2   qn 1
R 

K (P, Q) 

n 0
R
2
B
R 
 
R 
2
B
2
n 1
R 
 
 rr' 
2
n 1
Pn (cos  ) 
1
R 4B  (rr') 2  2rr' R 2B
C.C.Tscherning, University of Copenhagen, 2008-09-10
22
Hilbert Space with Reproducing Kernel
• Everything like in an n-dimensional vector
space.
GM
• COORDINATES: a 2 Cnn
• ANGLES  between two
f , g
cos(  ) 
gf
• functions, f, g
• PROJECTION f ON g:
• IDENTITY MAPPING:
C.C.Tscherning, University of Copenhagen, 2008-09-10
g
 f , g
g
T(P)  T(Q), K(P, Q)
23
Data and Model
In a (RKHS) approximations T from data for which the
associated linear functionals are bounded.
• The relationship between the data and T are expressed
through functionals Li,
yi = Li (T) + ei + A  X
T
i
yi is the i'th data element,
Li the functional, ei the error,
Ai a vector of dimension k and
X a vector of parameters also of dimension k.
C.C.Tscherning, University of Copenhagen, 2008-09-10
24
Data types
GEOCOL codes:
11
12
13
16
17
• Also gravity gradients,
• along-track or area
mean values.
C.C.Tscherning, University of Copenhagen, 2008-09-10
25
Test data
GRAVSOFT data from the New Mexico Area to be used in LSC geoid examples.
Type
Format (see Forsberg & Tscherning,
Error
File
2008)
name
Free-air gravity
(Fig. 2)
#, latitude, longitude, altitude, Δg
0.2 mgal nmfa
Deflections
#, latitude, longitude, ξ, η
0.5
arcsec
nmdfv
Height anomalies
#, latitude, longitude, altitude, ζ
0.02 m
nmzeta
5m
nmdtm
Digital height model Grid-label, data in N-S, E-W
C.C.Tscherning, University of Copenhagen, 2008-09-10
26
Linear Functionals, spherical approximation
T
g = r
T 2T
g = r r
1 T
 =r 
1
T
=r cos(  ) 
C.C.Tscherning, University of Copenhagen, 2008-09-10
27
Best approximation: projection.
Ti pre-selected base functions:
C.C.Tscherning, University of Copenhagen, 2008-09-10
28
Collocation
LSC tell which functions to
select if we also require that
approximation and
observations agree and
how to find projection !
Suppose data error-free:
We want solution to agree with
data
We want smooth variation
between data
C.C.Tscherning, University of Copenhagen, 2008-09-10
29
Projection
Approximation to T using error-free data is obtained
using that the observations are given by, Li(T) = yi
C.C.Tscherning, University of Copenhagen, 2008-09-10
30
LSC - mathematical
• The "optimal" solution is the projection on the ndimensional sub-space spanned by the so-called
representers of the linear functionals, Li(K(P,Q)) =
K(Li,Q).
• The projection is the intersection between the
subspace and the subspace which consist of
functions which agree exactly with the
observations, Li(g)=yi.
C.C.Tscherning, University of Copenhagen, 2008-09-10
31
Collocation solution in Hilbert Space
~
T (Q) 
n

bi K ( Li , Q)
i 1
Normal Equations
{bi } = {K Li , L j  } { y j }
-1
Predictions:
n
L(Tˆ) =  bi K( Li , L)
i=1
C.C.Tscherning, University of Copenhagen, 2008-09-10
32
Statistical Collocation Solution
We want solution with smallest “error” for all
configurations of points which by a rotation around the
center of the Earth can be obtained from the original
data. And agrees with noise-free data.
~
Li (T )  Li (T )
We want solution to be linear in the observations
~ ( P) 
T
n

i 1
n
 i yi 

i 1
C.C.Tscherning, University of Copenhagen, 2008-09-10
 i L i (T)
33
Mean-square error - globally
~ (P)  T(P) 
T
n
  L (T)  T(P)
i 1

1
2 
8  

/2
 / 2

2
0
i
i
n
(T(P)    i L i (T))
2
i 1
d PPi cos( )d d
C.C.Tscherning, University of Copenhagen, 2008-09-10
34
Global Covariances:
CPi  COV(T(P), L i (T)) 
1
2  T( P) L i ( T) cos  d d d
8
C0 (T(P))  COV(T(P), T(P)) 
1
2
T
(
P
)
cos  d d d
2 
8
Cij  COV(L i (T), L j (T)) 
1
2  L j ( T) L i ( T) cos  d d d
8
C.C.Tscherning, University of Copenhagen, 2008-09-10
35
Covariance – series development
COV(P, Q)  COV(T(P), T(Q) 


i 2
2 i 1
R 
 i   Pi (cos  ),
 rr'
2
 GM 
i   
Cij  , DEGREE  VARIANCES
j i R
i
C.C.Tscherning, University of Copenhagen, 2008-09-10
36
Collocation Solution
~
T ( P) 
n

i 1
bi COV (T ( P), Li (T ))
b   COV (( L (T ), L (T ))  y 
1
i
i
j
i

   COV (( L (T ), L (T ) COV (T ( P), L (T ))
1
j
i
j
C.C.Tscherning, University of Copenhagen, 2008-09-10
i
37
Noise
• If the data contain noise, then the elements σij of the variancecovariance matrix of the noise-vector is added to K(Li,Lj) =
COV(Li(T),Lj(T)).
• The solution will then both minimalize the square of the norm of T
and the noise variance.
• If the noise is zero, the solution will agree exactly with the
observations.
• This is the reason for the name collocation.
• BUT THE METHOD IS ONLY GIVING THE MINIMUM
LEAST-SQUARES ERROR IN A GLOBAL SENSE.
C.C.Tscherning, University of Copenhagen, 2008-09-10
38
Minimalisation of mean-square error
The reproducing kernel must be selected equal to the
empirical covariance function, COV(P,Q).
This function is equal to a reproducing kernel with the
degree-variances equal to
 
2
n
n

m _ n
 GM  2  R 
 2  Cnm  
 R 
a
2
C.C.Tscherning, University of Copenhagen, 2008-09-10
2n 2
39
Parameter estimate
• A simultaneous estimate of T and of the
parameters X are obtained as
T
1
T
~
T ( P)  {C Pi } C ( y  A X )
T 1
1
T 1
~
X  ( A C A  W ) ( A C y)
• where W is the a-priori weight matrix for the
parameters (Generally the zero matrix).
C.C.Tscherning, University of Copenhagen, 2008-09-11
40
Error-estimates
• The associated error estimates are with
H  {COV(L, Li )}T C 1
• the mean square error of the parameter vector
2
mX
 (A C
T
1
A  W)
1
• and the mean square error of an estimated
quantity L(T~ ).
m2L   L2  H{COV ( L, Li )}  HAm2X ( HA) T
C.C.Tscherning, University of Copenhagen, 2008-09-11
41
N0 estimation
• Height system bias may be determined.
• The determination of a datum-shift requires data
which covers a large area,
• If the area is not large,
this will be reflected in large
2
error estimates m X .
• This will be illustrated using 2920 gravity data and 20
height anomalies from the New Mexico test area.
C.C.Tscherning, University of Copenhagen, 2008-09-11
42
Covariance Propagation
• The covariances are computed using the "law" of
covariance propagation, i.e.
• COV(Li,Lj) = Li(Lj(COV(P,Q))),
• where COV(P,Q) is the basic "potential"
covariance function.
• COV(P,Q) is an isotropic reproducing kernel
harmonic for either P or Q fixed.
C.C.Tscherning, University of Copenhagen, 2008-09-10
43
Covariance of gravity anomalies
ev P (T)  T(P)
T
2
 g(P)  
 ev P (T)
P
r
r
Appy the functionals on
K(P,Q)=COV(P,Q)

  2
COV( g(P), g(Q)) =  - - ev P  

 r r
 
  2
- evQ    i
 i=2
 r  r '
i+1
 R2 
  Pi cos 
 r r 
i -1  R 

 
=  - - ev   
  P cos 

 r
r '  r r 
(i -1 )  R 
= 
  P ( cos )
r r  r r 
I
2
i
P
i+1
i
i=2
2

i
2
i+1
i
i=2
C.C.Tscherning, University of Copenhagen, 2008-09-10
44
Evaluation of covariances
The quantities COV(L,L), COV(L,Li) and
COV(Li,Lj) may all be evaluated by the sequence
of subroutines COVAX, COVBX and COVCX
which form a part of the programs GEOCOL and
COVFIT.
C.C.Tscherning, University of Copenhagen, 2008-09-10
45
Remove-restore.
If we want to compute height-anomalies from
gravity anomalies, we need a global data
distribution.
If we work in a local area, the information outside
the area may be represented by a spherical
harmonic model. If we subtract the contribution
from such a model, we have to a certain extend
taken data outside the area into account.
(The contribution to the height anomalies must later be
restored=added).
C.C.Tscherning, University of Copenhagen, 2008-09-10
46
Change of Covariance Function
C.C.Tscherning, University of Copenhagen, 2008-09-10
47
Homogenizing the data
We obtain a minimum mean square error in a very
specific sense:
• as the mean over all data-configurations which by
a rotation of the Earth's center may be mapped into
each other.
• Locally, we must make all areas of the Earth look
alike.
• This is done by removing as much as we know,
and later adding it. We obtain a field which is
statistically more homogeneous.
C.C.Tscherning, University of Copenhagen, 2008-09-10
48
Homogenizing (II)
• 1.We remove the contribution TEGM from a known
spherical harmonic expansion like the EGM96 to N=360
or EGM2008 to N=2190.
• 2. We remove the effect of the local topography, TM, using
Residual Terrain Modelling (RTM): Earths total mass not
changed,
• but center of mass may have changed !!!
• We will then be left with a residual field, with a
smoothness in terms of standard deviation of gravity
anomalies between 50 % and 25 % less than the original
standard deviation.
C.C.Tscherning, University of Copenhagen, 2008-09-10
49
Residual quantities
yir  yi  Li (TEGM )  Li (TM ) 
Li (T )  Li (TEGM )  Li (TM )  ei  A X
T
i
C.C.Tscherning, University of Copenhagen, 2008-09-10
50
Exercise 1.
• Compute residual gravity anomalies and residual
height anomalies using the EGM96 spherical
harmonic expansion and the New Mexico free-air
gravity anomalies in data/nmfa and the height
anomalies in data/nmzeta.
• The Python module GEOEGM may be used to
subtract the contribution from EGM96. Result
files: data/nmfa-egm96.dat and data/nmzetaegm96.dat.
C.C.Tscherning, University of Copenhagen, 2008-09-10
51
Input-file: data\nmfa
C.C.Tscherning, University of Copenhagen, 2008-09-10
52
GEOEGM with nmfa
C.C.Tscherning, University of Copenhagen, 2008-09-10
53
Output-files
Output from run on screen (last part) and in geoegm.log:
http://cct.gfy.ku.dk/geoidschool/appendix2.txt
COMPARISON OF PREDICTIONS AND OBSERVATIONS
DATA TYPE = 13
NUMBER: 2920
OBSERVATIONS PREDICTIONS DIFFERENCE
MEAND
9.181986
12.113828
-2.931842
ST.DEVI.
30.405342
23.100298
21.283325
MAX
162.500000
77.947778 126.430437
MIN
-58.700000 -28.049925 -74.792095
DISTRIBUTION OF DIFFERENCES, UNITS: 5.000000
17 21 40 66107169225261317309324279221149110 73 59 29 25 24 19
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 OUTSIDE
76
GEOCOL TERMINATED AT:
Wed Jul 23 11:18:40 2008
C.C.Tscherning, University of Copenhagen, 2008-09-10
54
Output-file contoured with GMT.
C.C.Tscherning, University of Copenhagen, 2008-09-10
55
Exercise 2. Residual topography removal
The RTM contribution may be computed and
subtracted using the TC module.
– First a reference terrain model must be constructed
using the program TCGRID with the file data/nmdtm
as basis,
– The result should be stored in files with names nmfaegm96-tc.dat.
C.C.Tscherning, University of Copenhagen, 2008-09-10
56
Residual topography removal
C.C.Tscherning, University of Copenhagen, 2008-09-10
57
Residual topography removal
C.C.Tscherning, University of Copenhagen, 2008-09-10
58
Smoothing or Homogenisation
C.C.Tscherning, University of Copenhagen, 2008-09-10
59
Consequences for the statistical model.
• The degree-variances will be changed up to the
maximal degree, N, sometimes up to a smaller
value, if the series is not agreeing well with the
local data (i.e. if no data in the area were used
when the series were determined).
• The first of N new degree-variances will depend
on the error of the coefficients of the series. We
will here suppose that the degree-variances at least
are proportional to the so-called error-degreevariances,
C.C.Tscherning, University of Copenhagen, 2008-09-10
60
Error-degree-variances

err
l
 GM 
 

 a 
2
l
 (
ml
EGM
lm
a
)  
R
2l  2
2
The scaling factor α must therefore be determined from
the data (in the program COVFIT, see later).
Note that the error-degree variances refer to the mean –
earth sphere.
C.C.Tscherning, University of Copenhagen, 2008-09-10
61
Covariance function estimation and representation.
The covariance function to be used in LSC is equal to
COV(P,Q) =
2 /2 2
1
T(P)  T(Q) d cos  d d



8

• where α is the azimuth between P and Q and φ, λ are
the coordinates of P. The spherical distance ψ is fixed.
• This is a global expression, and that it will only
dependent on the radial distances r, r' of P and Q and
of the spherical distance ψ between the points.
2
0 - /2 0
C.C.Tscherning, University of Copenhagen, 2008-09-10
62
Global-local evaluation
• In practice it must be evaluated in a local area by
taking a sum of products of the data grouped
according to an interval i of spherical distance,
 i - /2  < i + /2
• Δψ is the interval length (also denoted the
sampling interval size).
C.C.Tscherning, University of Copenhagen, 2008-09-10
63
Covariance
• In spherical approximation we have already
derived
•
COV( g(P), g(Q))
R 
 n - 1
= 
 n  
 r r 
n=2  R 

2
2
n+2
Pi ( cos )
• where R is the mean radius of the Earth.
C.C.Tscherning, University of Copenhagen, 2008-09-10
64
Exercise 3. Hand calculation of covariances.
The following data must be used, with format: number, latitude, longitude,
altitude and gravity anomaly in mgal.
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
56.0
56.1
56.2
56.3
56.4
56.5
56.6
56.7
56.8
56.9
57.0
57.1
57.2
57.3
57.4
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
4.0
2.0
0.0
-2.0
-4.0
-6.0
-8.0
-9.0
-7.0
-5.0
-3.0
-1.0
1.0
5.0
4.0
C.C.Tscherning, University of Copenhagen, 2008-09-10
65
Exercise 3.
• Compute the empirical gravity anomaly
covariance function using the program EMPCOV
for the New Mexico area both for the anomalies
minus EGM96 and for the anomalies from which
also RTM-effects have been subtracted (input files
nmfa-egm96.dat and nmfa-egm96-tc.dat).
C.C.Tscherning, University of Copenhagen, 2008-09-10
66
Exercise 3.
C.C.Tscherning, University of Copenhagen, 2008-09-10
67
Exercise 3: Part of empcov.log.
2920 VALUES INPUT FROM FILE data/nmfa-egm96-tc.dat
NUMBER OF OBS 1= 2920 MEAN = 0.3066 VAR. = 173.53
HISTOGRAM, USING BIN SIZE= 5.000
0 0 0 0 3 9 34 49118248372887357343233127 59 35 18 18 10 0 0
OUT-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9' 10OUT
PSI COVA( 1, 1) PROD. STDV OF COV..
O M (UNIT)**2 NUMB (UNIT)**2
0 0.00 175.509014
3132
4.9 0.330449
0 2.00 145.042029
4762
3.4 20.163033
0 4.00 119.317243
9487
2.3 45.463833
0 6.00 90.458667
14003
1.8 75.421001
0 8.00 67.998329
17995
1.5 101.052568
0 10.00 48.805606
21731
1.3 119.723639
0 12.00 31.054065
25167
1.2 136.020523
0 14.00 18.040267
28540
1.1 149.844158
C.C.Tscherning, University of Copenhagen, 2008-09-10
68
Empirical and analytic Covariances
200.0
Empirical and analytic fitted covariance functions for New Mexico
free-air gravity minus EGM96 and residual topography
150.0
Empirical covariance
Covariance (mgal**2)
Analytic covariance
100.0
50.0
0.0
-50.0
0.00
0.20
0.40
Spherical distance (degrees)
C.C.Tscherning, University of Copenhagen, 2008-09-10
0.60
69
Degree-variances
We see here, that if we can find the gravity anomaly
degree-variances, we also can find the potential
degree variances.
However, we also see that we need to determine
infinitely many quantities in order to find the
covariance function
C.C.Tscherning, University of Copenhagen, 2008-09-10
70
Model-degree-variances
• Use a degree-variance model, i.e. a functional dependence
between the degree and the degree-variances.
• In COVFIT, three different models (1, 2 and 3) may be
used. The main difference is related to whether the
(potential) degree-variances go to zero like n-2, n-3 or n-4.
The best model is of the type 2,
•
2n2
A
 RB 
 


(n  1)( n  2)( n  B)  R 
2
n
• where RB is the radius of the Bjerhammar-sphere, A is a
constant in units of (m/s)2, B an integer.
C.C.Tscherning, University of Copenhagen, 2008-09-10
71
COVFIT
• The actual modelling of the empirically
determined values is done using the program
COVFIT. The factors a, A and RB need to be
determined (the first index N+1 must be fixed).
• Instead of A the gravity variance C0 at zero is
used.
• The program makes an iterative non-linear
adjustment for the Bjerhammar-sphere radius, and
linear for the two other quantities
C.C.Tscherning, University of Copenhagen, 2008-09-10
72
Divergence ?
Unfortunately the iteration may diverge (e.g. result
in a Bjerhammar-sphere radius larger than R).
• This will normally occur, if the data has a very
inhomogeneous statistical character.
• Therefore simple histograms are always produced
together with the covariances (in EMPCOV) in
order to check that the data distribution is
reasonably symmetric, if not normal.
C.C.Tscherning, University of Copenhagen, 2008-09-10
73
Exercise 4.
Compute using COVFIT an analytic representation
for the covariance function.
Input are: empirical error-degree variances from
EGM96 in data/egm96.edg
Covariance table from EMPCOV:
data/covtable_nmfa.txt
Mean altitude and variance of gravity anomalies
Area boundaries and data spacing.
C.C.Tscherning, University of Copenhagen, 2008-09-10
74
Exercise 4.
C.C.Tscherning, University of Copenhagen, 2008-09-10
75
Exercise 4 - bottom of covfit.log.
TAU(J) USED IN THE CX MATRIX 0.10E+01 0.10E+01
0.10E+01
RESULTS IN VARIANCE OF GRAVITY ANOMALIES:
1'TH ROW OF INVERSE MATRIX 0.4132E-01 -0.1243E-01 0.3080E-01
2'TH ROW OF INVERSE MATRIX -0.1243E-01 0.1211E-01
0.4449E-01
3'TH ROW OF INVERSE MATRIX -0.3080E-01 0.4449E-01
0.1782E+00
STD.DEV.
0.576779E-01 0.732173E+05 0.334610E+03
STD.DEV.*RMS 0.329951E-01 0.418845E+05 0.191416E+03
RESULTS IN VARIANCE OF GRAVITY ANOMALIES: 334.36
MGAL**2.
• N
RATIO AA
A
RE-RB VARG IT
• 360 0.5721D+00 0.2837 0.6654D+06 -792.72 334.36 10
C.C.Tscherning, University of Copenhagen, 2008-09-10
76
Table of model covariances at h=1700 m.
Ψ (deg)
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
1.10
1.20
ζ,ζ
Δg,Δg
Δg,ζ
ξ,ξ
0.0476
0.0463
0.0430
0.0387
0.0342
0.0298
0.0260
0.0201
0.0167
0.0150
0.0141
0.0133
0.0120
0.0102
0.0081
0.0061
174.15
139.00
90.20
53.43
27.60
10.08
-1.15
-10.61
-9.71
-4.53
0.91
4.40
5.20
3.75
1.11
-1.51
2.058
1.885
1.535
1.167
0.837
0.566
0.358
0.112
0.038
0.062
0.115
0.152
0.152
0.117
0.061
0.003
3.878
2.749
1.307
0.318
-0.309
-0.678
-0.860
-0.857
-0.593
-0.264
0.004
0.153
0.174
0.100
-0.014
-0.121
C.C.Tscherning, University of Copenhagen, 2008-09-10
ξ,Δg
0.000
-11.167
-13.884
-13.132
-11.188
-8.886
-6.597
-2.730
-0.233
0.885
0.939
0.377
-0.353
-0.910
-1.117
-0.961
ξ,ζ
0.000
-0.092
-0.146
-0.167
-0.167
-0.153
-0.132
-0.084
-0.044
-0.021
-0.014
-0.019
-0.028
-0.036
-0.038
-0.035
77
Conversion of geoid heights to height anomalies
• In 3D LSC all data must refer to points outside the
surface of the Earth.
• Geoid heights must be converted to height
anomalies.
• GRAVSOFT module N2ZETA can be used. It
implements the equation:
N P0   P 
g B
0
C.C.Tscherning, University of Copenhagen, 2008-09-10
H
78
LSC geoid determination from residual data.
• We now have all the tools available for using
LSC: residual data and a covariance model.
•
•
•
•
1.establish the normal equations,
2.solve the equations, and
3. compute predictions and error estimates.
This may be done using GEOCOL.
C.C.Tscherning, University of Copenhagen, 2008-09-10
79
Equations
• However we have to solve a system of equations as large
as the number of observations. GEOCOL has been used to
handle 100000 observations simultaneously.
• This is one of the key problems with using the LSC
method. The problem may be reduced by using means
values of data in the border area.
• Globally gridded data can be used (sphgric)
C.C.Tscherning, University of Copenhagen, 2008-09-10
80
Necessary data density (d)
• Function of correlation length of the covariance function.
• We want to determine geoid height differences with an
error of 10 cm over 100 km. This corresponds to an error
in deflections of the vertical of 0.2".
• This is equivalent to that we must be able to interpolate
gravity anomalies with2
2
(d
0.3
/
)



C
0
e
1
d
a mean error of 1.2
mgal. The
where ed is the st.dev. ,
rule-of-thumb is
C 0 the gravity variance,
 1 the correlatio n distance
C.C.Tscherning, University of Copenhagen, 2008-09-10
81
Exercise 5. Data density.
• Use the residual gravity variance C0, and the
correlation distance determined in exercise 4 for
the determination of the needed data spacing.
• Then use the program SELECT for the selection
of points as close a possible to the nodes of a grid
having the required data spacing, and which
covers the area of interest.
C.C.Tscherning, University of Copenhagen, 2008-09-10
82
Data selection.
• The area covered should be larger than the area in
which the geoid is to be computed. Data in a
distance at least equal to the distance for which
gravity and geoid becomes less than 10 %
correlated, cf. the result of exercise 3.
• When data have been selected (as described in
exercise 5) it is recommended to prepare a contour
plot of the data. This will show whether the data
should contain any gross-errors. LSC may also be
used for the detection of gross-errors.
C.C.Tscherning, University of Copenhagen, 2008-09-10
83
Exercise 6.GEOCOL INPUT.
An input file for the program GEOCOL must then
be prepared,
In order to compute height-anomalies at terrain
altitude, a file with points consisting of number,
latitude, longitude and altitude must be prepared.
This may be prepared using the program GEOIP and
GLIST, and input from a digital terrain model
(nmdtm). A file data/nm.h2 has been prepared.
C.C.Tscherning, University of Copenhagen, 2008-09-10
84
Exercise 6.
• Prepare a file named nm.h covering the area
bounded by 33.0o and 34.0o in latitude and -107.0o
and -106.0o in longitude consisting of sequence
number, latitude, longitude and height given in a
grid with 0.1 degree spacing.
• Use the program GEOIP with input from nmdtm.
This will produce a grid-file. This must be
converted to a standard point data file (named
nmh2) using the program GLIST.
C.C.Tscherning, University of Copenhagen, 2008-09-10
85
GEOCOL INPUT/SPECIFICATIONS.
• the coordinate system used (GRS80),
• the spherical harmonic expansion subtracted (and
later to be added),
• the constants defining the covariance model
• the input data files (nmfa-egm96-tc.data and
nmzeta-egm96-tc.dat)
• the files containing the points in which the
predictions should be made (nmzeta-egm96tc0.dat and nm.h2).
C.C.Tscherning, University of Copenhagen, 2008-09-10
86
GEOCOL OPTIONS
• Several options must be selected such as whether
error-estimates should be computed or whether we
want statistics to be output.
• The Cholesky-reduced normal equations can be
re-used if the program is run with the same
observations as in an earlier run.
C.C.Tscherning, University of Copenhagen, 2008-09-10
87
Exercise 7.
• Run the program GEOCOL (geocol17) with the
selected gravity data for the prediction of height
anomalies and their errors in the points given by
nmzeta.
• Compare input and predicted values. Remember
to use negative data codes (-13 and -11) since
EGM96 has already been subtracted.
C.C.Tscherning, University of Copenhagen, 2008-09-10
88
Exercise 8.
• Make a second run where the reduced height
anomalies are used as a second data-set.
• Re-use the 2920 reduced normal equations.
• This will determine the height-system bias N0. (If
the Python module GEOCOL is used, see
following).
C.C.Tscherning, University of Copenhagen, 2008-09-10
89
C.C.Tscherning, University of Copenhagen, 2008-09-10
90
Exercise 7.
In geocol.log we see the height system bias N0:
ELEMENTS OF (AT*C**-1*A)**-1.
0.0085
CORRELATION MATRIX:
1.0000
PARAMETER TYPE ESTIMATE
ERROR
ESTIMATE (FOR TILT: ZERO POINT
1
11111 -0.783996623 0.092367316
C.C.Tscherning, University of Copenhagen, 2008-09-10
91
Exercise 7.
C.C.Tscherning, University of Copenhagen, 2008-09-10
92
Exercise 7.
C.C.Tscherning, University of Copenhagen, 2008-09-10
93
Exercise 7+8 continued.
In geocol.log is a comparison of the input height anomalies
with the predicted.
COMPARISON OF PREDICTIONS AND OBSERVATIONS
DATA TYPE = 11
NUMBER:
20
OBSERVATIONS PREDICTIONS DIFFERENCE
MEAND
-0.897275
-0.897275
0.000000
0.184852
ST.DEVI.
0.159149
0.156162
0.015303
0.000101
MAX
-0.632800
-0.635039
0.038347
0.185066
MIN
-1.267500
-1.246579
-0.028575
0.184632
DISTRIBUTION OF DIFFERENCES, UNITS: 0.010000
0 0 0 0 0 0 0 1 3 1 10 2 2 0 1 0 0 0 0 0 0
0
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 OUTSIDE
.
C.C.Tscherning, University of Copenhagen, 2008-09-10
94
Exercise 9. RESTORE.
When the LSC-solution has been made, the RTM
contribution to the geoid must be determined.
Use tc with the file nm.h defining the points of
computation.
The LSC determined residual geoid heights and the
associated error-estimates are shown above.
C.C.Tscherning, University of Copenhagen, 2008-09-10
95
Exercise 8.
• Compute the RTM contribution to the geoid using tc and
add the contribution to the output file from exercise 7,
nm.geoid.
• If mean gravity anomalies, deflections or GPS/levelling
determined geoid-heights were to be used, they could
easily have been added to the data.
• It would not be necessary to recalculate the full set of
normal-equations.
• Only the columns related to the new data need to be added.
Likewise, an obtained solution may be used to calculate
such quantities and their error-estimates.
C.C.Tscherning, University of Copenhagen, 2008-09-10
96
Error detection:
• LSC filters the data .
• We may as done in geodetic network
adjustment inspect the residuals by using LSC
for the prediction and comparison with the data
used to determine the model.
• Large difference – possible gross error.
Outlier
Predicted
C.C.Tscherning, University of Copenhagen, 2008-09-11
97
Exercise 10.
Detection of suspected gross-errors may be done by
comparing the differences between observed
quantities and predicted quantities to the
estimated error. Use GEOCOL (cf. Exercise 6) for
this purpose by predicting reduced gravity
anomalies (data/nmfa-egm96-tc0.dat) and
comparing these with values predicted from an
identical file, but named data/nmfa-egm96-tc.dat.
A file name for a file to hold suspected gross
errors must be input.
C.C.Tscherning, University of Copenhagen, 2008-09-10
98
Conclusion (I)
• We have now went through all the steps from data
to predicted geoid heights.
• The exercises describes the use of gravity data and
height anomalies,
• and the determination of a height-system bias.
C.C.Tscherning, University of Copenhagen, 2008-09-10
99
Conclusion (II)
• The difficult steps in the application of LSC is the
estimation of the covariance function and subsequent
selection of an analytic representation.
• The flexibility of the method is very useful in many
circumstances, and is one of the reasons why the method
has been applied in many countries.
• If the reference spherical harmonic expansion is of good
quality, only a limited amount of data outside the area of
interest are needed in order to obtain a good solution.
•
C.C.Tscherning, University of Copenhagen, 2008-09-10
100
Conclusion (III)
• But if this is not the case, data from a large
border-area must be used so that a vast
computational effort is needed to obtain a
solution.
• This may make it impossible to apply the method.
• A way out is then to use the method only for the
determination of gridded values, which then may
be used with Fourier transform techniques or Fast
Collocation.
C.C.Tscherning, University of Copenhagen, 2008-09-10
101