Physical geodesy: gravity, geoid '' with basics of airborne gravity

1
Physical geodesy gravity, geoid .. with basics
of airborne gravity

• Lectures and computational GRAVSOFT exersizes
• OCTAS workshop
• Hurdal, Norway, Jan 2005
• by
• Rene Forsberg, Geodynamics Dept., Danish National
  Space Centre
• rf_at_spacecenter.dk
• Overview
• Basics of physical geodesy
• Basics of satellite and airborne gravity
• Geoid determination by Fourier methods and
  collocation
• Application examples Nordic area, Malaysia
• GRAVSOFT excersize Arctic Ocean geoid and ICESat

2
Basic physical geodesy
Background material Torge Gravimetry (de
Gruyter, 1989) Good condensed review Heiskanen
and Moritz Physical Geodesy, 1967 (Freeman out
of print). The basic book. Lecture Notes IAG
International Geoid Schools (1994-) IGeS,
Milano

Newtonian potential and gravity G
Gravitational constant ? density
Physical potential of earth W V ?, ? is
potential of centrifugal force ? ½ ?2(x 2 y
2) in geocentric coordinates
?²V 0 outside the masses Laplace
equation ?²V -4?G? inside masses Poisson
equation
3
Basic physical geodesy (2)

• Normal gravity field
• Find radial symmetric solution to Laplace
  equation which fullfill
• U(r,?,?) on reference ellipsoid (e.g., GRS80)
  constant U0
• U0 in principle selected so that U0 W0 ..
  average potential at sea-level
• Select closest reference ellipsoid to approximate
  equipotential surface of W corresponding to
  global average sea level the geoid
• Closed formulas and series expansions for normal
  potential and its derivatives in space Heiskanen
  and Moritz, chap. 2. Determined solely from (a,
  f, M, ? - semi-major axis, flattening, earth mass
  and rotation rate).

Example GRS-80 normal gravity ??GRS80
978032.7(1.0053024sin2?-.0000058sin22?) mgal,
on ellipsoid Conventional unit 1 mgal 10-5
m/s2
4
Basic physical geodesy (3)

• Anomalous potential (non-ellipsoidal potential)
• Full-fill Lapace equation ?2T 0 gt classical
  potential field theory can be used ..
• spherical harmonic expansions solid spherical
  harmonics r-nYnm(cos?) eigenfunctions ..
• boundary value problems gt integral equations
  (e.g., Stokes formula)
• Gravity field quantities become functionals of T
• Geoid
• Quasi-geoid
• Gravity anomaly Deflections
• of vertical

5
Geoid and heights
Geoid Actual Equipotential Surface
H Orthometric Height
Ellipsoid Reference Model Equipotential
Surface
Q
N Geoid Height
gQ
Unmodeled Mass
P
?p

• GRAVITY ANOMALY ? g g Q - ? P

• H (Orthometric Height) h (Ellipsoid Height)
  N (Geoid Height)

6
Definition of the geoid
Geoid Level surface of global undisturbed oceans
.. Complications Geoid is global
equipotential surface in geopotential (W
W0), but heights refer to local reference -
Mean Sea Level (MSL) land - Lowest Astronomic
Tide (LAT) sea The cm-geoid always refers
to a specific datum .. and is not necessarily a
true geoid
Height determination with GPS Horthometric
hellipsoidal N Whether to use geoid N or
quasigeoid depends on national height system
7

• Difference geoid and quasigeoid Nordic area

Approximative formula for Helmert orthometric
heights
8
Basic physical geodesy (4) spherical harmonics

• Expanding T in spherical harmonic functions
• Pnm are associated Legendre functions of degree n
  and order m
• Geoid heights multiply coefficients by 1/?.
• Gravity anomalies multiply by (n-1)/R
• Global solutions
• Satellite-only solutions (GRACE GGM02S, nmax
  160)
• Combination solutions (EGM96, nmax 360), uses
  global ½º average surface ?g-values
• Next-generation EGM EGM05, nmax 2160, due
  ultimo 2005 need global 5 ?g

9
EGM96 global geopotential solution
(Courtesy S. Kenyon)
10
EGM96 geopotential solution

• NIMA/NASA-GSFC joint effort
• 30x30 terrestrial gravity anomalies from the
  NIMA data archives
• 30x30 gravity anomalies derived from the
  GEOSAT Geodetic Mission altimeter data
• Satellite tracking to over 30 satellites
  (including new satellites tracked by SLR, TDRSS,
  GPS)
• Direct satellite altimetry (TOPEX/POSEIDON,
  ERS-1, GEOSAT) ocean gravity anomalies
  inverted from altimetry by least-squares
  collocation planar logarithmic covariance
  model Forsberg, 1987)
• Formal inversion solution for coefficients
• Evaluation of model by International Association
  of Geodesy (International Geoid Service, M. G.
  Sideris)
• See www.nga.mil or http//164.214.2.59/GandG/wgs-8
  4/geos.html

11
EGM96 surface gravimetry (1)
12
EGM96 surface gravimetry (2)
13
EGM ocean gravimetry satellite altimetry
14
EGM96 geoid errors
15
Satellite gravity missions CHAMP, GRACE, GOCE

• CHAMP CHAllenging Microsatellite Payload-
  German satellite, lauched 2000. Accelerometer,
  GPS, magnetometer ..
• GRACE Geopotential Reseach and Climate
  Experiment- NASA/GFZ satellite, launch 2002.
  Precision tracking between two spacecraft
  looking for gravity change (monthly solutions)
• GOCE Global Ocean Circulation Explorer- ESA
  mission, lauch 2006. Gravity gradiometer.

16

• GRACE science objectives
• Temporal change of earths gravity field
  monthly spherical hamonic solutionsUltra-precise
  geoid models to be used for
• deep ocean current changes (ocean bottom
  pressure),
• large-scale evapotranspiration,
• soil moisture changes,
• mass balance of ice sheets and glaciers,
• changes in the storage of water and snow of the
  continents,
• mantle and lithospheric density variations,
• postglacial rebound or
• solid Earth's isostatic response.

Monthly fields to n 100 need careful science
to avoid aliasing
17
GRACE Improvement of spherical harmonics
18
GRACE/CHAMP STAR accelerometer
Resolution 3 x 10-9 m/s2 range 10 3 m/s2
only for space use .. Proof mass 70 g / 4 cm
Star camera (TUD)
19
GRACE Ku-band range and range-rate

• Single horn antenna for transmission and
  reception of the dual-band (24 and 32 GHz) k- and
  ka-band microwave signals,
• Ultra-stable oscillator (USO) serving as a
  frequency reference,
• Microwave assembly for  up-converting the
  reference frequency , down-converting the
  received phase from the other satellite to
  approximately 2 MHz and for amplifying and mixing
  the received and the reference carrier phase
• Instrument processing unit (IPU) used for
  sampling and digital signal processing of the
  k-band carrier phase signals and the data of the
  GPS space receiver, the accelerometer (ACC) and
  the star camera assembly (SCA).

20
GRACE and GPS tracking

• GPS Space Receiver (GPS)
• The GPS TurboRogue Space Receiver receiver
  assembly provided by JPL serves for
• Precise Orbit Determination (POD) with
  cm-accuracy,
• Coarse Positioning (lt50m) for real time use by
  the atttitude/orbit control system
• Time Tagging of all payload data and
• Atmospheric and Ionospheric Profiling.

21
GRACE spherical harmonic models

• University of Texas
• - GGM02S static field, nmax 160 (real nmax
  ..120) 1 years data
• Monthly solutions from varying quality,
  improving with time
• GFZ/Potsdam Eigen-GRACE models 2-3 months data
• JPL Validation model 6 months data
• All available at http//podasc.jpl.nasa.gov

GRACE orbit change lower heights gt higher
resolution
22
GRACE geopotential secular change results

• Analysis of 22 available monthly GRACE spherical
  harmonic models Apr 2002 Aug 2004
• Models by CSR Texas provided to European GRACE
  science team
• Computations of monthly geoid (N) grids, with
  and without J2-term
• Grid regression for bias, trend and yearly
  seasonal terms (gracefit/GRAVSOFT)
• NGRACE a bt c
  cos(t) d sin(t)

Global geoid trend (b-term mm/yr), to spherical
harmonic degree 30
23
GRACE geopotential change 2002-2004 (cont.)
Secular trends in geoid without J2 left to
degree 30 right to degree 15. Colour scale 5
to 5 mm/year
Seasonal terms in geoid left amplitude
(mm/yr) right phase (fraction of year)
24
GRACE geopotential trends in Greenland

• Overall mass loss from GRACE epoch results
• - Convert gravity change to equivalent ice loss
• Restrict comparison to ice sheet area
  gt Current Greenland total mass loss 84 km3/yr

Estimate from glaciological model ice loss 96
km3 per degree warming Estimate from US PARCA
project 50 km3/yr (B. Thomas et al. -
equivalent global sea rise 0.13 mm/yr)
Geoid change from GRACE to degree 15 (-6 to 2
mm/yr)
25
US Parca project 1993-99 (Thomas et al.)
26
Example of geoid differences between GRACE
solutions
JPL GGM01S
JPL Eigen-Grace2
27
Validation surface data example northern Russia
5 mean-free-air data provided by Gleb Demianov,
Tsniigaik, for ArcGP project
28
Validation surface data example northern
Russia (2)
GPS-levelling data Moscow-St. Petersburg
Gravity comparison / lowpass filter statistics
table (mgal)
Difference observed gravity - GGM01S (mgal)
Difference JPL
29
Airborne gravity Arctic region
Arctic gravity project Compilation of airborne
gravity, submarine, surface data ..
30
Arctic region airborne gravity data validation
Naval Research Lab Arctic Ocean airborne data
1992-99
NRL 1992-99
KMS/Norwegian airborne data of Greenland/Svalbard/
North Atlantic 1998-2003
KMS/SK 1998-03
Typical airborne error estimates 1.5-2.0 mgal
r.m.s.
31
Integral Formulas space domain

• Stokes Formula
• Relating N with gravity observations.
• Stokes Kernel
• w is weight function used to limit influence of
  low harmonics

32
Stokes kernel modification controlling
long-wavelength errors

• Stokes kernel modification (i.e. modified
  Wong-Gore method without sharp edges)
• where
• The use of modified kernels means that the
  satellite information (e.g., from GRACE) is used
  for the long wavelengths, and the surface gravity
  data at shorter wavelengths
• Terrestrial gravity data often biased due to
  errors (e.g. lacking atmospheric correction),
  uncertainties in reference systems, aliasing due
  to terrain

33
Fourier Transform (1)

• Method requires planar approximation gridded
  data on same levelgt Spherical harmonic model to
  be removed before using Fourier methods
• Basic definition of 2-D Fourier transform
• kx and ky are called wavenumbers (like frequency
  in 1-D time domain) defined on infinite x-y
  plane

• Advantage of Fourier transforms convolution
  theorem
• Convolutions must faster in frequency domain
  than space domain many geodetic integrals can be
  expressed as convolutions

34
Fourier Transform (2)

• Derivates of Fourier transform
• Vertical derivates from upward continuation
  formula
• Anomalous potential relationships follows from
  these (allow the direct determination of geoid
  transform filter inverse transform!

35
Fourier transformation and
Stokes integral

• Stokes integral can conveniently be evaluated
  using FFT methods (Strang van Hess, 1990).
• This is convolution form if cos? is considered
  constant (simple spherical FFT) and the
  sin-formula is used for ? GRAVSOFT program
  SPFOUR
• Stokes formula in planar approximation gives

GRAVSOFT planar FFT program GEOFOUR
36
Fast Fourier Transform FFT

• Extremely fast algoritm Cooley and Tukey ca.
  1960
• 1-D FFT algoritm discrete Fourier transform
• Has special effectsGridded data gt spectrum
  periodic Data given on interval X gt data
  assumed periodic with period X Periodic data gt
  spectrum discrete wavenumber spacing 2?/X
• Approximation of continous transform with FFT

Data
FFT Spectrum negative wavenumbers at left due
to periodicity
37
2D-FFT considerations

• Zero padding minimize periodicity
  errors(drawback increase array size)
• Realize negative parts of spectrum folds by
  periodicity a real function will always have a
  radially symmetric spectrum

38
Spherical FFT - improvements

• Multiband spherical 2D FFT.
• Virtually exact method for handling data in very
  large areas.
• The spherical FFT is obtained by smoothly merging
  smaller bands of planar 2D-FFT solutions with
  different reference parallels (Forsberg and
  Sideris, 1993).
• Advantage of this method is that it is very fast,
  and that data only has to be transformed once.
• GRAVSOFT program spfour.for
•  
• One dimensional spherical Fourier transformation
• The exact one-dimensional spherical Fourier
  transformation devised by Haagmans et al (1993).
• FFT is only applied row-vice in the longitude
  direction along each fixed parallel (?l) in the
  grid with a summation in the latitude bands
• Disadvantage. Much slower than multiband 2D FFT.
• GRAVSOFT sp1d.for

39
2D-FFT need for interpolation

• GEOGRID (GRAVSOFT) can perform interpolation
  using collocation (mode 1) or weighted means
  (mode 2)
• Collocation second order Gauss-Markov covariance
  function
• s is the distance, C0 is the signal variance, a
  is the correlation length
• Weighted means prediction (power 2) quick and
  dirty
• In practice select closest neighbours (e.g.
  5/quadrant)

40
Geoid determination taking into account
the terrain

In mountainous areas a major part of the
variations of the gravity field is due to
terrain The Bouguer anomaly (topographically
reduced anomaly) is much more smooth than the
free-air anomaly ?g (c is the classical terrain
correction and typically much smaller than 2?G?h)
Sognefjord gravity and terrain
41
Global source of DEM data SRTM 3 and 30
data

• SRTM Global radarmission Feb 2000
• US/German/Italian mission to map
• global topography
• Global 30 SRTM DEM data
• General release by NIMA 2003-04

Manila Bay SRTM sample data
42
Remove-restore geoid determination
terrain

General principle compute mass effects for
terrain assuming density known
Complete (Bouguer) reduction of terrain not OK
for geoid .. Nterrain way too big ..
Alternative Residual terrain model (RTM)
smaller effects but removes height correlation
43
Terrain effect computations the prism building
block

• The rectangular prism of constant density is a
  useful "building block" for numerical
  integrations of the basic terrain effects
  gravity and geoid formulas

• Implementation in practice GRAVSOFT TC
  program
• - Input of height data Sequence of DEMs -
  Speed up prism formulas by approximative formulas
  at large distance - Supplement space domain
  formulas with Fourier domain formulas

44
Example of data reducton and geoid Nordic area
Gravity data from 12 countries Geodetic marine
gravimetry Oil company data Airborne gravity
surveys
45
The basic composite DEM (0.02 x 0.04 )
46
Reference DEM

• Filtering 79
• moving averages
• across 0.1º x 0.2º
• mean height grid
• Corresponds
• approximately
• to 75 km resolution.

47
Statistics of data reductions remove steps
48
Quasi-geoid from FFT (contribution of residual
gravity data)
20001600 FFT grid (100 zero padding), gridding
of reduced data by collocation
49
RTM terrain quasi-geoid effect (the restore
effect)
50
Final (quasi-)geoid sum of all parts
51
Another method Least Squares Collocation

• The ability in-corporate randomly spaced data of
  various types
• Predicting related geodetic quantities taking
  into account the different statistics of the
  input data.
• No interpolation is required.
• Csx , Cxx covariance matrices
• Covariance matrices Dxx contain the noise
  variance
• Error estimates may be obtained
• Very expensive computationally (solve large
  linear equation system)

52
Least Squares Collocation (2)

• Prediction of different types Covariance
  propagation need to know model for potential
  covariance function cov(T,T) K(s) .. all other
  functions follow from this

• Spherical earth The covariance between the
  anomalous potential is

siTT are called degree variances
53
Tsherning-Rapp model covariance functions
A, RB .. fitted to empirical covariances
. GRAVSOFT GEOCOL, EMPCOV, COVFIT
54
Self-consistent planar covariance model

• Planar domain ok for downward continuation of
  airborne data
• Requirement spatial analytical covariance
  function model- e.g. Tscherning-Rapp model-
  e.g. planar logaritmic model (Forsberg, 1987)

• Model fitted to empirical data by three
  parameters C0, D, T
• D corresponding to Bjerhammar sphere depth
• T is a long-wavelength attenuation compensation
  depth- Complete formulas for gravity, geoid,
  2nd order gradients in Forsberg (1987)
• GRAVSOFT GPCOL and GPFIT