利用空載，地面重力與測高資料計算台灣大地起伏: 研究向上/向下延續與地形效應計算

全文

(1)國. 立. 交. 通. 大. 學. 土木工程研究所博士論文. 利用空載，地面重力與測高資料計算台灣大地起伏: 研究向上/向下延續與地形效應計算 Modeling Taiwan Geoid Using Airborne, Surface Gravity and Altimetry Data: Investigations of Downward/Upward Continuations and Terrain Modeling Techniques. 研究生：蕭宇伸指導教授:黃金維. 中華民國九十六年六月.

(2) 利用空載，地面重力與測高資料計算台灣大地起伏: 研究向上/向下延續與地形效應計算 Modeling Taiwan Geoid Using Airborne, Surface Gravity and Altimetry Data: Investigations of Downward/Upward Continuations and Terrain Modeling Techniques. 研究生：蕭宇伸. Student: Yu-Shen Hsiao. 指導教授：黃金維. Advisor: Cheinway Hwang. 國立交通大學土木工程研究所博士論文. A Thesis Submitted to Department of Civil Engineering College of Engineering National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Doctor in Civil Engineering June 2007 Hsinchu, Taiwan, Republic of China.

(3) 利用空載，地面重力與測高資料計算台灣大地起伏: 研究向上/向下延續與地形效應計算研究生：蕭宇伸. 指導教授：黃金維. 國立交通大學土木工程研究所. 摘要本論文的內容是結合地面、船載、測高與空載重力資料計算台灣與周邊海域的大地起伏模型。空載重力資料是在平均高度5156公尺下利用LaCoste and Romberg (LCR) System II 空載/船載重力儀所測得。為了得到最佳的大地起伏模型，本文研究兩個主要的課題。第一，考慮三種計算剩餘地形效應的方法，這三種方法分別為快速傅立葉轉換、柱狀體法與高斯求積法。在柱狀體法中，將考慮二維地質密度模型的影響。第二，快速傅立葉轉換與最小二乘配置法應用於向下延續的計算。在快速傅立葉轉換計算時，高斯與維納濾波將用平滑向下延續的重力值。最小二乘配置法則分為直接與間接大地起伏計算方法。此外，本文大地起伏計算策略為去除回覆法並搭配最小二乘配置法。空中重力異常與地表重力比較後發現，兩者間較大的差值分部於高山地區，其主要原因為此區域缺乏地面重力資料。在交叉點分析方面，在bias-only改正前後的交叉點差值的均方根分別為4.92 和 2.88 mgal。在重複分析比較方面，150 秒的濾波寬度是平滑空載重力值的最佳濾波寬度。在剩餘地形效應的研究方面，用快速傅立葉轉換計算此效應的大地起伏模型具有最佳的精度。此外，雖然考慮地質密度變化後，大地起伏面會比僅考慮地質密度常數的大地起伏模型有著4公分的變化量，但對改善大地起伏精度卻非常有限。在向下延續分析方面，先把重力向下延續到海水面(包括利用高斯與維納濾波的快速傅立葉轉換與最小二乘配置法)，再計算大地起伏的方法，所表現出的大地起伏模型很相似。然而採用最小二乘配置法直接計算大地起伏所得到的模型與其他方法所計算的比較，在某些區域有著30公分的差值。大致上來說，結合地面與空載重力所計算得到的大地起 I.

(4) 伏，其精度要比僅用地面重力所計算得的要佳，在部分山區可達到10公分以內的精度。. II.

(5) Modeling Taiwan Geoid Using Airborne, Surface Gravity and Altimetry Data: Investigations of Downward/Upward Continuations and Terrain Modeling Techniques Student: Yu-Shen Hsiao. Advisor: Cheinway Hwang. Department of Civil Engineering National Chiao Tung University. Abstract This dissertation is aimed at geoid modeling over Taiwan and the surrounding seas by land-based, shipborne, altimeter, and airborne gravity data. Airborne gravity data was obtained from an airborne gravity survey over Taiwan using a LaCoste and Romberg (LCR) System II air-sea gravimeter at an average altitude of 5156 m. In order to model the best geoid, two main topics are studied. First, three computational methods of the residual terrain model (RTM) effects are considered. The three methods are the fast Fourier transform (FFT), prism, and Gaussian quadrature methods. A 2-D density model of terrain is used in the prism method. Second, the FFT and least squares collocation (LSC) methods are adopted for the computation of the downward continuation (DWC). Both Gaussian and Wiener low-pass filters are used to smooth the downward-continued data by using FFT. Direct and indirect geoid computations are studied in LSC DWC. The methodology of the geoid modeling is mainly based on the remove-compute-restore (RCR) procedure by using LSC. The airborne gravity anomalies are compared with the surface values. Large discrepancies are found to occur over high mountains due to the sparse surface gravity data coverage. The RMS crossover differences before and after a bias-only adjustment are 4.92 and 2.88 mgal. A filter width of 150 s is the optimal width for filtering the airborne gravity data, according to a repeatability analysis. In the investigation of the RTM, the FFT method in the RTM-derived effect computation produces the best geoid accuracy. Although the density variation considered in the geoid modeling yields a 4-cm change in the geoid surface from that using a geological constant, the III.

(6) improvement in the geoid accuracy is extremely small. In the DWC analysis, the methods of DWC to sea level, including FFT with the Gaussian and Wiener filters and LSC, perform similar in geoid modeling. The method of direct geoid determination by LSC provides an obviously different geoid result due to the 30-cm differences of geoid surface from the other geoid models over some areas. Generally, the accuracies of the geoid models from the surface and airborne gravity data outperform the surface-gravity-only geoid models. The improvement in geoid accuracy reaches 10 cm over some high mountainous areas.. IV.

(7) Table of Contents Abstract (in Chinese)………………………………………………………….……Ⅰ Abstract…………………………………………………………….………………..Ⅲ Table of Contents…………………………………………………………….……...Ⅴ List of Tables.…………………...…………………………………………………...Ⅷ List of Figures…………………...………………………………………………..…Ⅸ Chapter 1 Introduction…………………...………………………………………….1 1.1 Background…...…………………………………………………………………...1 1.2 Literature Review…...……………………………………………………………..3 1.3 Outline of Thesis…...……………………………………………………………...6 Chapter 2 Principles of Geoid Determination and Upward/Downward Continuations…...……………………………………..……………………………...9 2.1 Introduction..……………………………………..…………………………...…...9 2.2 Methodologies of Geoid Determination..………………………………………….9 2.2.1 Spherical Harmonic Representation of Gravity Field..……………………..…...9 2.2.2 Stokes Integration………………………………………………………………11 2.2.3 Least Squares Collocation……………………………………………………...12 2.3 Remove-Compute-Restore Procedure……………………………………………15 2.3.1 Long-Wavelength Reference Geopotential Model……………………………..16 2.3.2 Residual Terrain Model………………………………………………...............17 2.4 Quasi-Geoid Correction……………………………………………….................18 2.5 Upward and Downward Continuations…………………………………………..19 2.5.1 Continuation by Fast Fourier Transform……………………………………….20 2.5.2 Continuation by Least Squares Collocation……………………………………23 Chapter 3 Data for Geoid Modeling…...…………………………………………..24 3.1 Introduction…………………………………………………………………..…..24 3.2 Surface Gravity…………………………………………………………………..24 3.2.1 Land Gravity…………………………………………………………….……..24 3.2.2 Shipborne Gravity……………………………………………………………...24 3.3 Altimeter-Derived Gravity……………………………………………………….25 3.4 Geopotential Model………………………………………………………………25 3.5 Digital Elevation Model………………………………………………………….26 V.

(8) 3.6 Density Model……………………………………………………………………26 3.7 GPS/Leveling Points for Evaluation……………………………………………..27 Chapter 4 RTM Effects in Geoid Modeling: Comparison of Three Methods…..33 4.1 Introduction……………………………………………………………………....33 4.2 RTM Effects by FFT……………………………………………………………..33 4.3 RTM Effects by Prism……………………………………………………………36 4.4 RTM Effects by Gaussian Quadrature……………………………………………39 4.5 Design of Experiments…………………………………………………………...41 4.6 Results……………………………………………………………………………44 4.6.1 Results from RTM-Derived Gravity and Geoid……………………………..…44 4.6.2 Results of Geoid Modeling…………………………………………………….45 Chapter 5 Airborne Gravity Data of Taiwan……………………………………...55 5.1 Introduction………………………………………………………………….…...55 5.2 Data Reduction in Airborne Gravity………………………………………...…...55 5.2.1 Gravity Reduction……………………………………….…..............................55 5.2.2 Aircraft Positioning……………………………….…........................................58 5.3 A Taiwan Airborne Gravity Survey……………………………….…..................59 5.3.1 Survey Campaign…………………………….…...............................................59 5.3.2 Data Processing…………………….…..............................................................62 5.4 Results of the Airborne Gravity Survey.................................................................64 5.5 Accuracy Assessment.............................................................................................64 5.5.1 Repeatability Analysis.........................................................................................64 5.5.2 Crossover Analysis..............................................................................................65 5.5.3 Comparison with Surface Gravity Data..............................................................66 Chapter 6 Geoid Modeling Using Combined Airborne and Surface Gravity Data…………………………………………………………………………………..72 6.1 Introduction………………………………………………………………………72 6.2 Continuation to Sea Level and Merging with Surface Gravity…………………..72 6.2.1 DWC by FFT with Gaussian Filter…………………………………………….72 6.2.2 DWC by FFT with Wiener Filter………………………………………………73 6.2.3 DWC by LSC…………………………………………………………………..74 6.3 Direct Use for Geoid Modeling…………………………………………………..74 6.4 Design of Experiments…………………………………………………………...75 6.5 Evaluating Downward-Continued Airborne Gravity Data……….........................80 VI.

(9) 6.5.1 Results of DWC by FFT………..........................................................................80 6.5.2 Results of DWC by LSC……….........................................................................81 6.6 Results of Geoid Modeling....................................................................................82 Chapter 7 Summary, Conclusions, and Recommendations.................................100 7.1 Summary..............................................................................................................100 7.2 Conclusions..........................................................................................................101 7.3 Recommendations for Future Work.....................................................................103 Reference...................................................................................................................105 Curriculum Vitae......................................................................................................112. VII.

(10) List of Tables 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………………42 Table 4-2 Statistics for RTM-derived gravity anomalies (mgal) ……………………44 Table 4-3 Statistics for RTM-derived geoids (m) …………………………………...45 Table 4-4 Statistics regarding the differences (m) between the observed and modeled geoidal heights at four leveling routes……………………….....................................47 Table 5-1 Overview of the L & R Air-Sea Gravity System II…………….................60 Table 6-1 Statistics of the differences between surface gravity anomalies and DWC (FFT) gravity anomalies…………...............................................................................81 Table 6-2 Statistics of the differences between surface gravity anomalies and DWC (LSC) gravity anomalies…………..............................................................................82 Table 6-3 Statistics of differences (in meter) between the observed and modeled geoidal heights at four leveling routes………….........................................................84 Table 6-4 The original average geoid errors (in meter) in the east route………….....85. VIII.

(11) List of Figures Fig. 1-1 Terrain and bathymetry around Taiwan………………………………………8 Fig. 2-1 Potential at point P due to the earth mass…………………………………...10 Fig. 2-2 Three different wavelengths of geoid undulation…………………………...16 Fig. 2-3 Residual terrain model (RTM).……………………………………………...17 Fig. 2-4 Physical surface of the earth………………………………………...............18 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………………28 Fig. 3-2 (a) Gravity anomalies and (b) geoid heights globally and over Taiwan obtained from the EIGEN-GL04C coefficients………………………………………29 Fig. 3-3 DEMs used in the geoid modeling………………………………………….30 Fig. 3-4 Density model over Taiwan…………………………………………………31 Fig. 3-5 Four leveling routes for evaluating the geoid accuracy……………………..32 Fig. 4-1 Geometry of the RTM-derived effects in the prism method………………..38 Fig. 4-2 Computational inner and outer zones at point P…………………………....38 Fig. 4-3 Geometry depicting an RTM-derived effect in the Gaussian quadrature…...39 Fig. 4-4 Flowchart for the geoid modeling procedure……………….........................43 Fig. 4-5 RTM-derived gravity anomalies. (a) case 1, (b) case 2, (c) case3, and (d) case 4……………................................................................................................................48 Fig. 4-6 RTM-derived geoids. (a) case 1, (b) case 2, (c) case3, and (d) case 4……...49 Fig. 4-7 Covariances for (a) surface gravity-surface gravity covariance matrix, (b) geoid-surface gravity covariance matrix, and (c) geoid-geoid covariance matrix…...50 Fig. 4-8 Residual gravity anomalies. (a) case 1, (b) case 2, (c) case3, and (d) case 4…................................................................................................................................51 Fig. 4-9 Residual geoids. (a) case 1, (b) case 2, (c) case3, and (d) case 4...................52 Fig. 4-10 Geoid models obtained in (a) case 1, (b) case 2, (c) case 3, and (d) case 4....................................................................................................................................53 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.................................................................54 Fig. 5-1 (a) Airborne gravity survey lines and GPS tracking stations for precise aircraft positioning. (b) The L&R Air-Sea Gravity System II gravimeter. (c) The King-Air Beechcraft-200 aircraft. (d) Inside of the King-Air Beechcraft-200............61. IX.

(12) Fig. 5-2 Flow chart of the airborne gravity data process implemented by NCTU.......63 Fig. 5-3 Gravity anomalies at the average flight altitude of 5156 m...........................68 Fig. 5-4 Standard deviation of the differences in the gravity anomalies obtained from two repeat flights………………………………………………………......................68 Fig. 5-5 Distribution and histogram of the crossover differences of gravity anomalies……………………………………………………………………………..69 Fig. 5-6 Zero-padding (100%) of the used gravity field G ( f x , f y ) in UWC and DWC………………………………………………………………………………….70 Fig. 5-7 The 1D (left) and 2D transfer functions of UWC…………………………...70 Fig. 5-8 (a) Bouguer anomalies at an altitude of 5156 m, (b) Bouguer anomalies of surface gravity data on the grid, (c) Surface-upward-continued Bouguer anomalies, (d) differences between surface and airborne Bouguer anomalies………………………71 Fig. 6-1 (a) All the input gravity data for geoid computation. (b) Zoomed-in view of the black rectangular area in (a)…..……………………………………………….....77 Fig. 6-2 Covariances: (a) Surface gravity–airborne gravity covariance matrix. (b) Airborne gravity–airborne gravity covariance matrix. (c) Airborne gravity–geoid covariance matrix…..………………………………………………...........................78 Fig. 6-3 Flow chart of geoid modeling.……………………………………………...79 Fig. 6-4 Transfer functions for 1-D (left-hand-side) and 2-D cases of DWC……......85 Fig. 6-5 1D frequency response of Gaussian filter………………………..................86 Fig. 6-6 2D frequency response of Gaussian filter……………………......................87 Fig. 6-7 Downward-continued Bouguer anomalies by Gaussian filter……................88 Fig. 6-8 Differences between the Bouguer anomalies of the surface and downward-continued airborne data by Gaussian filter…….........................................89 Fig. 6-9 1D frequency response of Wiener filter……..................................................90 Fig. 6-10 2D frequency response of Wiener filter……................................................91 Fig. 6-11 Downward-continued Bouguer anomalies by Wiener filter…….................92 Fig.. 6-12. Differences. between. the. Bouguer. anomalies. of. surface. and. downward-continued airborne data by Wiener filter……............................................93 Fig. 6-13 (a)Free-air anomalies of airborne data at 5156 m. (b) RTM-derived effects of airborne data at 5156 m. (c)Residual gravity anomalies of airborne data at 5156 m. (d) Downward-continued residual gravity anomalies by LSC method…....................94 Fig. 6-14 Differences between the residual gravity anomalies of surface and. X.

(13) downward-continued data by LSC method…..............................................................95 Fig. 6-15 The geoid models. (a) case A. (b) case B. (c) cases C. (d) case D (e) case E…...............................................................................................................................96 Fig. 6-16 Differences between all the geoid models. (a) Between cases A and B. (b) Between cases A and C. (c) Between cases A and D. (d) Between cases A and E…...97 Fig. 6-17 Differences (in m) between the observed and modeled geoidal heights (cases A~E) along four leveling routes…....................................................................98 Fig. 6-18 The adjusted geoid errors for (a) case A. (b) case B. (c) cases C. (d) case D (e) case E…..................................................................................................................99. XI.

(14) Chapter 1 Introduction 1.1 Background Taiwan experiences a large amount of seismic activity because it is located over the junction of the Eurasia plate and the Philippine Sea plate. The uplift and subsidence of land are created due to the collision of these two plates. Most areas of Taiwan Island are subjected to northwest-southeast compression at an average rate of 8.2 cm/year (Yu et al., 1997). The subduction of the Philippine Sea plate into the Eurasia plate creates a deep trench and large negative gravity anomalies to the east of Taiwan. On the other hand, it also creates the Central Range with a high terrain and huge positive gravity anomalies on land. The maximum altitude is at the Central Range, reaching 3952 m, which corresponds to the highest peak in East Asia. In eastern Taiwan and the surrounding sea, the mountains and the seabed, which are only several km away from the coast, reach heights of approximately 2000 m and –5000 m, respectively. Due to the extremely rough terrain and bathymetry (Fig 1-1), geoid modeling over Taiwan Island and its surrounding marine areas is quite a challenge for geodesists and geophysicists. Geoid determination with high accuracy is a primary goal for geoscientists. The importance of the geoid for geodesists is that it is a reference surface for orthometric heights. Once a reference surface is established, orthometric heights referred to the local vertical datum are obtained. In addition, it is feasible to determine the orthometric heights by using GPS. A high accuracy geoid is the key factor for obtaining orthometric heights without leveling. If we have a high quality geoid model, orthometric heights can be efficiently and economically computed using GPS-derived ellipsoidal heights. For oceanographers, the geoid is useful for the investigation of currents, tides, and sea surface topography. For geophysicists, the geoid can be used to understand the characteristics of the earth’s interior sources. Besides geodetic purposes, the geoid is also applied in mapping, photogrammetry, and remote sensing. This is why most countries around the world are making efforts to compute their own precise geoid models. The estimation of the topographic effect is important for geoid determination, especially over mountainous areas. This estimation can be used to calculate the effects 1.

(15) of high frequency on gravity and the geoid; these cannot be calculated using the geopotential model and local gravity data. In most investigations related to the topographic effect, it is assumed that the density of the topographic mass is constant. However, several studies in recent years have taken into consideration the influence of the density variation of the topographic mass. In addition, the chosen method and digital elevation models (DEM) used in topographic effect computations also need to be focused upon in order to obtain a more precise result in an efficient manner. Airborne gravimetry is a method to determine the gravity field by measurements from an aircraft. Based on this method’s feature, airborne gravity data are valuable for areas with sparse gravity data, such as high mountains where data are always collected along the roads in valleys. Airborne gravity data are also useful for coastal regions wherein the gravity data coverage, especially over shallow water areas, obtained from satellite altimetry and land gravimetry data is of poor quality. Therefore, airborne gravimetry is suitable for Taiwan Island where over 75% of the terrain comprises hills and high mountains and 70% of the coast is near shallow water areas. Poor gravity data coverage results in poor accuracy of geoid modeling. Airborne gravity surveys with equally spaced tracks can make data coverage denser and bring improvements in the geoid computation. Another interesting topic in the recent years has been how to combine different kinds of gravity data to compute a precise geoid. These data include terrestrial, shipborne, airborne, and altimeter-derived gravity. The combination of different types of gravity measurements is a challenging task for geodesists due to their different resolutions and characteristics. Airborne gravity data have an unusual property in that the gravity field level is different from that of other types of data. Thus, the technique of downward continuation (DWC) is important for airborne data to press aerial gravity field to the level which we are interested in. The objective of this thesis is to determine the most accurate geoid over Taiwan and its surrounding sea area by the use of surface, altimeter-derived and airborne gravity. Based on this objective, there are several main issues to be investigated in this dissertation: (1) Which is the best method for topographic effect computation? (2) Is the consideration of the density variation of the topographic mass necessary in topographic effect computation? Can it be ignored? (3) What is the quality of the airborne gravity data used in the geoid modeling over Taiwan? (4) What is the best DWC method that can be applied to airborne data? (5) What is the ideal geoid model 2.

(16) by combining all types of gravity data? All these topics are important and have been investigated in this geoid modeling study. In brief, this dissertation focuses on how to obtain the most accurate geoid over Taiwan by selecting the best (1) topographic effect computation method, (2) DWC technique, and (3) geoid determination method.. 1.2 Literature Review Geoid determination has been of interest to geodesists for more than a century. Basically,. two. types. of. methods. are. usually. used. for. local. geoid. determination—Stokes integration and least-squares collocation (LSC). They are deterministic and stochastical methods, respectively. Stokes integration can be performed very quickly using the fast Fourier transform (FFT) on gridded data. On the contrary, LSC requires a larger computational effort. Stokes integration generally only uses one data type with uniform noise. However, LSC can accept hybrid data with individual noises. Stokes integration is usually used for the continental areas and geoid models over several regions around the world (e.g., Boziane, 1996; Denker et al., 1997; Forsberg et al., 1996; Smith and Milbert, 1999; Sideris, 1995). The application of LSC in physical geodesy has been discussed in detail by Moritz (1980). The first centimeter geoid was computed for an area around Hannover in Germany in 1987 (Denker and Wenzel, 1987) by LSC. This method was subsequently used in many countries and was met with success (e.g., Sevilla, 1997; Hwang, 1997; Tscherning et al., 2001). Compared to Stokes integration, LSC gives error estimates and error covariances that reflect the data distribution and quality. For modeling the local gravity and geoid field at present, the LSC method has been proven to be a powerful technique. Geoid modeling using Remove-Compute-Restore (RCR) procedure by Stokes integration over Taiwan was first investigated by Tsuei (1995). In subsequent years, a number of studies based on RCR procedure by LSC were carried out, e.g. Hwang (1997, 2001, 2003, 2005), Hwang et al. (2006a, 2006b, 2007b). Most of these results show a geoid accuracy of several centimeters over the west plain but of 1~2 decimeters over high mountains. The topographic effect in geoid determination has also been studied for many years, especially in rough terrains. The relevant investigations can be divided into two main parts: (1) the effects of the residual terrain model (RTM) and (2) Helmert’s. 3.

(17) second method of condensation. There are several methodologies to determine the topographic effect. An earlier research containing the complete computation of this effect can be found in Forsberg (1984). It presented the FFT and prism methods to calculate the RTM-derived effect. Sjoberg (2000) used Helmert’s second method of condensation to reduce the topography. Omang and Forsberg (2000) investigated three different methods of dealing with topography in geoid modeling: the isostatic, Helmert condensation, and RTM methods. Other studies about the topographic effect in geoid modeling include Forsberg (1985), Nahavandchi and Sjoberg (2001), and Flury (2006). Furthermore, we usually assume that the density of the topographic mass is constant (2.67 g/cm3) while computing the topographic effect, but recently, several investigations have been performed to study the impact of more realistic density variations of the topographic masses. Martinec (1998) showed that the geoid can be changed approximately to the decimeter level by considering the lateral density variation of the topographical masses. Pagiatakis et al. (1998) reported that the effect of lateral density variations on the geoid can reach nearly 10 cm in the Skeena region in Canada, where the terrain is hilly. Huang (2002) showed that the total density variation effect on the geoid heights ranges from –7.0 cm to 2.8 cm in the Canadian Rocky Mountains. It is evident that the use of the digital topographical density model will significantly improve the accuracy of the geoid. Other studies about density variations can be found in Huang et al. (2001), Hunegnaw (2001), Smith (2002), Kuhn (2003), and Sjoberg (2004). Airborne gravity surveys have been performed for over forty years, but geodesists have been paying more attention to them recently due to the advancement of the methodology, improvement in instrumentation, and development of the precise kinematic GPS in the past decade. Due to the recent and rapid development of these techniques, 1~2 mgal and half wavelengths of 3–4 km can be achieved by airborne gravity surveys (Schwarz and Li, 1997). The first test of the airborne gravity survey was made by Thompson and LaCoste (1960). The main objective of this flight is to show that gravity measurement from a flying aircraft is feasible. The first large-scale airborne gravity experiment was performed over Greenland (Brozena, 1992). In addition, several airborne gravity surveys have also been performed in places whose terrains are similar to that of Taiwan, such as the Rocky Mountains, the Alps, and Malaysia. Airborne gravimetry tests are often conducted in the Rocky Mountains due to the 4.

(18) complex topography. In 1995, an airborne gravity survey (Wei and Schwarz, 1998) was carried out over the Rocky Mountains. The gravity system includes an inertial navigation system (INS) and two GPS receivers on the aircraft. The survey lines contain four flights with the same trajectory, which has an east-west profile of 250 km. The flying altitude and speed were 5.5 km and 430 km/h, respectively. The gravity result shows that the repeatability standard deviation is about 2 mgal with a filter length of 120 s and about 3 mgal with a filter length of 90 s. The standard deviation of the difference between the airborne gravity and upward continued ground gravity is about 3 mgal for both filter lengths. In the next year, another airborne gravity survey was conducted in the Rocky Mountains again (Glennie and Schwarz, 1999). The mission was carried out over a single 100 × 100 km2 area with a line spacing of 10 km. The analyses of the crossover differences showed a root mean square (RMS) agreement at the level of 1.6 mgal. In 1998, an airborne gravity survey was carried out over the Alps (Verdun et al., 2003). The mission consisted of 18 NS and 16 EW survey lines with a line spacing of 10 and 20 km, respectively. The gravimeter, which is a LaCoste & Romberg relative air/sea gravimeter (type SA), was mounted in a DeHavilland Twin Otter aircraft flying at a constant altitude of 5100 m and a mean ground speed of about 280 km/hr. Seven ground based GPS reference stations were used to determine the positions of the aircraft. The accuracies of the Bouguer anomaly are determined from the crossover analysis (15.34 mgal RMS before adjustment and line selection) and the ground upward continuation (UWC) (7.68 mgal RMS for a spatial resolution of 8 km). Airborne gravimetry in Malaysia was carried out by National Land Survey and Cadastre (KMS) of Denmark. The airborne gravity survey over the entire peninsula and Brunei was conducted with a 5 km line spacing, using an An-38 aircraft. More than 600 hours were flown with LaCoste and Romberg gravimeters (models S-93 and S-99) to collect the basic airborne gravity data at a flight speed of 150–250 km/hr and an aircraft altitude of less than 4500 m, which typically corresponds to a height of 300–1000 m above the topography, depending on the weather conditions. The resulting gravity anomalies of the crossover analysis were higher than 2 mgal RMS. Airborne gravity surveys are still conducted in many countries and regions. These areas include North Carolina (Brozena and Peters, 1988), Skagerrak (Kearsley et al., 1998), Azores Islands (Bastos et al., 1998), Antarctica (Bell et al., 1999), the 5.

(19) Nordic/Baltic area (Forsberg and Solheim, 2000), Greenland/Svalbard (Forsberg et al., 2003a), the Arctic sea (Childers et al., 2001 and Forsberg et al., 2003b), Lincoln Sea and Wandel Sea (Olesen et al., 2003), Baltic Sea, the Great Barrier Reef, Crete Island, and Mongolia. These missions yielded an average RMS error of 2 mgal based on the crossover comparisons and an average interior geoid accuracy of 5 cm based on reliable GPS/leveling data. The precision of geoid modeling has improved in recent years due to the development of airborne gravimetry. Forsberg et al. (2000) showed that the routine accuracy of airborne gravimetry is at the 2 mgal level, which may translate into a geoid accuracy of 5–10 cm on a regional scale. Kearsley et al. (1998) indicated that the gravity field determined from an aircraft with flight separations of 5 to 10 km can be used to evaluate precise (2 cm) relative geoid heights over north Jutland. Schwarz and Li (1996) pointed out that a centimeter geoid can be obtained if the minimum wavelength resolved is about 14 km in flat areas and 5 km in mountainous areas. In the airborne gravimetry in Malaysia as mentioned above, the data contributed to a geoid accuracy of smaller than 5 cm. Combining the different types of gravity data for geoid determination is also an interesting topic. Novak et al. (2003) reported that the first geoid model computed using the combination of airborne and global gravity data had a difference standard deviation of 5.5 cm; this is comparable to the reference geoid computed only from the ground gravity data. The second geoid model, based on the combination of the airborne and ground gravity data, had a difference standard deviation of 4.7 cm by comparison of the same reference geoid. Jekeli and Kwon (2002), and Serpas and Jekeli (2005), used the horizontal components of airborne gravity observations and also reported a sub-decimeter precision in the determination of the relative local geoid. Other research about geoid determination using airborne gravity data include, among others, Kern et al. (2003), Novak (2003), Olesen (2003), Li (2000), Serpas (2003), Bayoud and Sideris (2003), and Olesen et al. (2002).. 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. 6.

(20) 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.. 7.

(21) 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.. 8.

(22) 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). 9.

(23) b = G ∫∫∫ M. r′ − r r′ − r. 3. dm. (2-1). where r ′ and r are the geocentric position vectors of the element mass dm and the attracted point P.. P. Earth surface r. dm r′. Earth center Fig. 2-1 Potential at point P due to the earth mass. The corresponding potential V and the earth gravitation b have the relationships. b = grad V. (2-2). Thus, the gravitational potential of the earth can be given by. Vearth = G ∫∫∫ M. 1 ρ dm = G ∫∫∫ dv r′ − r r′ − r M. (2-3). where ρ and dv are the earth’s density and volume element, respectively.. lim V = 0 . V is harmonic outside the spheroid and can be determined by using a r →∞. spherical harmonic function given by (Torge, 1991). ∞ GM ⎛a⎞ (1 + ∑ ⎜ ⎟ V = r l =2 ⎝ r ⎠. n n. ∑. (C nm cos mλ + S nm sin mλ ) P nm (cos ϑ )). m =0. 10. (2-4).

(24) 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 C nm and S nm are fully normalized spherical harmonic coefficients, which are mass integrals that represent the mass distribution within the central body. P nm 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). ∆g long =. GM R2. ). (2-5). cos mλ + S nm sin mλ P nm (sin φ ). (2-6). N. ∑ (n − 1) n=2. ∑ (C n. nm. cos mλ + S nm sin mλ P nm (sin φ ). m =0. and N. N long = R∑ n=2. ∑ (C n. nm. ). m =0. 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. T (φ , λ ) =. R 4π. ∫∫σ S (ψ )∆gdσ. (2-7). where ∆g is gravity anomaly. σ represents the unit sphere and dσ denotes the element of solid angle. S (ψ ) is Stokes’ kernel and is expressed as. S (ψ ) =. 1 sin(ψ 2). − 6 sin. ψ 2. + 1 − 5 cosψ − 3 cosψ ln(sin. 11. ψ 2. + sin 2. ψ 2. ). (2-8).

(25) Using Bruns’ formula, we can obtain the geoid undulation as. N=. R 4πγ. ∫∫σ S (ψ )∆gdσ. (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). s = C sl (C ll + D) −1 l. (2-10). where s and l are sets of signals and observations, respectively. C ll is the covariance matrix of l and C sl is the covariance matrix between s and l. D is the matrix of the 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) E ss = C ss − C sl C ll−1C ls. (2-11). where Ess denotes the error covariance matrix. In the case of geoid determination by LSC, the formulae of signal and error are. (. N = (C Ng ) C gg + Dg. ) (∆g ) −1. (2-12). and 12.

(26) ⎞ ⎛ 1 ⎛ ⎞ δn 2 = s⎜ C NN − (C Ng )⎜ C gg + D g ⎟ (C Ng )T ⎟ ⎟ ⎜ −1. ⎝. ⎝. s. ⎠. ⎠. (2-13). where vector ∆g contains gravity anomalies. C Ng and C gg are covariance matrices for geoid-gravity and gravity-gravity, respectively. D g is a diagonal matrix containing error variances of gravity anomalies. In (2-13), s is a scale factor, which can be determined by i =m ⎛ i =m ⎜ ∑ (∆g i − ∆g ) 2 ∑ δg i2 − i =1 s = ⎜ i =1 ⎜ m m ⎜ ⎝. ⎞ ⎟ ⎟ ⎟ ⎟ ⎠. c gg. (2-14). where ∆g is the average gravity anomaly, δg 2 is error variance of gravity anomalies, c gg is error variance of the gravity anomalies derived from a geopotential 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. σ n2 (∆g, ∆g ) =. A(n −1 ) s n+2 (n − 2 )(n + B). (2-15). s can be expressed as. 13.

(27) ⎛R ⎞ s=⎜ B ⎟ ⎝ rr ′ ⎠. 2. (2-16). 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 n > 360 . 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 ∞. 360. C ∆g∆g ( P, Q) = ∑ σ~n (∆g , ∆g ) s n + 2 Pn (cos ΨPQ ) + n=2. n = 361. ∞. 360. CTT ( P, Q ) = ∑ σ~n (T p , TQ ) s n +1 Pn (cos ΨPQ ) + R 2 n=2. A(n-1 ). ∑ (n-2 )(n + B) s A. n+ 2. Pn (cos ΨPQ ). ∑ (n-1 )(n-2 )(n + B) s. n = 361. n +1. (2-17). Pn (cos ΨPQ ) (2-18). 360 ⎛ R2 ⎞ ∞ A ⎟⎟ ∑ s n +1 Pn (cos ΨPQ ) (2-19) C ∆gT ( P, Q) = ∑ σ~n (T p , ∆g Q ) s n +1 Pn (cos ΨPQ ) + ⎜⎜ r (n2 )(n B) + n=2 ⎠n =361 ⎝. C nn ( P, Q) =. 1. γ2. CTT ( P, Q). (2-20). where σ~n (∆g , ∆g ) , σ~n (Tp , TQ ) , and σ~n (T p , ∆g Q ) 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. Pn is the Legendre polynomial of degree n and ΨPQ is the spherical distance between P and Q. γ 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). 14.

(28) 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 ∆g can be expressed as. ζ = ζ long + ζ res + ζ short. (2-21). and ∆g = ∆g long + ∆g res + ∆g short. (2-22). where ζ long and ∆g long are the long-wavelength height anomaly and gravity anomaly, respectively; ζ short and ∆g short , the short-wavelength height anomaly and gravity anomaly, respectively; and ζ res and ∆g res , the residual height anomaly and 15.

(29) 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.. N long. N long + N res. N long + N res + N short. Ellipsoid Fig. 2-2 Three different wavelengths of geoid undulation. N long , N res , and N short. 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,. 16.

(30) 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.. P. True elevation surface RTM. Mean elevation surface 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. V = G ∫∫∫ RTM. dm ( x − x p ) + ( y − y p ) 2 + ( z − s) 2 2. (2-23). where dm is a mass element of RTM and ( x p , y p , z p ) and (x, y, z) are the 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 17.

(31) 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.. Sea level. HN. Topography HO. h Geoid. Quasi-geoid. ζ. N. Ellipsoid Fig. 2-4 Physical surface of the earth. h, H O , H N , N, and ζ 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). ζ −N ≈−. ∆g B. γ. H. (2-24). 18.

(32) where ∆g B is the Bouguer anomaly, γ is normal gravity, and H is the topographic height. Eq (2-24) can also be written as. N ≈ζ −. 2πGρ. γ. H2. (2-25). where ρ is the density of the terrain mass and G is the gravitational constant.. 2πGρ 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:. d = A(m). (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. 19.

(33) (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 ∇g = 0. (2-27). The gravity potential at some level z = 0 is assumed as U ( x, y, z ) = U ( x, y,0) = f ( x, y ). (2-28). where f ( x, y ) is some known function. If the problem is to calculate the potential from z = 0 to any other level z = h, it is called an UWC of the gravity potential. In contrast to UWC, if the problem is to compute the potential from z = h to z = 0, it is 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. ⎧ U ( x, y,0) = A[U ( x, y, h)] upward continuation ⎪ ⎨ ⎪⎩U ( x, y, h) = A[U ( x, y,0)] downward continuation. (2-29). 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. 20.

(34) implemented by using FFT in the frequency domain. Let the vertical component of the gravity field in the z = 0 and z = h planes be g ( x, y, z ) | z =0 = g ( x, y, z = 0). (2-30). and g ( x, y , z ) | z = h = g ( x, y , z = h ). (2-31). 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 = h) can be written as (Buttkus, 2000). g ( x, y , z = h ) =. h 2π. ∞. g (α , β ,0). ∞. ∫ ∫ [h −∞ −∞. 2. + (x − α )2 + ( y − β )2. ]. 3. dαdβ. (2-32). 2. We can use a convolution integral to represent Eq (2-32) as g ( x, y, z = h) = wupward ( x, y ) | hz=0 *g (α , β ,0). (2-33). where wupward ( x, y ) | hz=0 =. h 1 ( 2 ) 3 2 2π (h + x + y 2 ) 2. (2-34). wupward ( x, y ) |hz=0 is the impulse response function for UWC from the z = 0 plane to the. z = h plane. On the other hand, the two-dimensional Fourier transform is given by. W ( fx, fy ) = ∫. ∞. ∫. ∞. −∞ −∞. w( x, y )e. − i 2π ( f x x + f y y ). dxdy. (2-35). where w( x, y ) is a nonperiodic function of real variables x and y. W ( f x , f y ) represents w( x, y ) in the two-dimensional wavenumber domain. f x and f y denote the numbers of cycles per unit distance. If w( x, y ) is substituted in Eq (2-34), the corresponding wavenumber response function becomes 21.

(35) − i 2π ( f x + f y ). x y he dxdy 3 − ∞ ∫− ∞ 2π (h 2 + x 2 + y 2 ) 2. WUWC ( f x , f y ) | zz ==0h = ∫ =e. where f r =. ∞. ∞. − 2πh. f x2 + f y2. (2-36). f x2 + f y2 . Therefore, the UWC from the z = 0 plane to the z = h plane. can be expressed in the wavenumber domain as G ( f x , f y ) | z = h = e −2πhf r G ( f x , f y ) | z =0. (2-37). In contrast to UWC, the wavenumber response function of DWC from the z = h plane to the z = 0 plane is given by. G ( f x , f y ) | z =0 =. G( f x , f y ) | z =h e. − 2πhf r. = e 2πhf r G ( f x , f y ) | z = h. (2-38). 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 G ( f x , f y ) | z =0 = e 2πhf r G ( f x , f y ) | z = h S ( f x , f y ). (2-39). where S ( f x , f y ) is a low-pass filter in the wavenumber domain. If. fr. approximates to infinity, G ( f x , f y ) | z =0 approximates to zero such that. lim = e 2πhf r G ( f x , f y ) | z = h S ( f x , f y ) = 0. (2-40). f r →∞. In this study, UWC by FFT will be used to compare airborne and surface gravity 22.

(36) 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. (. ∆g h2 = C g h1 g h2. ) (C. g h1. + Dg h1. ) (∆g ) −1. (2-41). h1. where C g h1 g h2 is the covariance matrix for gravity at level h1 and level h2 and C g h is the covariance matrix for gravity at level h1 . Dg h is the variance of noise of 1. 1. the gravity data obtained at level h1 . In this study, C g h1 g h2 and C g h1 are both 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.. 23.

(37) 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 this study, the data was only considered for the locations between 119.2–122.8 E and 24.

(38) 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 25.

(39) 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 90′′ × 90′′ 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 26.