Chapter 1 Introduction
1.3 Outline of Thesis…
In chapter 2, the main principles of geoid modeling, including the spherical harmonic function, Stokes integration, and LSC, are presented. In addition, the RCR procedure and the theorem of UWC and DWC are also described.
All data except the airborne gravity data used in geoid modeling are introduced
in chapter 3. These data include the surface gravity, altimetry data, global geopotential model (GGM), DEM, density model, and GPS/leveling points for evaluating the geoid accuracy.
Geoid modeling results using three RTM-derived effects methods are presented in chapter 4. The three methods are FFT, prism and Gaussian quadrature. Besides, the influence of the density variation of the topographic mass on the geoid is also discussed.
In chapter 5, a description of the airborne gravimetry theorem and an airborne gravity survey over Taiwan are presented. Furthermore, three methods to evaluate the accuracy of airborne gravity data are mentioned.
In chapter 6, the application of airborne gravity data is discussed. The data are processed by the DWC technique and are used for the geoid computation. A comparison between the downward-continued airborne gravity data by the FFT and LSC is analyzed and described. Two low-pass filters used in the frequency domain and two types of geoid modeling by the LSC DWC are also investigated.
A summary, conclusions, future research, and suggestions are presented in the final chapter.
Fig. 1-1 Terrain and bathymetry around Taiwan (Hwang et al., 2007b). The inset is a tectonic map of Taiwan from Angelier et al. (1997). The Philippine Sea plate moves towards the Eurasia plate at a speed of 8.2 cm/year.
Chapter
2
Principles of Geoid Determination and Upward/Downward Continuations
2.1 Introduction
The strategy of geoid modeling used in this study is based on the remove-compute-restore (RCR) procedure, which is useful for high-resolution local gravity or geoid determinations. Geoid modeling takes into account information regarding three parts of the gravity field, namely, the long-, intermediate-, and short-wavelength parts. In this study, the long-wavelength part is determined from the global geopotential model by using a spherical harmonic function; the intermediate-wavelength part, from local gravity observations by using least squares collocation (LSC); and the short-wavelength part, from the high-resolution digital terrain model.
Upward/downward continuation (UW/DWC) is a method that can be used to transform the gravity potential on a surface into that on a higher/lower surface (Blakely, 1995). In other words, UWC and DWC are performed in order to obtain gravity functional from one level surface to another. It is important to apply both continuations to airborne gravity data in order to calculate the geoid by using gravity data at a different surface level.
2.2
Methodologies of Geoid Determination
On the global scale, the geoid can be represented in terms of a spherical harmonic expansion. On local and regional scales, a geoid model based on gravity can be obtained by using Stokes integration and LSC. The spherical harmonic representation and Stokes integration are deterministic, while LSC is stochastic.
2.2.1 Spherical Harmonic Representation of Gravity Field
According to Newton’s law of gravitation, the earth gravitation at point P can be expressed as (Fig 2-1) (Torge, 1989)
r dm
where r′ and r are the geocentric position vectors of the element mass dm and the attracted point P.
Fig. 2-1 Potential at point P due to the earth mass.
The corresponding potential V and the earth gravitation b have the relationships
V grad
b= (2-2)
Thus, the gravitational potential of the earth can be given by
rdv
where ρ and dv are the earth’s density and volume element, respectively.
0 lim =
∞
→ V
r . V is harmonic outside the spheroid and can be determined by using a spherical harmonic function given by (Torge, 1991)
))
where a is the semimajor axis of the ellipsoidal earth model and GM is the geocentric gravitational constant with respect to the total mass. λ , ϑ , and r are spherical coordinates and Cnm and Snm are fully normalized spherical harmonic coefficients, which are mass integrals that represent the mass distribution within the central body. Pnm is the associated Legendre function with degree n and order m.
The gravity anomaly and geoid undulation in the spherical harmonic function can be expressed as (Heiskanen and Moritz, 1967)
(
cos λ sin λ) (
sinφ)
where R is the radius of the earth. The long-wavelength features of the earth’s external gravity field are determined by using satellite gravimetry and are modeled as a series of solid spherical harmonics truncated at the maximum values of n and m. The spherical harmonic function is usually used along with the spherical harmonic coefficients to determine the global long-wavelength geoid or gravity field.
2.2.2 Stokes Integration
As shown in Heiskanen and Moritz (1967), the disturbing potential, T, can be determined by Stokes integration as
∫∫
∆Using Bruns’ formula, we can obtain the geoid undulation as
∫∫
∆=
σ
σ πγ Sψ gd
N R ( )
4 (2-9)
where γ denotes the normal gravity. In theory, Stokes integration can simply be calculated by using global gravity data coverage. However, in a geoid computation task, the RCR procedure is required in order to determine the geoid surface more accurately. Stokes integration is usually calculated rapidly in the frequency domain by using a fast Fourier transform (FFT) technique. On a sphere, a rigorous implementation of FFT can apply the spherical FFT or multi-band FFT technique (Forsberg and Sideris, 1993).
2.2.3 Least Squares Collocation
LSC can be used to determine an anomalous gravitational field by using different combinations of geodetic observations. The basic principle of LSC is given by (Moritz, 1980)
l D C C
s= sl( ll + )−1 (2-10)
where s and l are sets of signals and observations, respectively. C is the covariance ll matrix of l and C is the covariance matrix between s and l. D is the matrix of the sl noise vector, which functions as a filter and weight in LSC computations. To estimate the error of signal s, the error covariance matrix is computed as (Moritz, 1980)
ls ll sl ss
ss C C C C
E = − −1 (2-11)
where E denotes the error covariance matrix. In the case of geoid determination by ss LSC, the formulae of signal and error are
( )
C(
C D) ( )
gN = Ng gg + g −1 ∆ (2-12)
( ) ( )
⎟⎟ matrices for geoid-gravity and gravity-gravity, respectively. Dg is a diagonal matrix containing error variances of gravity anomalies. In (2-13), s is a scale factor, which can be determined bygg anomalies, c is error variance of the gravity anomalies derived from a geopotential gg model, and m is the number of gravity data points.
The covariance function provides the covariance between two signals, between two observations, and between a signal and an observation, and it is used in LSC to predict those signals that are of interest to us. The key factor for a precision geoid model by LSC is covariance functions. Thus, it is essential to find a suitable covariance function for use in LSC computation. In this study, for up to 360 degrees, we adopt the error anomaly degree variances of a geopotential model ; for higher degrees, we adopt the Tscherning-Rapp anomaly degree variance model 4. The Tscherning-Rapp model is generated from an empirical covariance function developed by Tscherning and Rapp (1974). The Model 4 of anomaly degree variance is expressed as
2
where n denotes the selected degree. A and B are both free parameters whose values are adopted to be 425.28 mgal2 and 24 in this study and RB is the radius of the Bjerhammar sphere. r and r′ are the distances of points P and Q from the earth’s center.We can determine the covariance between two points by using data obtained at different levels, such as airborne and surface gravity data. Eq (2-15) is used only for
>360
n . Based on the combination of the geopotential model and the Tscherning-Rapp degree variance model 4, the covariance functions between two gravity anomalies, between two disturbing potentials, and between a gravity anomaly and a disturbing potential for points P and Q can be expressed as
)
σ are the error variances between two gravity anomalies, between two disturbing potentials, and between a gravity anomaly and a disturbing potential, respectively; these error variances are associated with the corresponding geopotential model coefficients. P is the Legendre n polynomial of degree n and Ψ is the spherical distance between P and Q. PQ γ denotes normal gravity. More covariance functions such as those between two geoid gradients, between a geoid gradient and a gravity anomaly, and between a geoid gradient and a disturbing potential can be found in Tscherning and Rapp (1974).
The rapid developments in LSC over the past few years clearly demonstrate that LSC is being used as the primary technique for local geoid determination because it can accurately estimate the signals of interest to us by using heterogeneous data having different resolutions. Due to the multi-resolution characteristic of LSC, we select LSC for the primary geoid modeling methodology in this study.
2.3 Remove-Compute-Restore Procedure
The RCR procedure is one of the most well-known strategies used for regional geoid determination. The RCR procedure is also called the remove-restore technique.
In theory, geoid determination can only be performed for gravity data having a global coverage; however, a global gravity field model may represent data far beyond the area of interest. If the RCR procedure is used, gravity field data beyond the area of interest can be removed. In areas with complex topographies, it is very important to remove and subsequently restore the potential of the topography. For these areas, terrestrial gravity values are usually available locally at accessible spots; the remove procedure makes these values more smooth and representative. For many years, because of the valuable characteristics of the RCR procedure, considerable attention has been focused on the application of the RCR procedure to geoid modeling.
Furthermore, when performing geoid determination, long-wavelength and short-wavelength errors may arise if the RCR procedure is not properly applied.
The geoid and gravity field can be divided into three parts: long-wavelength (low-frequency), intermediate-wavelength (intermediate-frequency or so-called residual), and short-wavelength (high-frequency) parts. Therefore, both the height anomaly ζ and the gravity anomaly ∆ can be expressed as g
short res
long ζ ζ
ζ
ζ = + + (2-21) and
short res
long g g
g
g=∆ +∆ +∆
∆ (2-22)
where ζlong and∆glong are the long-wavelength height anomaly and gravity anomaly, respectively; ζshort and ∆gshort, the short-wavelength height anomaly and gravity anomaly, respectively; and ζres and ∆gres, the residual height anomaly and
gravity anomaly, respectively. Fig 2-2 shows geoid undulations at three different wavelengths. In the RCR procedure, the long- and short-wavelength parts are attributed to geopotential-derived and residual terrain model (RTM)-derived effects, respectively. Local gravity observations subtracted from the two gravity effects can be used in Stokes integration or LSC to determine the intermediate-wavelength geoid.
Subsequently, the geopotential-derived and RTM-derived geoid effects can be restored to obtain the final geoid.
Fig. 2-2 Three different wavelengths of geoid undulation. Nlong, N , and res Nshort denote the long-, intermediate- (residual), and short-wavelength parts of geoid undulation, respectively.
In this study, the long-wavelength gravity and geoid are based on a global geopotential model, and the intermediate geoid is obtained by local gravity data by LSC.
2.3.1 Long-Wavelength Reference Geopotential Model
The global geopotential model (GGM) is a model that can represent the earth’s potential field. This model is important for regional geoid determination because it takes care of the long-wavelength part of geoid.
For geopotential-derived gravity, the higher the degree and order used for the geopotential coefficients, the smaller is the area required with local gravity data, but errors in high-degree coefficients can be a problem if not carefully modeled. The factors influencing the accuracy of the GGM include the amount and quality of local gravity and satellite tracking data and the maximum degree of the model. In addition,
short res
long N N
N + +
res
long N
N +
Nlong
Ellipsoid
the GGM usually yields an absolute geoid height error (so-called long-wavelength error) of the order of a few decimeters due to biases. However, the relative geoid height is often accurate because the biases at two computational points will largely be canceled out when differential geoid height computations are performed.
2.3.2 Residual Terrain Model
The RTM represents the residual part between the true and mean elevation surfaces (Fig 2-3). For determining the short-wavelength geoid in high mountainous regions, it may be insufficient to use only the geopotential model and local gravity observations. This is due to the signal contribution of the topography, which is particularly strong at short wavelengths for a rough terrain. The effect of the RTM can represent these short-wavelength signals appropriately.
Fig. 2-3 Residual terrain model (RTM), which represents the difference between the true and mean elevation surfaces.
The RTM-derived effect can be expressed as the difference between two surface-derived effects. In a planar approximation, the potential of point P due to an RTM mass is
∫∫∫
− + − + −=
RTM x xp y yp z s
G dm
V 2 2 2
) ( ) (
) (
(2-23)
where dm is a mass element of RTM and (x , p y , p z ) and (x, y, z) are the p coordinates of point P and every mass element dm, respectively. It is important to use both the true and mean elevation surfaces in geoid computation. The true elevation surface should be represented by a digital elevation model (DEM) containing detailed
True elevation surface
Mean elevation surface P
RTM
information in order to take into account high-frequency signals. The mean elevation surface should be selected in such a manner that it represents the global distribution of the regionally varying signal characteristics as far as possible. The practical computational methods for the RTM-derived effects are described in chapter 4. Three such methods used for computing the effects are investigated.
2.4 Quasi-Geoid Correction
The difference between the geoid and a quasi-geoid is that the geoid corresponds to a datum of orthometric height and the quasi-geoid to that of normal height (Fig 2-4).
By considering the normal gravity gradient with respect to the surface of the mean reference ellipsoid, the quasi-geoid is defined as a function of the normal height (Vanicek et al., 1999). In practice, when orthometric heights are used for determining the vertical datum, a quasi-geoid correction is applied to the fundamental formula of physical geodesy in order to accurately determine the geoid.
Fig. 2-4 Physical surface of the earth. h, H , O H , N, and N ζ denote the ellipsoid height, orthometric height, normal height, geoid undulation, and quasi-geoid undulation, respectively.
The relationship between the height anomaly ζ and the geoid undulation N is expressed as (Heiskanen and Moritz, 1967)
g H
N B
ζ
− ≈ − ∆γ
(2-24) Sea levelGeoid
Ellipsoid Quasi-geoid
h
HO
HN
ζ N
Topography
where ∆gB is the Bouguer anomaly, γ is normal gravity, and H is the topographic height. Eq (2-24) can also be written as
2 2
G H
N
γ
ρ ζ
−π
≈ (2-25)
where ρ is the density of the terrain mass and G is the gravitational constant.
ρ πG
2 is the Bouguer term. The difference between the geoid and the quasi-geoid is minute over moderate topographies, but it can reach several decimeters over high mountainous areas. Thus, the quasi-geoid correction cannot be ignored over rough terrains.
2.5 Upward and Downward Continuations
UW/DWC is employed to calculate the potential at any point above/below a planar surface having a known potential. It is important to apply UWC and DWC to airborne gravimetry for assessing the quality of airborne gravity data and for computing geoid undulation. However, the characteristics of the two continuation operations are different. UWC is a smooth operation that is characterized as a well-posed problem, whereas DWC is an unstable operation that is characterized as an ill-posed problem.
An inverse problem is expressed as the solution of an operator equation by the following expression:
) A(m
d = (2-26)
where m is a function obtained from a metric space of model parameters, d is an element obtained from a metric space of data sets, and A is an operator. According to the classical theory of inverse problems, there are three definitions for well-posed and ill-posed problems (Zhdanov, 2002). A well-posed problem must satisfy the following conditions.
(1) Solution m of Eq (2-26) exists.
(2) Solution m of Eq (2-26) is unique.
(3) Solution m depends continuously on the left-hand side of the equation, i.e., on d.
The problem in Eq (2-26) is ill-posed if one of the three conditions fails. The gravity potential outside the mass of earth satisfies the Laplace equation
=0
∇g (2-27)
The gravity potential at some level z = 0 is assumed as
) called a DWC of the gravity potential. We can write an operator equation of the relationship between the potentials at z = h and z = 0 as
[ ]
where A is an operator used for calculating the UW/DWC of the gravity potential.
UWC is usually used to assess the accuracy of airborne gravity observations.
These airborne data can be compared with the surface gravity data that are upward continued to the flight altitude. DWC plays a key roll in geoid determination when using airborne gravity data. On the other hand, the estimation of downward-continued data is sensitive to noise. Therefore, some types of noise suppression operations are required to enhance the data quality.
In this study, two UW/DWC methods, FFT and LSC, are taken into consideration. Both methods have been applied to UWC and DWC for many years.
2.5.1 Continuation by Fast Fourier Transform
UW/DWC by the FFT method is based on the integral Poisson formula. If an airborne gravity survey is carried out at a constant altitude, DWC can be readily
implemented by using FFT in the frequency domain. Let the vertical component of
where z is the altitude of gravity field g. For the three-dimensional condition, the relationship between g(x,y,z=0) and g(x,y,z= can be written as (Buttkus, h)
We can use a convolution integral to represent Eq (2-32) as
) z = h plane. On the other hand, the two-dimensional Fourier transform is given by
dxdy (2-34), the corresponding wavenumber response function becomes
y dxdy can be expressed in the wavenumber domain as
0
In contrast to UWC, the wavenumber response function of DWC from the z = h plane to the z = 0 plane is given by
DWC by FFT is essentially a high-pass filtering operation that will amplify short-wavelength noise in data processing. Therefore, the DWC procedure used for airborne gravity is a very unstable process, and it will result in a rapid increase in noise, particularly at high flight altitudes. To reduce the noise, a filtering or smoothing technique should be applied to the FFT downward-continued method. Thus, Eq (2-38) becomes
In this study, UWC by FFT will be used to compare airborne and surface gravity
data, and DWC by FFT will be applied to geoid modeling. These investigations are described in chapter 5 and chapter 6, respectively.
2.5.2 Continuation by Least Squares Collocation
UW/DW C can also be performed by LSC in either the spectral domain or the spatial domain (Sideris, 1995). Although processing by LSC is not performed as rapidly as that by FFT, the advantage of LSC is that it provides a scheme that can combine airborne gravity data with surface gravity or other heterogeneous data. The equation for the case in which the gravity field at level h1 UW/DWC to level h2 can be expressed as
(
1 2) (
1 1) ( )
1the gravity data obtained at level h1. In this study,
2 determined by using the combination of GGM and the Tscherning-Rapp degree variance model 4. DWC by LSC in spatial domain will be used for investigating geoid modeling in chapter 6. Eq (2-41) is just one of the LSC downward continuation methods used in this study. Another method that involves direct use for geoid determination is also introduced in chapter 6.
Chapter 3
Data for Geoid Modeling
3.1 Introduction
The data used for geoid modeling in this study mainly include the local gravity data, GGM, and DEM. They are used in the calculation of the residual long-wavelength and short-wavelength gravity or geoid. In addition, the
The data used for geoid modeling in this study mainly include the local gravity data, GGM, and DEM. They are used in the calculation of the residual long-wavelength and short-wavelength gravity or geoid. In addition, the