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 h₁ UW/DWC to level h₂ can be expressed as

(

1 2

) (

1 1

) ^{( )}

1

the gravity data obtained at level h₁. 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 altimeter-derived data and a density model are also considered in geoid modeling. In this study, the airborne gravity data is the most important; it has been discussed in detail in chapter 5. To evaluate the gravimetric geoid models, 38 high-quality GPS/leveling points are employed to assess the geoid accuracy.

3.2 Surface Gravity 3.2.1 Land Gravity

Land data (Fig 3-1(a)) were collected during 1980–2003 by Academia Sinica, Base Survey Battalion and Ministry of Interior (MOI), Taiwan (Yen et al., 1990; Yen et al., 1995; Hwang, 2001; Chen, 2003), using LaCoste&Romberg gravimeters (LCR, 1997) tied to some absolute gravity stations. These data were mainly measured along roads at intervals of 2 km between two observations and on geodetic control points.

The average data accuracy of Hwang (2001) and Chen (2003) are about 0.04 mgal;

they are both based on the adjustments of the relative gravity networks. The total number of land data is 3641. Most land gravity measurements are performed on the west plain. There are only a few gravity points over the Central Range due to the difficulty in performing the survey. Gravity anomalies over flat regions are moderate;

however, they become large over the high mountains, reaching values of approximately 200–300 mgal.

3.2.2 Shipborne Gravity

A part of the shipborne gravity data was surveyed by the National Central University (NCU) using the gravimeter R/Vl’ Atalante KSS30 in 1996 (Hsu et al., 1998) and the other part was obtained from the National Geophysical Data Center data set of the National Oceanic and Atmospheric Administration (NOAA), USA. In

21.2–25.8 N (Fig 3-1(b)). Most shipborne data are located over the Pacific Ocean and Bashi Channel. However, fewer data are located over the Taiwan Strait. The standard deviation of the crossover analysis by the NCU and the total shipborne data are 2.6 and 11.2 mgal, respectively. Some bad-quality shipborne data were removed and not subsequently used in geoid modeling. The total number of shipborne data after eliminating the outliers was 4084. There is an obvious local low near the eastern coast, reaching approximately –250 mgal, and the other shipborne data show moderate gravity anomalies.

3.3 Altimeter-Derived Gravity

Recently, altimeter-derived data has assumed more importance in marine geoid computations due to the major developments in satellites with altimetry missions and a rapid increase in the altimeter-derived data coverage. Although the altimeter-derived data usually provides lower accuracy than shipborne gravity data, it is sometimes more useful than shipborne data in geoid modeling; this is because obtaining a considerable amount of data for marine gravimetry is time-consuming and expensive.

In this study, we select the data from the KMS02 model for the geoid modeling investigations. KMS02 gravity field was modeled according to the GEOSAT mission and ERS using the DGM-E04 and JGM-3 orbit models (Anderson et al., 2003). The gravity model improved the quality and coverage of the altimetric height observations, particularly in the coastal regions. The region located between 119.2–122.8 E and 21.8–25.8 N was selected with a 2-min grid spacing (Fig 3-1(c)). Some outliers, particularly near the coast and over shallow water, were removed to enhance the geoid accuracy. As compared to shipborne gravity, the KMS02 data exhibited better coverage over the Taiwan Strait; therefore, it could compensate for the lack of shipborne gravity data and enhance the geoid accuracy over this area.

3.4 Geopotential Model

We adopted the EIGEN-GL04C coefficients for the computations of the long-wavelength geoid and gravity and of the error anomaly degree variances from 2°

to 360° in the LSC. The EIGEN-GL04C coefficients were determined by GFZ Potsdam and GRGS Toulouse. These coefficients were determined from the GRACE and LAGEOS missions (from 2003 to 2005) and the 0.5° × 0.5° gravimetry and

altimeter-derived data (GFZ, 2006); with regard to spherical harmonic coefficients, degree and order 360 up to a wavelength of 110 km was developed.

The EIGEN-GL04C model significantly improves the knowledge of the gravity field of the Earth. As compared to the other geopotential models (e.g., EGM96, EIGEN-CG01C, and EIGEN-CG03C), the EIGEN-GL04C model exhibits an improvement of approximately 3–10 cm in the geoidal heights obtained from the GPS/leveling points over USA, Canada, and Europe. Fig 3-2 shows the gravity anomalies and geoidal heights of EIGEN-GL04C up to degree and order 360 globally and over Taiwan. Both gravity anomalies and geoid heights over Taiwan significantly vary from 200 to –200 mgal and 12 to 28 m, respectively, because of the complex terrain. Both these parameters are higher over the Central Range and lower at the Ryukyu arc.

3.5 Digital Elevation Model

Three DEMs with different resolutions—9′′×9′′, 90′′×90′′, and 6′×6′—are used in the RTM investigation (Fig 3-3). Because bathymetry is not considered in this study, all the elevations at sea level in the three DEMs are zero. The 9′′×9′′ and

0 9 0

9 ′′× ′′ models, which are both considered to be true elevation surfaces, are applied to the inner and outer zone computations. The 6′×6′ model is considered to be the mean elevation surface. The reason for the division of the RTM computation task into two zones has been described in chapter 4.

All these DEMs were sampled from a high-resolution DEM, which is formed on a 3′′×3′′ grid (with a horizontal resolution of approximately 80 m) using photogrammetry by the Aerial Survey Office belonging to the Forest Bureau (Hwang et al., 2003a), Taiwan. The accuracy of the 3′′×3′′ DEM is approximately 4 m rms determined by comparing with the hundreds of benchmarks with precise elevations.

3.6 Density Model

The density data used in this study were provided by Chiou (1997). According to the distribution of rocks over Taiwan, the density data were obtained by associating each type of rock with an average density and stored in a 5′×5′ grid. The density model has been validated by reliable seismology data. Fig 3-4 shows a color map of the density over Taiwan. The densities are relatively low and are mostly below 2.0 g

cm^–3 over the west plain. Over the high mountains, the densities are much higher, and the highest density can reach approximately 3.0 g cm^–3. In Fig 3-4, the average density on land is 2.35 g cm^–3. Therefore, the rock density over Taiwan has an obvious variation and cannot be assumed to be the global density constant—2.67 g cm^–3.

3.7 GPS/Leveling Points for Evaluation

The general method for the evaluation of a geoid is a comparison with external data. Geoid height differences can be compared with the differences obtained from the GPS/leveling points. According to this method, the gravimetric geoid models are compared to the available GPS/leveling benchmarks with the observed geoidal heights. An observed geoidal height is the difference between the GPS-derived ellipsoidal height (from 24-h observations and at cm-level accuracy) and the precision leveling-derived orthometric height (at mm-level accuracy). Rigorous orthometric corrections have been incorporated into these GPS/leveling routes (Hwang and Hsiao, 2003; Hwang et al., 2007a). These GPS/leveling benchmarks can be divided into four routes (Fig 3-5). The north route is located at the northern coast of Taiwan; the east route lies in a valley, and the center and south routes are situated from the hills to the mountains and plains to the mountains, respectively. Geoid variation is moderate along the north (approximately 1 m) and east routes (approximately 3 m); however, it is considerable along the center and south routes (approximately 8 m).

Fig. 3-1 Distributions and free-air gravity anomalies of surface and altimeter-derived gravity. (a) Land data. (b) Shipborne data. (c) Altimeter-derived data. The total number of land, shipborne, and altimeter-derived data are 3641, 4084, and 10228, respectively.

Fig. 3-2 (a) Gravity anomalies and (b) geoid heights globally and over Taiwan obtained from the EIGEN-GL04C coefficients.

Fig. 3-3 DEMs used in the geoid modeling. The resolutions of the DEMs are (a) 9 s, (b) 90 s, and (c) 6 min. The elevations at the sea level are zero.

Fig. 3-4 Density model over Taiwan (unit: g/cm³). Data are stored in a 5-min grid.

The average density on land is 2.35 g/cm³.

Fig. 3-5 Four leveling routes for evaluating the geoid accuracy. Circles represent the benchmarks along the north leveling route, which lies along the north coast; stars denote the center route, which is spread from the hills to the high mountains; triangles represent the south route, which is located from the plains to the high mountains;

squares correspond to the east route, which lies at a valley. The colors denote the topography.

Chapter 4

RTM Effects in Geoid Modeling: Comparison of Three Methods

4.1 Introduction

We investigate three different methods—FFT, prism, and Gaussian quadrature—for the computation of RTM-derived effects in order to determine the most appropriate one for geoid modeling. Among these, the FFT method is a gridwise computation technique and the other two are pointwise computation methods. In the prism method, a density model of a topographic mass is taken into consideration.

4.2 RTM Effects by FFT

Although various methods are available for RTM-derived effects computation, the most commonly used method is the FFT technique due to its computational speed.

The main characteristic of this method is that it uses gridded information and returns the RTM-derived effect values for all the points on a grid. RTM-derived effects can be considered as the difference between two Bouguer reductions of true and mean topographic surfaces. Thus, the computation in this method requires at least two DEMs representing the two surfaces. The approximate expression for the RTM-derived effect on gravity can be expressed as follows (Forsberg, 1984):

(

p p

) (

ref

) (

p p

)

RTM x y G h h c x y

g , =2 − − ,

∆ π ρ (4-1)

where c

(

xp,yp

)

is the terrain correction at point P; h and h_ref, the elevations of the true and mean DEMs, respectively; G, the gravitational constant; and ρ, the mass density. c

(

xp,yp

)

can be computed in frequency domain. The terrain correction term in Eq (4-1) can be expressed in the convolution form as follows (Schwarz et al., 1990):

(

^{x y}

)

^G

[

^{h f h h f h g}

]

c _{p p} ² 2 _p( ) _p²

2

, = 1 ρ ∗ − ∗ + (4-2)

where 1₃

f = r , ⁼

∫

E dxdy g r1₃

, r= x² +y² , and E is the domain of integration in the X-Y plane. Further, h is the elevation of point P; E, the domain of integration in _p the X-Y plane; and ∗ , the convolution operator. If t₁ =h²∗ f and t₂ =h ∗ f , they can be expressed by Fourier transform as follows:

)

Subsequently, t₁ and t₂ can be obtained by inverse Fourier transform as follows:

(

^{( 2) ( )}

)

It is necessary to introduce the equation for deriving the RTM gravitational potential at point P to model RTM-derived geoid effects. This equation can be expressed as follows:

According to Bruns formula, which is given by N =T γ , the RTM-derived effect on geoid yields

where γ implies normal gravity. If we assume that the terrain effect on the geoid is

Eq (4-10) is a linear expression and its higher-order terms are ignored. It can be expressed in the convolution form as follows:

(

^{x y}

)

^G

(

^{h h}

)

^g can be obtained by inverse Fourier transform as follows:

(

^{x ,y}

)

^{G F 1}

⁽

^{F(h ) F(g)}

⁾

N_{RTM p p} = ⁻ _res ∗

γ

ρ (4-13)

On comparison with other algorithms, the obvious advantage of FFT is its rapid computation, but the unavoidable edge effects and cyclic convolution should be eliminated carefully by 100% zero-padding.

4.3 RTM Effects by Prism

The prism method is a simple technique to estimate RTM-derived effects by pointwise computation. The RTM-derived effects on geoid N_RTM and gravity

gRTM

∆ at point P due to the residual terrain mass can be expressed as

( ) _{∫∫ ∫}

These effects can be decomposed into a combination of many prisms. We can calculate the effect at point P(xp, yp, hp) as the sum of the mass of each prism (Fig 4-1) and add all the prism-derived effects within a selected zone. Thus, the RTM-derived gravity or geoid effects can be obtained. The equations for calculating these two effects by the prism method can be expressed as follows:

( ) _{∑ ∑}

true DEM and ∆ , its residual elevation. Residual elevation is the difference between z the mean and true elevations. The total number of prisms selected is n. In this study, in addition to the case where ρ is constant, we also consider the case where ρ is variable. Subsequently, Eq (4-16) and (4-17) become

( ) _{∑ ∑}

Because the prism method is a pointwise computation technique, it requires a considerable amount of time to complete the computational task. To reduce the time-consuming calculation, the most efficient strategy is to split the computational area into two parts—an inner zone with a fine elevation grid and an outer zone with a coarse elevation grid. Fig 4-2 shows the decomposed prisms of the inner and outer zones. The inner zone comprises thinner prisms from the detailed DEM and outer zone comprises thicker ones from the coarser DEM. Theoretically, the prism method is considered the most accurate, but it may be not suitable for the computation of high-resolution output due to its inefficiency.

Fig. 4-1 Geometry of the RTM-derived effects in the prism method. This method calculates the potential at point P due to the mass of a prism.

Fig. 4-2 Computational inner and outer zones at point P. The black thin and the grey thick prisms shown on the right hand side in the above figure belong to the inner and outer zones, respectively. The residual height denotes the difference between the elevations obtained from the true and mean DEMs.

∆x ∆y

∆ z

Mean elevation surface True elevation surface

P(x , _p y , _p h ) _p l

O(x, y, z)

4.4 RTM Effects by Gaussian Quadrature

The Gaussian quadrature method is a useful technique to obtain the integration of a function over a domain. The Gaussian quadrature theorem is based on a weighted sum of the function values at specified points within the integration domain. This method also belongs to pointwise computation methods. It was successfully employed by Hwang et al. (2003b) in the study of terrain correction.

Fig. 4-3 Geometry depicting an RTM-derived effect in the Gaussian quadrature.

According to Eq (4-7), the RTM-derived effects of gravity ∆g_RTM and geoid NRTM due to the topographic mass above and below a point at P(x_p,y_p,h_p) (Fig 4-3) can be derived as follows:

( ) _{∫∫ ∫}

₌⁼

( ) _{∫∫ ∫}

The RTM-derived geoid is only considered as a linear effect term which is the same as FFT. For a given area bounded by X1 (west), X2 (east), Y1 (south), and Y2 (north), Eq (4-20) and (4-21) can be numerically integrated as follows (Hwang et al., 2003b):

(

^,

)

^{2 ( , ) ( )} and M and N, the numbers of the weighting coefficients and nodes along the x and y axes over the domains [X₁,X₂]and[Y₁,Y₂], respectively (Press et al., 1989). To obtain the highest possible precision, M and N should be the numbers of the given grids along the x and y directions. The values of the function c(y) at nodes x and _i y were interpolated using the Newton-Gregory forward polynomial (Gerald and j

Wheatley, 1994) from the evenly spaced function values on a given grid. For the interpolations required in Eq (4-24) and (4-25), we experimented with various

polynomial degrees and found that the use of degrees higher than six yields no further improvement in the interpolation accuracy. This computation of one-dimensional case was proven to be successful by Press et al. (1989). However, Hwang et al. (2003b) have reported the successful usage of the required two-dimensional Gaussian quadrature.

The main advantage of using the Gaussian quadrature method is the very high-order accuracy it provides with fewer points. This is useful in the cases wherein a function requires a long time to compute by the pointwise method. However, it is also a time-consuming task. Therefore, the Gaussian quadrature method also requires the segregation of the computational area into inner and outer zones to make the computation more efficient. Practical tests reveal that the computation time required by the Gaussian quadrature method is more than that required by the FFT method but less than that required by the prism technique.

4.5 Design of Experiments

The primary objective of this study is to determine the most suitable method for computing RTM-derived effects in geoid modeling. The local gravity data used in this study include land and shipborne data. The altimeter-derived data are not considered.

In order to compare the three methods stated above, three geoid models whose RTM-derived effects are created by these methods are compared to the GPS/leveling points to assess their accuracies.

The geoid modeling procedure employed in this study, which is also based on the RCR procedure, is shown in Fig 4-4. In this figure, while carrying out the remove and restore steps during the computation of short-wavelength gravity and geoid, the three methods are taken into account individually. The process of geoid modeling is divided into four cases. The only discrepancy between these cases is that their RTM-derived effects are delivered by the FFT method (case 1), prism method with a constant ρ (case 2), prism method with a variable ρ (case 3), and Gaussian quadrature method (case 4).

The selected radii of the inner and outer zones and the sizes of the output grids of the RTM-derived effects computation are summarized in Table 4-1. Because the FFT method is a rapid computational technique, it does not require the usage of an outer zone in practical calculation. On the other hand, since cases 2~4 are based on the

prism and Gaussian quadrature methods, they require inner and outer zones to reduce the time required for computation. Further, in order to make the computations more efficient, the resolutions of the output grids in cases 2~4 have to be stored in coarser grids the resolution of which is set to 1 min in this study. With regard to the radii chosen for the inner and outer zones, the gravity effect will decay more rapidly than the geoid effect with the increase in the distance between the RTM mass and point P.

Thus, the RTM-derived geoid computation requires longer inner and outer zones’

radii than the RTM-derived gravity computation. In cases 2~4, the radii of the inner and outer zones for the gravity computation are 15 and 100 km, respectively, and for geoid computation, 30 and 300 km, respectively. In case 1, the radii of the inner zone for the gravity and geoid computations are 50 and 100 km, respectively.

Table 4-1 Radii of the inner and outer computational zones for the RTM-derived effects and the resolutions of the output grids in the four case models

Case Effect Radius of the inner zone

Radius of the outer zone

Output grid resolution

RTM gravity 50 km -

Case 1 RTM geoid 100 km - 9 s

RTM gravity 15 km 100 km

Case 2 RTM geoid 30 km 300 km 1 min

RTM gravity 15 km 100 km

Case 3 RTM geoid 30 km 300 km 1 min

RTM gravity 15 km 100 km

Case 4 RTM geoid 30 km 300 km 1 min

During computation, the grid sizes of the long-, residual-, and short-wavelength geoid parts are varied taking into account the different resolutions of the GGM, local gravity data, and DEM. The long- and residual-wavelength geoid effects are stored in the 3-min and 1-min grids, respectively. In order to make the grid sizes of the different wavelength geoid effects the same, we employ the GMT package (Wessel and Smith, 1995) to sample all the grids and add individual geoid effects to obtain the

During computation, the grid sizes of the long-, residual-, and short-wavelength geoid parts are varied taking into account the different resolutions of the GGM, local gravity data, and DEM. The long- and residual-wavelength geoid effects are stored in the 3-min and 1-min grids, respectively. In order to make the grid sizes of the different wavelength geoid effects the same, we employ the GMT package (Wessel and Smith, 1995) to sample all the grids and add individual geoid effects to obtain the

在文檔中 利用空載，地面重力與測高資料計算台灣大地起伏: 研究向上/向下延續與地形效應計算 (頁 36-0)