Results of Geoid Modeling

4.6.2 Results of Geoid Modeling

In the residual geoid computation by LSC, the covariances of gravity-gravity anomaly, geoid-gravity anomaly, and geoid-geoid anomaly based on the combination of the EIGEN-GL04C model and Tscherning-Rapp degree variance model 4 are shown in Fig 4-7. The covariances containing shorter spherical distances have higher values. The patterns of the covariance values tends to be mild when the spherical distance exceeds 0.4º. In addition, the variations in the quality of the gravity data must be taken into account in order to determine the noise for different data types. We assign 0.1 and 1.0 mgal data noise to land and shipborne gravity anomalies, respectively. These values are empirical and yield the best results.

The residual gravity anomalies obtained by subtracting the EIGEN-GL04C- and RTM-derived effects from the original gravity anomalies are shown in Fig 4-8. The residual geoid effects obtained by LSC are shown in Fig 4-9. In these two figures, four distinct local low can be observed over the northern, central, and southern regions on land and the eastern region at sea. The minimum value of the residual gravity anomalies and geoids can reach over –150 mGal and –1.5 m, respectively. On the other hand, in the case where higher RTM-derived gravity anomalies are obtained,

smaller residual gravity anomalies are observed. The same result is obtained for the RTM-derived and residual geoids in the four cases considered in this study.

After restoring long and residual wavelength parts of the geoid and considering quasi-geoid correction, we can obtain the final geoid models for the four cases (Fig 4-10). Fig 4-10(a) obviously contains more detailed signals (or noises) over the Central Range due to the higher-resolution grids (9 s). Fig 4-11 shows the differences between the case 1 geoid model and the others. In most areas, the differences are less than ±0.1 m, but there exists a large difference area over the east coast, especially in Fig 4-11(c), reaching 0.5 m. Other large differences occur in the Central Range, which contains complex geoid variations.

Table 4-4 lists the statistics of the differences between the observed and modeled geoidal heights at the four leveling routes. The standard deviations of these differences in the four cases on the north and east routes are all within 8 cm. Most standard deviations on the center and south routes where the current land gravity data are sparsely distributed and the geoid variation is large are over 10 cm. However, we can conclude that case 1 presents the best accuracy as compared to the other cases as it offers the smallest standard deviation, especially in the center route; the standard deviation of 0.144 m in this case is 3~4 cm better than that in the other cases. A comparison of the mean values of case 2 and case 3 listed in Table 4-4 reveals that the maximum difference between the geoid surfaces obtained in these two cases is 4 cm.

However, their standard deviations are within 1 cm. This implies that the density variation considered in geoid modeling results in only a very limited improvement in the accuracy of the geoid over the Taiwan Island.

The FFT method is thus more suitable for geoid modeling over Taiwan due to its computational speed and the accuracy of the results obtained in this study, which was the best of all the three methods considered. In the subsequent investigations on geoid modeling, the RTM-derived effects are all delivered by the FFT method.

Table 4-4 Statistics regarding the differences (m) between the observed and modeled geoidal heights at four leveling routes

Method Leveling route Max Min Mean Std. dev.

North –0.034 –0.205 –0.130 0.062

East –0.186 –0.399 –0.312 0.068

Center –0.195 –0.577 –0.350 0.144

Case 1

South –0.330 –0.476 –0.379 0.046

North –0.046 –0.265 –0.154 0.071 East –0.213 –0.454 –0.321 0.080 Center –0.208 –0.775 –0.435 0.176 Case 2

South –0.235 –0.524 –0.399 0.139

North –0.027 –0.210 –0.138 0.067

East –0.294 –0.529 –0.365 0.082

Center –0.200 –0.799 –0.444 0.180

Case 3

South –0.227 –0.539 –0.381 0.149 North –0.140 –0.351 –0.249 0.070 East –0.292 –0.576 –0.433 0.079 Center –0.308 –0.884 –0.543 0.195 Case 4

South –0.233 –0.496 –0.377 0.141

Fig. 4-5 RTM-derived gravity anomalies. (a) FFT method (case 1), (b) prism method:

constant density (case 2), (c) prism method: variable density (case3), and (d) Gaussian quadrature method (case 4).

Fig. 4-6 RTM-derived geoids. (a) FFT method (case 1), (b) prism method: constant density (case 2), (c) prism method: variable density (case 3), and (d) Gaussian quadrature method (case 4).

Fig. 4-7 Covariances for (a) surface gravity-surface gravity covariance matrix, (b) geoid-surface gravity covariance matrix, and (c) geoid-geoid covariance matrix.

Fig. 4-8 Residual gravity anomalies. (a) FFT method (case 1), (b) prism method:

constant density (case 2), (c) prism method: variable density (case 3), and (d) Gaussian quadrature method (case 4).

Fig. 4-9 Residual geoids. (a) FFT method (case 1), (b) prism method: constant density (case 2), (c) prism method: variable density (case 3), and (d) Gaussian quadrature method (case 4).

Fig. 4-10 Geoid models obtained in (a) case 1, (b) case 2, (c) case 3, and (d) case 4.

Fig. 4-11 Geoid differences between the geoid models of (a) case 1 and case 2, (b) case 1 and case 3, and (c) case 1 and case 4.

Chapter 5

Airborne Gravity Data of Taiwan

5.1 Introduction

Airborne gravimetry can be used over regions with sparse gravity data coverage.

As mentioned in chapter 1, Taiwan’s terrain is complex, and it would be difficult to conduct a land gravity survey over its mountainous regions. In order to bridge the gaps in the existing ground gravity coverage (mainly in inaccessible areas), the Ministry of the Interior (MOI) of Taiwan sponsored an airborne gravity survey over the period from May 2004 to May 2005. The entire project including the survey work and software development was carried out by National Chiao Tung University, Taiwan, and National Survey and Cadastre (KMS), Denmark. The survey area covers the entire Taiwan Island and its surrounding seas.

5.2 Data Reduction in Airborne Gravity 5.2.1 Gravity Reduction

Airborne gravimetry can be classified into two types—scalar and vector types.

The scalar-type gravimetry employs a relative gravimeter and can be implemented in airborne and oceanic surveys. However, the vector-type component measurements employ the accelerometer of inertial measurement units (IMU), whose accuracy is generally lesser than that of the scalar-type but suitable for gravity measurements of 3D gravity components. The basic equation for a vector-type gravimetry is expressed as follows (Schwarz and Li, 1996):

(

2Ω P

)

v f v

g= & − + − (5-1)

where g represents the vector of gravity, and v and v& represent the vectors of velocity and acceleration, respectively, of the aircraft. f represents the specific force sensed by an IMU. Ω and P are both skew-symmetric matrices, which are expressed as (Olesen, 2003)

⎥ ⎥

where ω represents the angular velocity of the earth’s rotation in an inertial frame, and λ& , the angular velocity along the east–west direction of an aircraft in an Earth-fixed Cartesian frame. The term

(

2Ω+P

)

v is the Eto&&vo&&s effect (Harlan, 1968), which is attributed to the difference between the angular velocity of an aircraft and that of stationary objects. Thus, the gravitational attraction exerted on the aircraft slightly increases or decreases when the aircraft moves along the east or west direction. Eq (5-1) can be separated into three components; these components are expressed as north–south, and vertical directions of the aircraft. v& , _e v& , and _n v& represent the _u accelerations along the east–west, north–south, and vertical directions of the aircraft.

f , e f , and _n f are specific force components as observed by a gravimeter. _u v and _e v represent the velocities along the east–west and north–south directions,

respectively, of the aircraft. R and R represent the radii of curvatures along the _{N M} meridian and prime vertical, respectively. φ and h represent the latitude and flight height of the aircraft.

Because the airborne gravimeter is placed in a strap-down system, the navigation frame is not the same as the frame of the gravimeter sensors. The sensor frame should be converted into the navigation frame. This means that we need to consider the small difference between the horizontal accelerations recorded by the gravimeter and the GPS measurements. Thus, the rotation of specific forces in Eq (5-1) into the local coordinate frame results in (Olesen, 2003)

(

2Ω P

)

v Rf v

g= & − + − (5-7)

where R represents the rotation matrix that converts the sensor frame into a navigation frame. This transformation yields

⎥⎥ the gravimetry sensor and navigation system. c represents the conversion factor between f and q. If this correction is considered, Eq (5-6) becomes

⎟⎟

(

^sin ⋅ fx −^{cos sin} ⋅fy −⁽¹−^{cos cos )}⋅ fz

)

^}

represent the Eto&&vo&&s effect and tilt correction of the vertical components, respectively.

5.2.2 Aircraft Positioning

The position, velocity, and acceleration of the aircraft play an important role in airborne gravity surveys (Schwarz and Li, 1997, and Kennedy, 2002). The GPS positioning for airborne gravity used in this survey not only provides a precise flight trajectory position of the aircraft but, more importantly, precisely estimates the first and second derivatives with respect to time for computing the velocity and acceleration required for airborne gravimetry data processing. It is possible to achieve an accuracy of the order of centimeters for the aircraft position (Goad and Yang, 1997). Highly accurate velocities and accelerations can be obtained based on these positions by the use of a precise numerical technique. In this study, the trajectories of the aircraft are determined by using Bernese 5.0 (Beutler et al., 2004) with the IGS precise ephemeris of the GPS.

In the kinematic positioning using Bernese, a number of parameters included in the double-differenced phase observations were estimated together with the aircraft position. These parameters are grouped into two subsets in normal equations (Hwang et al., 2007b):

subset x₁ contains the ground station coordinates, tropospheric parameters, and phase ambiguities, while subset x₂ contains the epoch-by-epoch kinematic positions of the aircraft. Subset x₂ can be inverted as follows:

)

To obtain this solution, the standard stochastic model of the GPS-phase observables was used (Seeber, 2003). In this case, the double-differenced phase observables between the aircraft and the eight tracking stations are selected and used to obtain the final coordinates. The initial values of the kinematic positions (parameter subset x₂) are required for the linearization of the nonlinear GPS observation equations.

5.3 A Taiwan Airborne Gravity Survey 5.3.1 Survey Campaign

The survey lines are shown in Fig 5-1(a). These lines consist of 64 north–south, 22 east–west, 10 northeast–southwest, and 6 northwest–southeast oriented lines with a spacing of 4.5 km, 20 km, 5 km, and 30 km, respectively. The west–east and northwest–southeast lines are mainly used for crossover analyses. The survey area covers the whole of Taiwan Island and its offshore regions. The survey area is approximately 75,000 km² and the total distance covered by the survey lines is approximately 53,000 km.

A scalar-type gravimeter called LaCoste and Romberg (LCR) Air-Sea Gravity System II (serial number: S-133) (Fig 5-1(b)) mounted on a laser gyro-stabilized platform is used to record the airborne gravity data at 1 Hz. This gravimeter has a resolution of 0.01 mgal and an accuracy of 1 mgal, as seen from the shipborne test (LCR, 2003). It uses spring tension and beam velocity measurements to obtain the relative gravity variations. Additional information on the Air-Sea Gravity System II gravimeter is summarized in Table 5-1. The gravimeter is placed in a medium-size

aircraft, Beechcraft-200 (Fig 5-1(c)), flying at an average altitude of 5156 m (The spatial resolution is approximately 6 km (Torge, 1989)) at a mean ground speed of approximately 306 km/h. Both the airborne gravimeter and aircraft belong to the Ministry of the Interior, Taiwan.

The King-Air Beechcraft-200 is equipped with a Trimble 5700 GPS receiver (Fig 5-1(d)) that samples data at 2 Hz. For the kinematic positioning of the aircraft, eight ground-based GPS reference stations (Fig 5-1(a)) around Taiwan are used to determine the kinematic solutions. The eight stations are YMSM, SNAM, KDNM, PKGM, TMAM, FLNM, KMNM, and MZUM. The sampling rate of these reference stations is 2 Hz except that of SNAM station, which is 1 Hz. CCK, shown in Fig 5-1(a), is the Taichung airport, where aircraft take off and landing occurred.

The reference gravity value can be determined by a land gravimeter based on the absolute gravity reference points at the Taichung FG5 (Micro-g, 1999) absolute gravity station. The gravity value at the aircraft parking spot was recorded using a Graviton-EG gravimeter (LCR, 2002). The standard error of this gravity value is 0.04 mGal based on the relative gravity network adjustment. A number of gravity base readings of the airborne gravity system need to be obtained during the field survey period to obtain a smooth drift of the airborne gravimeter.

The airborne gravity survey was carried out from May 2004 to May 2005. The survey took 43 days, including 3 days of re-flights where bad data were found. The number of flight hours exceeds 200.

Table 5-1 Overview of the L & R Air-Sea Gravity System II Resolution 0.01 mGal

Accuracy <1.00 mGal Size 71 × 56 × 84 cm

Weight 116 kg

Power supply 240 W (avg), 450 W (max) Sampling rate 1 Hz

(a) (b)

(c) (d)

Fig. 5-1 (a) Airborne gravity survey lines and GPS tracking stations (solid circles) for precise aircraft positioning. The star represents the Taichung (CCK) airport, where the King-Air Beechcraft-200 is based. (b) The L&R Air-Sea Gravity System II gravimeter and (c) the King-Air Beechcraft-200 aircraft; the circle denotes the antenna. (d) Inside of the King-Air Beechcraft-200; the L&R Air-Sea Gravity System II and Trimble 5700 are mounted inside the aircraft.

5.3.2 Data Processing

Kinematic GPS solutions are obtained by using the combination of eight different GPS based stations and processed using Bernese 5.0 (Beutler et al., 2004) and the IGS precise ephemeris ( http://igscb.jpl.nasa.gov/ ). In order to determine the velocity and acceleration of the aircraft, we use the program DERIV of the International Mathematical and Statistical Library (IMSL) to perform the numerical differentiations. DERIV first computes the spline interpolants to the input functions (i.e., coordinate components x, y, and z) and then differentiates the spline interpolants to obtain their derivatives (Hwang et al., 2006b). Following the GPS procedure, two data processing techniques have to be considered: correction for time shift and filtering of raw gravity observations.

The time systems of the raw GPS and the gravimeter observations are inconsistent. The gravimeter time associated with a gravity reading is obtained from the clock of the computer attached to the gravimeter. Therefore, the gravimeter is not synchronous with the GPS clock and requires correction. In order to synchronize the two time systems, the time series of the raw gravity reading and vertical aircraft acceleration can be used (Olesen, 2003). Because gravity signals are much smaller than the vertical accelerations of the aircraft in common weather conditions, most raw readings recorded by the gravimeter are those of the vertical accelerations of the aircraft. Thus, the patterns of the raw gravity readings of the gravimeter and vertical acceleration readings of the GPS receiver are very similar. According to this characteristic, the shift between these two time series can be determined using a correlation analysis (Hwang et al., 2006b).

It is necessary to use an along track filter for raw airborne gravity data containing considerable noise due to turbulence. We use a Gaussian filter with a filter width of 150 s to eliminate high-frequency signals. The chosen filter width is a trade-off between noise reduction and gravity signal preservation and is proved to be the optimum width in the latter part of this study.

A description of the software that fulfils many of the requirements summarized above can be found in Shih (2004), Hwang (2005), and Hwang et al. (2006b and 2007b). The procedure for the airborne data processing is summarized in Fig 5-2.

Fig. 5-2 Flow chart of the airborne gravity data process implemented by NCTU. f , _x f , and n f are the three components of the gravimeter measurements. _u

ϕ

, λ, and h represent the latitude, longitude, and ellipsoidal height. ϕ_{& ,} λ& , h&, and ϕ&&, λ&&, h&&

represent the three components of the velocities and accelerations of the aircraft, respectively. g is the output gravity at the flight altitude.

Airborne gravimeter

GPS receiver

f ,x f ,_n f _u

ϕ

λ h

‧Bernese software

‧IGS precise ephemeris

ϕ_& λ& h&

ϕ&& λ&& h&&

Numerical differential

g

‧Time-shift adjustment

‧Eto&&vo&&s correction

‧Tilt correction

5.4 Results of the Airborne Gravity Survey

Fig 5-3 shows the free-air gravity anomalies at 5156 m using a Gaussian filter with a width of 150 s. The free-air gravity anomalies vary at an altitude of 5156 m over Taiwan and around the sea ranges from approximately –200 mGal over the east trench to 300 mGal over the high mountains (Fig 5-3). The values in the range of ±50 mGal nearly vary over the mild areas. Compared to the surface free-air gravity anomalies (Fig 3-1), the airborne gravity anomalies are much smoother.

5.5 Accuracy Assessment

The quality of the airborne gravity data can be evaluated by using three methods:

repeatability analysis, crossover analysis, and comparison with surface gravity data.

Repeatability and crossover analyses are used to evaluate the internal accuracy of the airborne data. External accuracy of airborne gravity can be obtained by comparison with surface gravity.

5.5.1 Repeatability Analysis

Repeatability analysis is a basic method to quantify the accuracy of airborne gravity measurements based on the gravity difference between two repeatable flight lines. Parts of Lines 26 and 55 (Fig 5-3) were flown over twice for repeatability analysis. The repeatability standard deviation is chosen as an index of the measurement precision and is calculated as the standard deviation of the differences at all repeat measurement points. Fig 5-4 shows the difference of the standard deviations of repeat lines 26 and 55 on different filter widths. For both the lines, the standard deviation decreases with increasing filter width but becomes flat beyond a certain filter width. At filter widths smaller than 75 s, the standard deviations of Line 26 are higher than those of Line 55; moreover, beyond 75 s, the standard deviation of Line 55 surpasses that of Line 26. The most important factor that the repeat flights of Line 55 lead to the larger repeatability standard deviation is on the GPS positioning. In the first flight over Line 55, some of the estimated aircraft coordinates appear to be erroneous due to changes in the number of visible GPS satellites and disturbances from unknown sources. The iterative Gaussian filter cannot eliminate these errors. As shown in Fig 5-4, at a filter width of 150 s, the repeatability standard deviation of Line 26 is approximately 3 mgal, and it does not decrease significantly as the filter

width increases. Since an increase in the filter width will eliminate detailed gravity information, it appears that a filter width of 150 s is a trade-off between noise reduction and gravity signal preservation.

5.5.2 Crossover Analysis

Crossover analysis can be used to assess the quality of an airborne gravity survey based on the differences of all the intersecting points. One method of crossover analysis is based on quality weighting assignments using a variance criterion (Mittal, 1984, and Wessel, 1989). However, our approach is to solve the bias and drift at each survey line and corrupt the observed gravity values. The basic expression for crossover analysis is

where g represents the observed gravity value at point r along survey line q and _r^q

q

g is corrupted by a bias, a drift, and a random error. r a and ^q b represent the bias ^q and the drift, respectively, pertaining to survey line q. e is the random error and _r^q

q

t is the time at point r relative to the beginning time of line q. r g is the true gravity _r^q value. At all intersecting points, the observation equations are given by

l

where x^kl_p represents the differenced gravity value at crossover point p pertaining to lines k and l; x^kl_p is the residual; i + m and n are the number of survey lines and crossover points, respectively. Eq (5-16) can be expressed by a matrix representation, such that

AX L

V+ = (5-17)

where V, L, and X are the vectors containing residuals, observations, and parameters (bias and drift), and A is the design matrix. In order to avoid the rank defect, at least one survey line must be fixed. This line should contain the most stable weather

condition of the flight and GPS positioning result.

Fig 5-5 shows the crossover differences (total 736) and a histogram of these values. Most differences are within ±7 mGal, except for a fewer large ones due to inaccurate GPS data and turbulences. If the large differences exceed 15 mgal, they are considered as outliers and not used for the subsequent analyses. The distribution of crossover differences approximately follows a normal distribution, suggesting that these differences are largely due to random noises. Before and after the bias and drift correction on each flight line, the standard deviations of the differences are 4.92 and 2.88 mgal, respectively. The data quality obtained by crossover analysis is consistent with those of other airborne gravity campaigns conducted in other parts of the world.

Fig 5-5 shows the crossover differences (total 736) and a histogram of these values. Most differences are within ±7 mGal, except for a fewer large ones due to inaccurate GPS data and turbulences. If the large differences exceed 15 mgal, they are considered as outliers and not used for the subsequent analyses. The distribution of crossover differences approximately follows a normal distribution, suggesting that these differences are largely due to random noises. Before and after the bias and drift correction on each flight line, the standard deviations of the differences are 4.92 and 2.88 mgal, respectively. The data quality obtained by crossover analysis is consistent with those of other airborne gravity campaigns conducted in other parts of the world.

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