台灣大地起伏模型之精進及其在DEM製作之應用
50
0
0
全文
(2) 台灣大地起伏模型之精進及其在 DEM 製 作之應用 An Improved Taiwan Geoid Model And Its Application To DEM Generation. 研 究 生：李佩珊. Student：Pei-Shan , Lee. 指導教授：黃金維. Advisor：Dr. Cheinway Hwang. 國立交通大學 土木工程學系 碩士論文. A Thesis Submitted to Department of Civil Engineering College of Engineering National Chiao Tung University in Partial Fulfillment to the Requirements for the Degree of Master In Civil Engineering June 2005. Hsinchu, Taiwan, Republic of China.
(3) 台灣大地起伏模型之精進及其在 DEM 製 作之應用 學生：李佩珊. 指導教授：黃金維. 博士. 國立交通大學土木工程學系. 摘要. 台灣重力大地起伏模型採用去除計算回復法，其將大地起伏分成長、中及短 波長三部分，此法使用全球大地位模式及台灣數值高程模型，殘餘大地起伏則以 最小二乘配置法計算。新增加一改正項大幅提昇大地起伏在高山區精度。長波長 部份採用最新地球重力模型 GGM01、GGM02、EIGEN-CG01C，經測試最佳成果為 GGM02C 結合 EGM96 計算之台灣重力大地起伏模型，在四條一等水準線用實際觀測值導得 之大地起伏檢核本文大地起伏模型，其差值的標準偏差均在 10 公分以內。採用 1920 個 GPS/水準點上之大地起伏模型與實際觀測大地起伏差值，製作一諧和面修正重 力大地起伏模型，經在四條水準線測試成果為標準偏差 2、3 公分左右。用此改善 後之大地起伏模型應用在航測及光達製作 DEM 時之高程轉換上，首先在新竹測試 區採用航測法生產橢球高再以大地起伏模型改算為正高，與以航測正高高程控制 之正高成果相比，因地面控制點之選擇不同造成兩法高程差值的標準偏差有 11~40 公分不等之成果；在屏東測試區比較光達橢球高經大地起伏模型改算之正高 DEM， 與水準正高之差值獲得標準偏差 15.8 公分。. i.
(4) An Improved Taiwan Geoid Model And Its Application To DEM Generation Student：Pei-Shan , Lee. Advisor：Dr. Cheinway Hwang. Institute of Civil Engineering National Chiao Tung University. Abstract A. gravimetric. geoid. model. of. Taiwan. is. computed. using. remove-compute-restore technique, which divides a geoid model into the long, medium and short wavelength parts. This technique requires a global geopotential model and an elevation model. The residual geoid modeling is done by least-squares collocation. An correction term. accounting for the. terrain effect is introduced, and it improves greatly the accuracy of modeled geoidal undulations, especially in the high mountain areas. For the long wavelength part, this research experiments with earth gravity models GGM01, GGM02 and EIGEN-CG01C. The combined GGM02C-EGM96 model yields the best result in geoid modeling. The modeled and observed geoidal undulations are compared at four first-order leveling lines, resulting in standard deviations of geoidal differences of less than 10 cm. A refined. ii.
(5) geoid model is created by applying a harmonic surface of geoidal differences between modeled and observed geoidal undulations at 1920 first-order benchmarks, where GPS-derived ellipsoidal heights are available. The refined geoid model yields standard deviations of geoidal differences of 2 to 3 cm at the same four leveling lines. The refined geoid model is used in the photogrammetry and Lidar methods to generate orthometric heights. In a testing area in Hsinchu County, ellipsoidal heights are first generated by the photogrammetry method, orthometric heights are then computed using the geoid model.. Comparing the orthometric heights generated by this. approach. the. and. conventional. orthometric-height-controlled. photogrammetry approach, the standard deviations of height differences range from 11 to 40 cm, depending on the scenarios of the ground control points for the aerial triangulation adjustment. At a testing line in Pintung County,. the. standard. deviation. of. the. differences. between. Lidar-geoid-derived orthometric height and levelling-derived orthometric heights is 15.8 cm.. iii.
(6) 誌謝. 在交大求學的這兩年我學到許多知識與經驗，也認識很多不同的人，更受到 了大家的照顧與幫助，感謝所有人給予我美好的回憶。. 非常感謝恩師黃金維教授在學業以及生活上的教導與照顧，尤其在研究過程 中觀念的啟發、研究方法的指導、做學問的態度都是學生們非常好的表率，在此 僅向恩師致上最誠摯的敬意與感激。感謝史天元教授、陳春盛教授與李振燾教授 於課業上的教誨。感謝口試委員楊名教授與崔國強教授的熱心指正與寶貴意見， 使本論文更臻完善。. 感謝博士班蕭宇伸學長在研究過程的大力協助，提供許多寶貴經驗、資料及 指點。感謝博士班學長姊宣昶、成機、廷融、榮寬、欣瑩、弘基、自強、俊廷、 鉅富、大綱、豫麒、志敏以及已畢業的學長姐米粉、宜珊、大雄、勇者、虎爛、 福利、阿福與介嵐哥在學業上的照顧。感謝學長 BOSS、小支及同班的好同學貓哥、 祐廷、展鵬、粘、達哥、嘉哥、印淞、青哥在生活上研究上互相陪伴、打氣與成 長。感謝研究室學弟子榜的幫忙，及學弟妹心瑜、家桂、紹禎、元俊、承昌與建 評，很高興認識你們。感謝男朋友世民的包容、支持、鼓勵與陪伴。也感謝助理 思研以及許多朋友。. 最後謹將本論文以及我最深的謝意獻給我摯愛的雙親。. iv.
(7) 目. 錄. 中文摘要……………………………………………………………………………. i. Abstract………..……………………………………………………………...…… ii 誌謝...…………………………………………………………………..……...…... iv. 目錄..…………………………………………………………………..………….... v. 表目錄…………………………………………………………………...…….…. viii. 圖目錄……………………………………………………………………..……….. x. 第一章. 前言…………………………………………………..…….………..…. 1. 1-1. 研究動機及目的……………………………………..…………………. 1. 1-2. 文獻回顧……………………………………………..…………………. 2. 1-3. 研究方法及流程………………………………………..………………. 3. 1-4. 論文架構………………………………………………..………………. 6. 第二章. 資料介紹..………………………………………………………………. 7. 2-1. 全球的地球重力模型………………………………………...…………. 7. 2-1-1 GRACE…………………………………….…….…………………. 8. 2-1-2 GGM01…………………………………………..…………………. 8. 2-1-3 GGM02……………………………………....…..…………………. 9. 2-1-4 GGM02S、GGM01S 與 EGM96……………….……….………… 11. 2-2. 2-1-5 EIGEN-CG01C…………….…………………..……………….…. 13. 重力資料………………………………….……………………………. 14. 2-2-1 陸測重力……………………….…………………………………. 14. 2-2-2 船測重力………………………………………….………………. 15. v.
(8) 2-2-3 衛星測高……………………………………………….…………. 15. 2-2-4 空載重力……………………………………………….…………. 15. 2-2-5 資料前處理……………………………………………….………. 16. 2-3. 台灣數值高程模型.................................................................................. 18. 第三章. 大地起伏計算原理…………………………………………….………. 19. 3-1. 重力大地起伏計算原理簡介…………………………………….……. 19. 3-2. 去除-計算-回復法.………………………………………………….…. 21. 3-3. 全球大地位模式.…………………………………………....…………. 24. 3-4. 剩餘地形模型理論.……………………………………………………. 26. 3-5. 最小二乘配置法…………………………………………………….…. 29. 3-6. 擬大地水準面………………………………………………...………. 33. 3-7. 重力大地起伏的修正…………………………………………….……. 36. 第四章. 大地起伏成果分析…………………….………………………………. 38. 4-1. 檢核路線…………………………….…………………………………. 38. 4-2. 長波長大地起伏…………………….…………………………………. 39. 4-3. 中、短波長大地起伏……………….…………………………………. 43. 4-4. 改正項.…………………………………………………………………. 45. 4-5. 台灣重力大地起伏模型…………….…………………………………. 47. 4-6. GPS/水準資料修正重力大地起伏….…………………………………. 50. 第五章. 大地起伏模型應用於 DEM 之製作……………….…….….…………. 54. 5-1. 5-2. DEM 製作介紹.………………………………………………...……… 55 5-1-1 航測空三……………………….…………………………………. 5 5. 5-1-2 光達…..……………………………..……………….……………. 57. 航測空三資料及測試.…………………………………………………. 57. 5-2-1 資料來源及分佈..………………………… ………..…………… 57. vi.
(9) 5-2-2 控制測量……………………….…………………………………. 58. 5-2-3 高程控制方法..…………………………………………...………. 60. 5-2-4 精度評估……………………….…………………………………. 62. 5-2-4-1. 控制點部分..…………………………………...…………. 62. 5-2-4-2. 連接點部分..………………………………...……………. 63. 光達資料及測試..……………………………………………...………. 72. 5-3-1 資料來源……………………….…………………………………. 72. 5-3-2 精度評估……………………….…………………………………. 72. 5-3-2-1. 檢核點..………………………………………...…………. 72. 5-3-2-2. DEM 資料..……………………………………..…………. 73. 結論與建議…………………………………………………………….. 76. 參考文獻……………………………………………………………….…………. 79. 作者簡歷……………………………………………………………….…………. 83. 5-3. 第六章. vii.
(10) 表. 表 4-1. 目. 錄. EGM96 之大地起伏檢核成果統計表……………………………...…….... 41. 表 4-2 GGM01S 結合 EGM96 之大地起伏檢核成果統計表…………………... 41. 表 4-3 GGM01C 結合 EGM96 之大地起伏檢核成果統計表………...……….... 41. 表 4-4 GGM02S 結合 EGM96 之大地起伏檢核成果統計表…………………... 41. 表 4-5 GGM02C 結合 EGM96 之大地起伏檢核成果統計表……...………….... 42. 表 4-6. CG01C 之大地起伏檢核成果統計表………………….…………………... 42. 表 4-7. EGM96 經 RCR 計算之重力大地起伏檢核成果統計……………..…….... 49. 表 4-8 GGM01S 結合 EGM96 計算之重力大地起伏檢核成果統計…..….……... 49. 表 4-9 GGM01C 結合 EGM96 計算之重力大地起伏檢核成果統計……...….... 49. 表 4-10 GGM02S 結合 EGM96 經 RCR 計算之重力大地起伏檢核成果統計…... 49. 表 4-11 GGM02C 結合 EGM96 經 RCR 計算之重力大地起伏檢核成果統計…... 50. 表 4-12 CG01C 經 RCR 計算之重力大地起伏檢核成果統計...……….…………. 50. 表 4-13. 台灣大地起伏模型檢核成果統計………………..…………………….... 52. 表 5-1. 23 個高控點之實測橢球高級正高資料………………………………….... 62. 表 5-2. 內差本文大地起伏模型得 23 個高程控制點位置之大地起伏植………. 63. 表 5-3. 實測大地起伏與本文大地起伏模型之差異統計值….…………………... 63. 表 5-4. 各模型使用之控制點及控制類別……………………….………………... 65. 表 5-5. 各模型使用之高程系統及平差方法………………………...……….….... 66. 表 5-6. 空中三角平差之 A、B 橢球高以大地起伏模式改算為正高與空中三角平差 之 C、D、E、F 正高成果比較……………………………………………….. 68. 表 5-7. 檢核點正高與檢核點橢球高經大地起伏改算之正高差值統計….……... 73. 表 5-8. 檢核點橢球高與 Lidar DEM 橢球高之差值統計………….……………... 74. viii.
(11) 表 5-9. 檢核點正高與 Lidar DEM 橢球高經大地起伏改算為正高之差值統計..... ix. 75.
(12) 圖. 目. 錄. 圖 1-1. 大地起伏計算流程…………………………………………………….…... 4. 圖 1-2. 航測 DEM 製作之高程控制流程…………………………………….…... 5. 圖 2-1. 估計全球重力模型 TEG4、EGM96、GGM01S 的階數變方及誤差階數變 方…………………………………………………………………………. 10. 圖 2-2. 估計 EGM96、GGM02S、GGM02C 的階數變方、誤差階數變方………... 10. 圖 2-3. GGM02C 與 GGM02S、EGM96 的頻譜差異……………………………... 10. 圖 2-4 GGM02S 模型在 70 階數次數前之大地起伏誤差………………………... 12. 圖 2-5. EGM96 與 GGM01S 在 70 階數次數前之大地起伏誤差............................. 12. 圖 2-6. 重力異常值……………………………………………………………….... 17. 圖 2-7. 台灣及周圍海域地形圖……………………….…………………………... 18. 圖 3-1. 重力法大地起伏計算理論發展…………………….……………………... 21. 圖 3-2. 本文大地起伏模型計算流程……………………….……………………... 23. 圖 3-3. 剩餘地形模型………………………………………….……. ……………. 27. 圖 3-4. telluroid 面…….……………………………………….……………….…... 34. 圖 4-1. 四條檢核路線上檢核點之分佈及其大地起伏值………………..……….. 38. 圖 4-2. 長波長大地起伏…………………………………………….……………... 40. 圖 4-3. 中波長大地起伏………………………………………………….………... 44. 圖 4-4. 短波長大地起伏值……………………………………………………….... 46. 圖 4-5. 改正項………………….………………………….………...……………... 46. 圖 4-6. 重力大地起伏模型…………………………….…………………………... 48. 圖 4-7. 幾何大地起伏與重力大地起伏差值之諧和面….………………………... 51. 圖 4-8. 台灣大地起伏模型……………………………….………………………... 52. x.
(13) 圖 5-1. 航空攝影測量法生產 DTM、DSM 流程………………………………….. 56. 圖 5-2. 航測空三連接點分佈圖………………………………………………….... 58. 圖 5-3. 平控點 53 點點號及分佈圖…………………………………….………….. 59. 圖 5-4. 高控點 23 點點號及分佈圖………………………………………………... 59. 圖 5-5. 航測空三之高程控制測量方法……….…………………………………... 61. 圖 5-6. A、B 模型使用之控制點…………..……….…………….……………….. 67. 圖 5-7. C、D 模型使用之控制點…….………………………………………….... 67. 圖 5-8. E、F 模型使用之控制點…..……………….…………….………………... 68. 圖 5-9. A、B 橢球高分別減去 C、D、E、F 正高所得之大地起伏值…………….... 70. 光達檢核點周圍一百公尺光達資料分布…..………………………….... 74. 圖 5-11 以 TerraScan 繪製之高程分佈圖……………………………………….... 75. 圖 5-10. xi.
(14) 第一章 前 言. 1-1. 研究動機及目的. 在物理大地中，大地水準面具有重要物理意義，除了將公分甚至公釐級精度 GPS 橢球高轉換為正高成為可能以外，在如研究長時期的垂直地殼變動也需要大 地水準面與水準測量資料、GPS 量測結合才能達成 （Kuroishi et al., 2002） 。 大地起伏的主要貢獻來自於全球大地位模式，依據 Tscherning /Rapp 模型，可 知大地起伏的頻譜（spectral component）訊號組成，0 階～36 階組成 99.479％，37 階～180 階組成 0.479％，181 階～360 階組成 0.018％ （ Basic′ and Rapp, 1992） 。. 即全球大地位模式主導了大地起伏的趨勢，因此為決定台灣區域性大地起伏模 型，一合適之地球重力模型具有舉足輕重的影響，本研究針對 EGM96 模型、目前 GRACE 衛星任務最新的 GGM01S（120 階） 、GGM01C（200 階） 、GGM02S（160 階） 、GGM02C（200 階）分別結合 EGM96 之高階項係數、及結合 CHAMP、GRACE 衛星任務資料之 EIGEN-CG01C（360 階）重力模型製作大地起伏模型，以探討全 球大地位模型與台灣區域性大地起伏模型精度之關係。 另外，航空攝影測量現行製作地形圖的方法為：測量控制點之橢球高及正高， 獲得控制點之大地起伏後，經模式或軟體計算一大地水準面，利用之改正所有連 接點之橢球高，成為適合一般民生使用之正高系統，再進行空中三角、立體測繪 後獲得正高系統數值高程模型（Digital Elevation Model, DEM） 。由於 DEM 的高程 需採正高系統，考量水準測量耗費人力、成本問題，及 GPS 測量的高精度與便利、 大地起伏模型亦已達到公分級精度等優勢，認為於橢球高系統製作 DEM，再經大 地起伏模式改算為正高系統 DEM，似乎已成為一可行且具經濟效益之作業方法（史 天元，2005）。. 1.
(15) 1-2. 文獻回顧. 大地起伏的計算方法主要分為三類，第一類為天文大地法（astrogeodetic） ，此 法用垂線偏差與大地起伏的關係（Heiskarnen and Moritz, 1979）求解，國內 Hwang （1998） 、Hwang et al.（2002）導得一 deflection-geoid 公式，即採此原理計算大地 起伏；第二類為重力法，此法應用重力異常與大地起伏關係 （Heiskarnen and Moritz, 1979） 求解，國內崔國強（1994）採用 Stokes 公式，Hwang（1997）、黃金維等 （2001）、黃金維（2003）、黃金維等（2005）採用最小二乘配置法，均獲得台灣 重力大地起伏模型；第三類為幾何高差法，為水準點上實測之高精度橢球高與正 高差值，國內陳國華（2004）採此原理，並應用於改進重力大地起伏模型精度。 重力大地起伏之研究在世界上行之已久，部分國家亦發展出其法定之區域大 地起伏模式，如美國 G96SSS、G99SSS 採用 Stokes’積分及一維快速傅立葉轉換 （fast fourier transformation, FFT）法計算，並分別以最小二乘配置法（least square collocation, LSC）與 GPS/水準結合為 GEOID96、GEOID99 模型；澳洲 AUSGeoid98 大地起伏使用 Stokes’積分；日本 JGEOID2000 採 Stokes’積分，並與 GPS/水準資料 利用 LSC 結合為 GSIGEO2000 等。另外如義大利 IGeS （International Geoid Services）網站蒐集全球及區域性的大地起伏模型供各界使用、並舉辦訓練及提供 軟體等服務。 區域重力大地起伏計算方法，最基本的理論為採用去除-計算-回復法，配合 deterministic 方法或 stochastic 方法計算大地起伏。另近十數年在大地文獻中提出改 善 deterministic 方法中的 Stokes 公式，主要分為兩類：一類為 deterministic modification of Stokes’法，目的為減少因忽略球帽外範圍之重力產生之誤差;另一類 為 stochastic modification of Stokes’ 法 ， 目 的 為 減 少 全 球 大 地 位 模 式 （ global geopotential model, GGM）及重力資料的誤差（Ellmann, 2005）。. 2.
(16) 1-3 研究方法及流程. 本文大地起伏計算流程如圖 1-1，乃採去除計算回復法（remove-compue-restore, RCR） ，並聯合陸測重力（gravity data）、船測重力（gravity data） 、衛星測高化算 資料（gravity data）、地球重力模型（earth gravity model, EGM）、及台灣數值高程 模型（DEM） ，以擾動位與重力異常及擾動位與大地起伏關係（GGM） 、剩餘地形 模型（residual terrain model, RTM）理論、最小二乘配置法（LSC）平差計算，可 獲得一擬大地水準面（quasigeoid） ，再經一與高程相關之地形效應（terrain effect） 改正項，獲得台灣區域性重力大地起伏模型（Taiwan local gravimetric geoid model），最後以 GPS/水準資料修正重力大地起伏模型成為成為台灣區域性大地起 伏模型（Taiwan local geoid model）。. 在 DEM 製作之應用方面，根據橢球高（h）、正高（H）及大地起伏（N）之 線性關係（Heiskanen and Moriz, 1979） h−H −N =0. （1-1）. 本文在航空攝影測量資料採用內政部五千分之一像片基本圖測製成果，包括 控制點之實測橢球高、正高資料，及連接點之空中三角平差成果；比較現行方法-正高系統測製 DEM 成果，其流程如圖 1-2 左圖，與建議方法--橢球高系統測製 DEM 成果再輔以大地起伏模型改算後之正高系統 DEM，其流程如圖 1-2 右圖，方形虛 線範圍內為橢球高高程值，橢圓虛線範圍內為正高高程值，以瞭解現行作業成果 與建議作業即輔以大地起伏模型進行高程轉換成果之關係。在光達方面，採用工 業技術研究院能源與資源研究所資料，包括檢核點實測橢球高、正高資料及光達 橢球高系統點雲資料。以檢核點資料進行大地起伏模型之檢核，並以檢核點資料 檢核光達成果，及將光達點雲資料進行橢球高 DEM 製作後輔以台灣大地起伏模型 以獲得之正高 DEM。. 3.
(17) 圖 1-1、大地起伏計算流程. 4.
(18) 圖 1-2、航測 DEM 製作之高程控制流程，左圖為現行之高程控制流程，右圖配合 採用大地起伏模型之高程控制流程. 5.
(19) 1-4 論文架構 本論文共分六章，各章之內容安排如下： 第一章：前言，說明本論文的研究動機及目的、文獻回顧、研究方法及流程及本 文架構。 第二章：介紹計算大地起伏所需資料，包括全球大地位模式所需之 GRACE 衛星任 務最新 GGM01、GGM02 模型及結合 CHAMP 及 GRACE 衛星任務資料 之 CG01C 模型、剩餘地形效應所需之台灣數值高程模型、及最小二乘配 置法計算所需之觀測重力資料，來源包括陸測、船測及衛星測高化算成 果。 第三章：說明本文大地起伏模型計算理論，首先簡介各種重力大地起伏計算理論， 接著介紹本文採用之去除計算回復法、全球大地位模式理論、剩餘地形 效應理論、最小二乘配置法原理、地形效應改正項及 GPS/水準資料之加 入計算。 第四章：展示不同長波長、中波長、短波長、重力大地起伏模型之差異，及模型 檢核統計。檢核資料來源由內政部提供，包括了海岸區、高山區、平原 區共計四條檢核路線。最後提出台灣最新區域性重力大地起伏模型及台 灣最新區域性大地起伏模型。 第五章：簡介國內現行航空攝影測量於正高系統製作地形圖方法，及其於橢球高 系統製作之 DEM 經大地起伏模型轉換為正高系統 DEM 成果差異；並採 用光達資料進行 DEM 製作及檢核。 第六章：總結本研究之結論與建議。. 6.
(20) 第二章 資料介紹. 大地起伏的計算，理論上需要全球重力資料，由於實務上執行並不容易，因 此計算區域性重力大地起伏時，就將高階數次數之球諧函數模型（spherical harmonic model, SHM） （即文中 2-1 節模型）代替計算區域周圍之重力場資料，並 考慮高頻的地形效應，此區域性大地起伏計算理論需要使用之資料來源敘述如下。. 2-1 全球的地球重力模型（Global Earth Gravity Model）. 自 1977 年以來，世界上已經發展出多套全球的地球重力模型，其發展請見崔 國強（1994） ，這些模型亦已運用於大地起伏計算，如澳洲前國家製圖部（Division of National Mapping） ，現為澳洲測量及土地資訊組織（Australian Surveying and Land Information Group），公布之全國重力大地起伏模型 AUSGeoid91 即採用 OSU89A、之後又使用 OSU91A 計算出 AUSGeoid93 模型（Featherstone et al., 2001），美國的 G96SSS、G99SSS 則採用 EGM96 模型（Smith and Milbert, 1999 ; Smith and Roman, 2001）。台灣亦利用 EGM96 於 2001 年（黃金維等，2001）配合. 30′′ × 30′′ 解析度 DEM，計算得一 30′′ 解析度大地起伏模型，並於 2003 年（黃金維， 2003）配合 3′′ × 3′′ 解析度 DEM，計算得一 3′′ 解析度大地起伏模型，較新的模型包括 EGM96（詳見 Lemoine et al., 1998）、GGM01、GGM02、EIGEN-CG01C。本文則針對 EGM96 及最新的 GGM01、GGM02、CG01C 模型進行計算比較及分析。. 另 外 ， 美 國 地 理 空 間 資 訊 社 （ National Geospatial-Intelligence Agency,NGA）正在籌備發展 EGM05，EGM05 將採用 CHAMP 與 GRACE 衛星的長波長資 訊、高程資料來自 NASA/NGA 的太空梭雷達製圖任務（shuttle radar topographic mission, SRTM）任務、NGA 持續蒐集的全球重力資料，及採用較佳模型及估計技. 7.
(21) 術，並希望提供一目前公佈模型尚未達到的高階數次數地球重力模型（Kenyou and Pavlis, 2004）。. 2-1-1 GRACE. 美國太空總署（National Aeronautics and Space Administration,NASA） 地球科學先導組（ESSP）在 1997 年 5 月計畫GRACE （gravity recovery and climate experiment）任務，該任務開始於 2002 年 3 月，利用兩個相距 220 公里的飛行器 在距地球 500 公里的高空上運行，生命週期預計 5 年，功能為精確地量測地球重 力場，預計獲得之成效含地球內部質量及環繞地球質量的分佈及流動 （http://www.csr.utexas.edu/grace/overview.html）。. 2-1-2 GGM01. GGM01 模型包含了 GGM01S、GGM01C 兩模型。 GGM01S由 2002 年 4 月至 11 月間共 111 天的GRACE的K波段射程方向(range rate)、姿態、及加速計資料所計算得。其成果沒有加上Kaula約制、其他衛星資 訊、表面重力資料、也沒有加上其他先處理的資訊。GGM01S場估計到 120 的階數、 次. 數. ，. 95. 階. 數. 前. 成. 果. 正. 確. 性. 高. （. 圖. 2-1. ，. http://www.csr.utexas.edu/grace/gravity/ggm01/），因此建議高於 95 階數要 謹慎使用。基礎的靜態場模型是TEG4、海潮模型是由CSR4.0 基於選擇TEG4 重力解 的共鳴共振潮、固體地球潮模型則與IERS96 彈性地球模型一致。同時也解算下列 參數：1、每天弧形的初始狀態，2、加速計每天的偏差及全球比例因子，3、KBR 偏差、GPS未定值、天頂延遲。 GGM01S與TEG4 資訊(information)方程（由歷史多衛星追蹤資料、表面重力 資 料 、 測 高 之 海 面 高 程 製 造 ） 結 合 產 生 階 數 和 次 數 到 達 200 的 GGM01C 模 型. 8.
(22) （http://www.csr.utexas.edu/grace/gravity/ggm01/）。 2-1-3 GGM02. G G M 0 2 模 型 包 含 了 G G M 0 2 S 及 G G M 0 2 C 兩 模 型 。 GGM02S由 2002 年 4 月至 2003 年 12 月間共 363 天的GRACE K波段射程方向、姿態、 及加速計資料估計得。其成果沒有加上Kaula約制、沒有其他衛星資訊、表面重力 資料、也沒有加上其他先處理的資訊。GGM02S估計到階數和次數到達 160 階，其 成果到 120 階數前均保持正確的訊號能量頻譜（圖 2-2， http://www.csr.utexas.edu/grace/gravity/ggm02/ ） 。同時也解算下列參數：1、每日弧形的初始狀態，2、加速計每日的偏差及 每週的比例因子，3、KBR 偏差、GPS 未定值、天頂延遲。 GGM02C是由GGM02S結合地面重力資訊（包括表面重力及平均海面） ，利用TEG4（到 200 階數）的協變方函數的將高階數部分約制到EGM96 的和諧係數。圖 2-3 （ http://www.csr.utexas.edu/grace/gravity/ggm02/ ） 可發現GRACE衛星資訊主導地球重力模型 110 階以內的解，並與EGM96 結合 120 階以上到 200 階的部分，因此GGM02C的 200 階數以上係數可採用EGM96 係數以擴充 到 360 階。（ http://www.csr.utexas.edu/grace/gravity/ggm02/ ）. 9.
(23) 圖 2-1、估計全球重力模型 TEG4、EGM96、GGM01S 的階數變方（degree variances） 及誤差階數變方（degree error variances），（單位：公釐）. 圖 2-2、估計 EGM96、GGM02S、GGM02C 的階數變方及誤差階數變方(單位:公釐). 圖 2-3、GGM02C 與 GGM02S、EGM96 的頻譜差異（spectral differences）. 10.
(24) 2-1-4 GGM02S、GGM01S 與 EGM96. GRACE任務的地球重力模型強調全球均勻的資料分佈，因此大地起伏誤差在陸 地. 及. 海. 域. 並. 沒. 有. 明. 顯. 界. 線. ，. 如. 圖. 2. -. 4. （ http://www.csr.utexas.edu/grace/gravity/ggm02/ ） 為GGM02S在 70 階前（空間解析度約 300 公里）的大地起伏誤差、圖 2-5 （http://www.csr.utexas.edu/grace/gravity/ggm01/）為GGM01S、EGM96 在 70 階前的大地起伏誤差，圖 2-5 中EGM96 的誤差即可看出明顯在海、陸的差異。GGM02S 在 70 階數次數大地起伏誤差的全球均方根誤差約 7 公釐，最大誤差為 9 公釐。 GGM01S在 70 階數次數時的大地起伏誤差約 2 公分，90 階數次數（空間解析度約 2. 0. 0. 公. 里. ）. 時. 的. 誤. 差. 約. 6. 公. 分. （ h t t p : / / w w w . c s r . u t e x a s . e d u / g r a c e / g r a v i t y / g g m 0 1 / ）， （ h t t p : / / w w w . c s r . u t e x a s . e d u / g r a c e / g r a v i t y / g g m 0 2 / ）。. 11.
(25) 圖 2-4、GGM02S 模型在 70 階數次數前之大地起伏誤差（單位：公分）. 圖 2-5、EGM96 與 GGM01S 在 70 階數次數前之大地起伏誤差（單位：公分）. 12.
(26) 2-1-5 EIGEN-CG01C. 目前德國波茨坦GFZ（GeoForschungsZentrum Potsdam）製作之EIGEN-CG01C 模型，利用GRACE衛星任務 200 天（2002 年 4 月、5 月、8 月、11 月及 2003 年 4 月、5 月、8 月、10 月、11 月）、CHAMP衛星任務 860 天（2000 年 10 月至 2003 年 6 月）的資料，加上 0.5°× 0.5° 地面重力資料計算而得之高精度高解析度之全球 重力場模型，解算至 360 階。結合方法採應用一個特別的帶寬組合方法以便保存 從衛星資料在低頻帶的高精度，並利用平滑技術取得自地面重力資料的高頻資訊 （王成機，2005）。用 130317 個球諧係數表示地球重力位，本模型之全球大地起 伏 及 自 由 空 間 重 力 異 常 的 空 間 特 徵 可 顯 示 到 100 公 里 的 全 波 長 範 圍 。 （http://www.gfz-potsdam.de/pb1/op/grace/results/grav/g003_eigen-cg01c.html#t1）. 另外EIGEN-CG03C亦使用與EIGEN-CG01C相同之CHAMP資料，及GRACE376 天的 資料（2003 年 2 月～5 月及 7 月～12 月、2004 年 2 月～7 月） ，EIGEN-CG03C亦展 開. 到. 360. 階. 數. 次. 數. （ http://www.gfz-potsdam.de/pb1/op/grace/results/grav/g004_eigen-cg03c .html）。本文並未採用此模型進行大地起伏計算。. 13.
(27) 2-2 重力資料. 重力資料來源包括：陸測重力、船測重力、衛星測高資料。. 2-2-1 陸測重力. 陸測重力資料來源涵蓋： 1、中央研究院地球科學研究所(Yen et al., 1990; Yen et al., 1995)，以 LCR-G 型 重力儀，歷時 7 年於 1987 年完成之 603 個重力點，其座標系統為 TWD67， 重力異常為 GRS67 系統。 2、中國測量學會，以 LCR-G 型重力儀，歷時 3 年於 1988 年完成之 276 個重 力點（黃金維等，1998），其中 2 點位於海面為明顯錯誤典故予以剔除， 並剔除金門、馬祖、澎湖等地區資料，其座標系統為 TWD67、重力異常為 IGF30 系統；及歷時 3 年於 1999 年完成之 747 個重力點，座標系統為 TWD67、重力異常為 GRS67 系統。 3、內政部，以 LCR-G 型重力儀，在一等一級水準點上施測，歷時 2 年於 2001 年完成之 6 個絕對重力點、10 個一等重力點、1020 個一等一級水準點上 重力值，其座標系統為 TWD97（黃金維，2001）；在一等二級水準點上， 歷時 2 年於 2003 年完成之 1092 個一等二級水準點上重力值（陳春盛， 2003），其重力異常為 GRS80 系統。 陸測資料共 3715 點。上述之重力資料由於採用之座標系統及參考之橢球系統 不同、系統偏差等誤差，因此需要改正為統一之基準，主要改正為將實測重力值 經自由空間約化至海水面、從 IGF30、GRS67 系統改正為 GRS80 系統 （Torge, 1989） 、大地改正、垂直及水平基準面不一致改正，改正部分詳見郭重言（1998）。. 14.
(28) 2-2-2 船測重力. 中央大學地球科學系教授許樹坤進行重力、磁力測量及在 Hsu et al.（1998） 中整理及修編，共約 4930 點。重力改正包括重力儀的飄移、重力參考場不一致、. 效應，改正部分詳見郭重言（1998）。 基站偏移、 Evotov. 2-2-3 衛星測高. 衛星測高來源有二，其一來自交通大學，為 Seasat、Geosat/GM、Geosat/ERM、 ERS-1/GM、ERS-1/35d、ERS-2/35d、TOPEX/POSEIDON 衛星資料，使用觀測量為海 水面高度（sea surface height, SSH） ，SSH 中包含了大地起伏、海水面地形（sea surface terrain, SST）、軌道徑向誤差、改正模式誤差、偶然誤差，改正部分見 郭重言（1998） 。利用大地水準面梯度理論消除軌道徑向誤差後，以去除回復法得 之殘餘大地梯度配合最小二乘配置法組成東西、南北向量，並利用 inverse Vening Meinesz 公式及 1D FFT 計算技巧求得重力異常。詳細資料及理論請見徐欣瑩（1997）; 其二來自丹麥國家測量局（National Survey and Cadastre, Denmark. Kort & Matrikelstyrelsen, KMS）在 2002 年計算得之重力資料 KMS02，詳細製作細節請 見 Andersen and Knudsen（1998） ，共約 10000 多點。重力改正包括地心偏移誤差、 每轉軌道誤差、海潮改正模式誤差、海水面地形等。. 2-2-4 空載重力. 內政部於 2005 年委託交通大學完成空載重力測量工作，空載重力資料位於空 中約 5156 公尺處，使用時採取向下延續法（downward continuation）將重力資 料計算到地表面或海平面再利用。空載重力之優點為不受台灣高山起伏之地形限 制，且可以快速獲取大範圍之重力資料，相對上可節省大量時間、金錢與人力。. 15.
(29) 本文計算大地起伏時並未採用空載重力資料，僅在陸測、船測重力異常資料 的前處理時作為比較基準，以找出陸測、船測資料的大錯點位予以剔除。. 2-2-5 資料前處理. 上述資料大部分位於台灣本島，小部分位於金門、馬祖、澎湖、綠島、小琉 球、龜山島等，另外海測資料、衛星測高資料均在北緯 21° ~ 26° 、東經 119° ~ 123° 範 圍 ， 由 於 第 3 節 之 台 灣 數 值 高 程 模 型 僅 涵 蓋 北 緯 21.5° ~ 25.5° 、 東 經. 119.5° ~ 122.5° 範圍，故實際使用資料為陸測重力 3559 點、船測重力 2933 點、交 通大學衛星測高資料 12575 點、丹麥國家測量局衛星測高資料 5471 點。. 本文參考內政部 2005 年空載重力測量工作之空載重力資料，消除陸測、船測 有明顯錯誤之重力點資料。利用衛星測高資料與海測資料相比，由於衛星測高資 料一般而言不會有大錯，因此可藉此趨勢相同的地區找出海測資料錯誤。並利用 空載重力資料與陸地資料相比較，亦是利用重力變化趨勢相同的原理找出陸測資 料之大錯。. 陸測、船測、衛星測高、空載測量原始重力異常值如圖 2-6，如圖 2-6（a） 與（c）相比可發現圖 2-6（a）在經度 121 度與緯度 23 度交叉處有一明顯黑點， 比對圖 2-6（c）認為其為大錯點，予以刪除。共計發現錯誤點計陸測重力 8 點， 船測重力 4 點，共 12 點。. 16.
(30) （a）. （b）. （c） 圖 2-6、重力異常值（單位：mgal） （a）陸測、船測重力異常值、 （b）衛星測高重 力異常值、（c）空載測量重力異常值. 17.
(31) 2-3 台灣數值高程模型. 本研究中使用之數值地形資料為規則網格，由行政院農委會林務局農林航空 測量所提供，原始資料為航空照片、平面解析度 40 公尺。經過萃取後成為平面解 析度 80 公尺網格資料，共計 6421075 點，座標系統為 TWD67。利用最小曲率法 （minimum curvature method）製作網格大小為 3 秒之數值高程模型，如圖 2-7 左，再利用一、二等衛星控制點及一等一、二級水準點採用諧和面法（harmonic surface method）修正，最後與國家海洋科學研究中心提供之 2′ × 2′ 海底地形網格 組成為以台灣為中心，北緯 21.5° ~ 25.5° 、東經 119.5° ~ 122.5° 範圍之 DEM 規則網 格資料。詳細資料介紹及資料處理步驟請參閱黃金維（2003） 、呂誌強(2004)，為 後續應用並將此 DEM 資料製作成 3′ × 3′ 參考 DEM (使用 GMT 軟體，採用 median filter, 50km)，如圖 2-7。. 圖 2-7、台灣及周圍海域地形圖（陸上用農航所 DEM，海上用國家海科中心 DEM） ， 左圖為原始 DEM，右圖為參考高程面（解析度 3′ × 3′ ），（單位：公尺）. 18.
(32) 第三章 大地起伏計算原理. 陸地上大地起伏模型的計算主要有三種方法：重力法（gravimetric）、天文大 地法（astro-geodetic） 、幾何高差（橢球高減去正高）法（Smith and Milbert, 1999）。. 重力法求區域大地起伏主要分為，一、deterministic method，代表性的方法為 採用 Stokes’ 積分配合一維球面快速傅立葉計算（1D. FFT）或二維平面快速傅立. 葉計算（2D FFT） ；二、stochastic method，代表性的方法為採用最小二乘配置法， 其優勢為可結合不同之重力資料來源作計算。. Deterministic 型中，自 1849 年 Stokes 發展 Stokes’ 積分公式以來，已廣泛應 用於大地起伏的計算；使用 Stokes’ 公式需要全球的重力資料，然而要獲得全球重 力資料實行上不甚容易，因此之後許多學者作了修正，請見 3-1 節。最小二乘配置 法為 stochastic 型的代表性方法，其最大的優點為能結合不同來源的資料，且不需 將觀測資料作成網格，缺點為協變方矩陣會因資料量太大造成求逆時的不穩定 （陳俊廷，1997）。. 為了計算區域大地起伏，去除計算回復法是最廣泛被採用之方法 （Featherstone et al., 2004），此法乃結合最高階數次數之全球大地位模式及區域重 力資料、高程資料作計算，其中全球大地位模式代表的是大地起伏的低頻或長波 長部分，而重力資料及高程資料分別代表中頻（中波長）及高頻（短波長）的大 地起伏效應。. 3-1. 重力大地起伏計算理論簡介. 圖 3-1 為重力法大地起伏計算理論發展的流程，主要分為 deterministic method 19.
(33) 及 stochastic method。因 Stokes 公式需要全球重力資料，後來分為由 Molodensky 及 Helmert 作了改善。Molodensky 公式乃決定擬大地水準面（quasigeoid）；Helmert 公式則決定大地水準面。其中 Molodensky 公式，又由 Moritz 及 Bjerhammar 又分 別將其導得近似公式使用。. 若具備全球重力資料可採 Stokes 公式計算大地起伏。然而實際計算上，僅用 以計算點為中心之球帽積分，而造成截去球帽範圍外積分之誤差。首先由 Molodensky 提出為減低該誤差將 Stokes 公式作修正（Ellmann, 2005） 。為了計算方 便，在應用 Molodensky 公式，Moritz 用解析延續（analytical continuation）法修改 之；Bjerhammar 亦採用一次（order）垂直梯度近似修改之（ Sjoberg , 2005）。. Helmert 第二壓縮理論乃將大地水準面外質量壓縮在大地水準面上，可決定大 地起伏高，後由 Sideris 導得一近似公式（ Sjoberg , 2005）。. 近數十年大地文獻主要針對 Stokes 公式的修改方法分為兩類，一類為 deterministic modification，此法可大幅度減低由於忽略球帽範圍外遠處區域的重力 所造成的影響，不過沒有結合大地位係數及重力資料的精度分析，雖然這兩者的 誤差其實會傳播到大地起伏成果；另一類為 stochastic modification，此法目的為減 低全球大地位模式及重力資料的誤差，並可發現在有些模型中亦降低了球帽範圍 外重力影響的誤差（Ellmann, 2005）。. 20.
(34) 圖 3-1、重力法大地起伏計算理論發展. 去除-計算-回復法（remove-compute-restore, RCR）. 3-2. 目前全球最廣為使用於決定區域重力法大地起伏模型的方法，即為去除計算 回復法（Featherstone et al., 2004） ，即利用觀測之重力異常，去除利用全球大地位 模式理論所計算得之長波長重力異常，再去除利用剩餘地形效應理論所計算得之 短波長重力異常，其值稱為殘餘重力異常（中波長重力異常） 。將殘餘重力異常經 最小二乘配置法計算後，獲得殘餘大地起伏值（中波長大地起伏值） ，最後再回復 長波長大地起伏及短波長大地起伏，若以擾動重力位表示為 T = TGGM + TRTM + Tres. （3-1）. 式中. 21.
(35) TGGM ：全球大地位模式的擾動重力位. TRTM ：剩餘地形模型理論的擾動重力位 Tres ：殘餘的擾動重力位，即未模型化的殘餘重力場. 本文使用 RCR 法之流程如圖 3-2，圖中上面的虛線部分為 RCR 法中的去除部 分，其將重力資料去除長（輸入資料為萬有引力常數、地球質量、參考橢球體之 長半徑及完全正規化之地球引力位球諧係數；用 gmod.f 程式，請見 3-3 節） 、短波 長重力異常（輸入資料為台灣數值高程模型及台灣數值參考高程模型；用 tcfour.f 程式，請見 3-4 節） ，接著使用 3-5 節中 Tscherning/Rapp 模型理論（輸入資料為萬 有引力常數、地球質量、參考橢球體之長半徑、完全正規化之地球引力位球諧係 數及其誤差；用 errvar.f 程式）、協變方矩陣（用 f300.f 程式），採用最小二乘配置 平差法（用 collocg.f 程式）計算得殘餘大地起伏。在圖中下面的虛線部分為 RCR 法的回復部分，其將殘餘大地起伏與長（用 gmod.f 程式，請見 3-3 節） 、短波長大 地起伏（用 tcfour.f 程式，請見 3-4 節）之網格資料用普通製圖工具（generic mapping tools, GMT）軟體相加而得一擬大地起伏模型，再經一與地形面高程相關之地形效 應改正項（請見 3-6 節） ，可獲得區域性重力大地起伏模型，最後採用 GPS 之橢球 高及水準測得之正高資料的差值，修正此區域性重力大地起伏模型（請見 3-7 節） ， 即可得台灣區域大地起伏模型。. 22.
(36) 圖 3-2、本文大地起伏模型計算流程 23.
(37) 3-3. 全球大地位模式（Global Geopotential Model, GGM）. 全球大地位模式是基於物理大地學門中之地球重力場理論，重力場理論是基 於一不隨時間變動的地球重力場的假設。 地球重力位 W(r)函式為地球引力位 V(r)與離心力位 Z(r)之和。 W(r)=V(r)+Z(r). （3-2）. 式中 V(r)= G ∫ ∫ ∫ v. Z(r)=. ρ (r / ) | r/ − r |. （3-3）. dv. w2 2 d 2. （3-4）. 式中 G：萬有引力常數，6.673x 10−11 m3kg −1s −2 ，. ρ (r / ) ：體積密度函數， r 、 r / ：地心位置向量，. v ：地球體積。 w ：地球角速度，7.292115x 10−5 rad s −1 ， d ：垂直於旋轉軸之距離。 由上述式子可知離心力位可以很容易地被計算得，但地球引力位受制於密度 函數是未知的，而無法獲得。 然而，全球地球物理密度函式只考慮徑向密度的改變，故當需要到區域性的 計算時，就需要獲得 ρ (r / ) 的資訊方為可能。為了達成此目的，在地球表面或地球 外部的由不同技術來取得觀測量才能導出地球引力位函式（Torge, 1989）。 由建構引力位在邊界面 S 的一階導數 Gauss’積分方程為. 24. ∂V 與二階導數 Vxx = ∂ 2V / ∂X 2 的關係 ∂ns.
(38) ∫∫ S. ∂V = ∫ ∫ ∫ ΔVdv ， v ∂ns. （3-5）. 式中. ΔV ：拉普拉斯微分運算子， 在 XYZ 空間為. ΔV = Vxx + Vyy + Vzz. （3-6）. 地 球 質 量 M= ∫ ∫ ∫ ρ (r / )dv ， 由 位 理 論 可 知 ∫ ∫ V. S. ∂V = −4π GM ， 並 同 時 考 量 ∂ns. Poisson’s 微分方程 ΔV = −4π G ρ (r / ). （3-7）. 及假設大氣質量為零，使得地球表面成為引力來源的地球內部與無質量的外 部空間之邊界面，可得拉普拉斯微分方程 ΔV=0，其解可將全球的重力場 V 展開為 球諧函數，在地球外空間之表示式如下 V (r ,θ , λ ) =. ∞ a n GM {1 + Σ ( ) n ∑ (Cn ,m cos mλ + S n , m sin mλ ) Pn ,m (cos θ )} n=2 r m =0 r. （3-8）. 式中 a：橢球的半長軸 Pn , m (cos θ ) ：聯合勒建德函數. n、m：階(degree)、次(order) Cn , m 、 S n ,m ：球諧係數 r 、 θ 、 λ ：球座標（地心距離、地心經度、極距）. 而擾動位 T 為真實及地球參考橢球體正常引力位 U 之差，可得 T (r ) = W (r ) − U (r ). （3-9）. 又由於地球外空間之拉普拉斯微分方程亦為零，即 ΔT = 0 ，並假設橢球體之 質量與地球質量相同可得擾動位之球諧展開式. 25.
(39) T (r ) =. GM ∞ a n n { ∑ ( ) ∑ (ΔC n ,m cos mλ + Δ S n ,m sin mλ ) P n ,m (cos θ )} r n = 2 r m =0. （3-10）. 式中 ΔC n ,m 、 Δ S n ,m ：真實與正常重力之完全正規化球諧係數之差. P n ,m (cos θ ) ：完全正規化勒建德多項式 利用 Bruns 公式（Heiskanen and Moritz, 1985），可得參考之長波長大地起伏 值為. N=. T. γ. =. GM rγ. nmax. a n n ( ) ∑(ΔC n ,m cos mλ + Δ S n ,m sin mλ ) P n ,m (co s θ ) ∑ n=2 r m =0. （3-11）. 式中. γ ：正常重力 利 用 擾 動 位 T 與 重 力 異 常 Δg 的 關 係 ， 可 得 參 考 之 長 波 長 重 力 異 常 值. Δg = − GM = 2 r. ∂T 2 − T ∂r r nmax. a 2 n n − ( 1)( ) ∑(ΔC n ,m cos mλ + Δ S n ,m sin mλ ) P n ,m (cos θ ) ∑ r m =0 n=2. （3-12）. 此所獲得之大地起伏值，由於球諧係數模型可能會造成系統性誤差（Denker et. al., 2000），另外若以空間解析度探討，空間解析度 Δs =. 180° ，式中 nmax 代表球諧 nmax. 係數展開的階數次數，就現有之 360 階數次數之全球大地位模式，其空間解析度 約 0.5° ，可發現計算出來之精度不足，考量此兩因素，無法僅以全球大地位模式精 確地決定區域大地起伏模型。. 3-4. 剩餘地形模型理論（Residual Terrain Model, RTM）. 由於陸測、船測、衛星測高重力分佈約 2′ ～ 3′ 一點，為獲得高於重力資料解 析度的區域大地起伏模型，我們就需要結合高程資料（Forsberg, 1984）。. 26.
(40) 大地起伏的決定當中，針對地形效應的處理主要有三種計算方法，一、地形 均衡改正，其利用觀測得之重力場參數描述地球深處質量如何均衡補償之主要趨 勢，但實際上，地球上很多區域，尤其在深海溝與中海脊處，在深層並不是一定 處於均衡狀態（Forsberg, 1984）；二、Helmert 壓縮法，其可視為均衡改正理論之 特例，乃將質量沿著垂線方向壓縮至大地水準面上之一表層，造成重力位及大地 水準面改變，此質量並非被移除，此法獲得之重力值是參考於大地水準面；三、 剩餘地形效應改正，則採用一平滑的參考面，利用以點資料製作之地形高程模型 與此參考地形高程模型之差值作計算，RTM 只考量真實地形之短波長效應，即殘 餘值接近零且變方值很小，因此很適合用於最小二乘配置法，需注意的是此法獲 得之重力值是參考於擬大地水準面（Omang and Forsberg, 2000）。. 本研究採用剩餘地形模型理論，目的為模型化地形的高頻效應。RTM 計算重 力異常及高程異常時需參考於一平滑參考平均高程面（如圖 3-3） ，且 RTM 地形改 正能估計相對於平均高程面的質量異常效應，即移去平均高程面上之高山質量及 填滿平均高程面下之山谷質量。. RTM 計算主要有四種方法： 1 、數值法求積分（蔡進雄， 1997 ）， 2 、捲積 （convolution）定理採用快速傅立葉（Fast Fourier Transformation, FFT）積分法 （Forsberg, 1984） ，3、最小二乘配置法（蔡進雄，1997） ，4、prism 積分法（Forsberg,. 1984）。本文用捲積定理採快速傅立葉積分法（Schwarz et al., 1990）計算 RTM 大 地起伏、重力異常。. 圖 3-3、剩餘地形模型. 27.
(41) RTM 之大地起伏與重力異常計算公式如下：. 一、RTM 大地起伏採用公式為(Hwang,1997)：. 經考量地球密度為點位位置之函數，可得. N RTM ( x p , y p ) =. G. γ. ∫. E. ρ ( x, y ) ( h ( x, y ) − href ( x, y ) ). (x − x ) +( y − y ) 2. p. =. 2. dxd y. p. (. ). 1 G 1 ⎡ ρ ( h − href ) ⎤ * = ρ1 ⎡ h − h ref * ⎤ ⎢ ⎣ ⎦ r γ r ⎥⎦ γ ⎣. G. （3-13）. 式中. G：萬有引力常數. γ ：正常重力 E：為 x-y 平面之積分定義域. ρ ：積分元素之地質密度 r = x2 + y2. x p 、 y p ：將求大地起伏之點位座標 ∗ ：捲積運算因子. h 、 href ：分別為積分元素之實際高程與參考面高程. ( h, h ) r. ⎧ ( h, hr ) ：陸地上 ⎪ = ⎨ρ ，稱為比例(scaled)高程，且在海洋上值為負 2 ⎪⎩ ρ1 ( h, hr ) ：海洋上. 二、RTM 重力異常公式為(Hwang,1997)：. 28.
(42) Δg RTM ( x p , y p ) = 2π G ρ1 ( h − href ) − c ( x p , y p ). （3-14）. 式中 c(x p , y p ) =. ⎡ ⎤ 1 1 1⎞ ⎛ Gρ1 ⎢hˆ 2 * 3 − 2hˆ p ⎜ hˆ * 3 ⎟ + hˆ p2 F ⎥ ：若在陸地上稱為地形改正，若 2 r ⎝ r ⎠ ⎣ ⎦. 在海洋上僅代表由真實海深及參考海深所包圍的布格片效應的改正。. 式中. (. ⎧ ( h, hp ) ：陸地上 ⎪ ， h , hp = ⎨ ρ ⎪ 2 ρ ( h.hp ) ：海洋上 1 ⎩. ). 1 dxd y E r3. F=∫. 3-5. 最小二乘配置法（Least Square Collocation, LSC）. 最小二乘配置法可以結合多種資料來源，為此法優點。本文中重力資料的使 用，在衛星測高資料來源，交通大學資料使用水準梯度形式，丹麥國家測繪局資 料使用重力異常形式。 應用在大地起伏之最小二乘配置法首先假設所有觀測量均以海水面為參考 面，假設大地水準面之外沒有質量存在，參考橢球體以球近似當作球體使用。可 得基本方程為. Sˆ = Cst (Ctt + Cnn ) −1 l. （3-15）. 式中. 29.
(43) S：預估值（於此為大地起伏值） Cst ：預估值與觀測值之協變方矩陣 Ctt ：觀測值之間的協變方矩陣 Cnn ：誤差量之間的協變方矩陣. l ：觀測量（於此為重力異常值） 使用最小二乘配置法的先決條件是要決定出觀測量與預估值的協變方矩陣， （3-15）式之 Cst 。一般假設兩點間的協變方矩陣只與球面距離 Ψ 有關。 其次針對（3-15）式之 Ctt 項，及殘餘重力異常的協變方矩陣作一探討。經分 解觀測重力異常的殘餘值， （原本此殘餘值應包括扣除短波長重力異常量，經將短 波長重力異常之誤差由程式 tcfour.f 處理，故假設其值為零）可分解為真實的殘餘 T 重力異常 Δg res 加上觀測重力異常的誤差 nΔgobs 減去長波長重力異常的誤差 nΔg L ，如. 下 obs Δg res. = Δg obs − Δg L = (Δg T + nΔgobs ) − (Δg LT + nΔg L ). （3-16）. = (Δg T − Δg LT ) + (nΔgobs − nΔg L ) T = (Δg res ) + (nΔgobs − nΔg L ). 式中 obs ：觀測重力異常的殘餘值 Δg res. Δgobs ：觀測重力異常. Δg L ：長波長重力異常 Δg T ：真實重力異常 nΔgobs ：觀測重力異常的誤差. Δg LT ：真實長波長重力異常. 30.
(44) nΔg L ：長波長重力異常的誤差 T ：真實的殘餘重力異常 Δg res. （3-16）式之協變方矩陣如下 obs obs cov(Δg res , Δg res ) T T = cov(Δg res , Δg res ) + cov(nΔgobs , nΔgobs ) + cov(nΔg L , nΔg L ). （3-17）. 式中 T T cov( Δg res , Δg res ) ：真實殘餘重力異常的協變方矩陣(訊號值). cov( nΔgobs , nΔgobs ) ：觀測重力異常之隨機誤差的協變方矩陣，即（3-15）式之 Cnn. 項(隨機誤差) cov( nΔg L , nΔg L ) ：長波長模型誤差的協變方矩陣(誤差值) T T （3-17）式即為（3-15）式中括弧內之 Ctt ＋ Cnn 項， Ctt 項為 cov( Δg res , Δg res )與. cov( nΔg L , nΔg L ) 之和， cov( nΔgobs , nΔgobs ) 項為 Cnn 項。. 為了推導（3-17）式，茲先進行理論值（真實重力異常）之推導，首先分解 T 真實重力異常 Δg T ＝ Δg LT ＋ Δg res ，其協變方函式為. cov(Δg T , Δg T ) ∞. = ∑ Cn S n + 2 Pn (cos Ψ ). （3-18）. n=2. ∞. nmax. = ∑ Cn S n + 2 Pn (cos Ψ ) + n=2. ∑ Cn S n + 2 Pn (cos Ψ ). n = nmax +1. 而本文要計算之（3-16）式之觀測之重力異常包含誤差，故將（3-18）式改為 cov(Δg T − Δg L , Δg T − Δg L ) T T = cov(nΔg L , nΔg L ) + cov(Δg res , Δg res ) nmax. = ∑ δ Cn S n +1 Pn (cos Ψ ) + n=2. （3-19）. ∞. ∑ Cn S n + 2 Pn (cos Ψ ). n = nmax +1. 31.
(45) 式中 Pn ：n 階勒建德多項式 Cn ：異常階數變方（anomaly degree variance）. δ Cn ：誤差異常階數變方（error anomaly degree variance） （3-19）式即為（3-17）式中之第一及第三項，亦為（3-15）式之 Ctt 項。 由於球諧係數展開僅到 360 階數，為了模式化 360 階數以上之球諧係數，. Tscherning/Rapp（1974）發展了 T/R 模型，並定義了 Cn 、 δ Cn 項( Basic′ and Rapp, 1992)，其中 Cn =. A(n − 1) (n − 2)(n + B). （3-20）. 式中 A = 425.28mgal 2 , B = 24. δ Cn = (. n GM 2 a ) (n − 1) 2 ( )2( n + 2) ∑ (δ Cn2,m + δ Sn2,m ) 2 m =0 RB RB. （3-21）. 式中. δ Cn2,m、δ Sn2,m：Cn,m、Sn,m 的標準偏差 RB ：Bjerhamar’s 球體半徑 因此，即推導得本文所需使用之觀測重力異常之殘餘量。 再經考慮本文採用交通大學衛星測高化算之水準梯度資料，與陸測、船測為 重力異常資料，將最小二乘配置法公式由（3-15）式改寫為（Hwang, 1997）. N res = (C NΔg. ⎛ CΔg + DΔg C Ne )⎜⎜ ⎝ CeΔg. −1. CΔge ⎞ ⎛ Δg res ⎞ ⎟ ⎜ ⎟ Ce + De ⎟⎠ ⎜⎝ eres ⎟⎠. （3-22）. 式中. C NΔg ：大地起伏、重力異常的協變方矩陣，意義近似（3-15）式之 Cst. 32.
(46) C Ne ：大地起伏、水準梯度的協變方矩陣，意義近似（3-15）式之 Cst C Δg 、 DΔg ：重力異常的協變方矩陣（意義近似（3-15）式之 Ctt ）及隨機誤差. 變方對角矩陣（意義近似（3-15）式之 Cnn ） C Δge 、 C eΔg ：重力異常、水準梯度的協變方矩陣（意義近似（3-15）式之 Ctt ，. 不同處僅為此項有兩類意義不同之觀測量） Ce 、 De ：水準梯度的協變方矩陣（意義近似（3-15）式之 Ctt ）及隨機誤差. 變方對角矩陣（意義近似（3-15）式之 Cnn ） Δg res = Δg − ΔgGGM − Δg RTM ：殘餘重力異常（意義近似（3-15）式之 l ） eres = e − eGGM − eRTM ：殘餘海洋表面梯度（Sea Surface Gradient, SSG） ， （意義. 近似（3-15）式之 l ）. e ：水準梯度. 3-6. 擬大地水準面（Quasigeoid）. Molodensky 提出擬大地水準面觀念，如圖 3-4 將 P 點依據 Helmert 定理投影. 到橢球，當 Q 點的正常位與 P 點的真實位相等時，P 與 Q 就會都投影於橢球上的 同一點，將 Q 點形成的面稱為 telluroid。Q Qo 間距離稱為正常高 H * ， H * 與大地位 數（geopotential number）的關係為線性式。 PQ 間距離稱為高程異常 ζ. ζ =. T. （3-27）. γ. 式中 T ：地面上的擾動位. γ ：telluroid 上的正常重力，由橢球上正常重力向上作自由空間改正獲得. 33.
(47) 將每一點的高程異常 ζ 值加於橢球面之外，所得之面 Molodensky 稱為擬大地 水準面，只是此面並非水準面，亦無物理意義。此面在海面上與大地水準面是一 致的，在陸地上則否，然而擬大地水準面與大地水準面仍舊是非常接近的 （Heiskanen and Moritz, 1979） 。. ζ. H* Qo. 圖 3-4、telluroid 面（Heiskanen and Moritz, 1979） 重力資料不論是在空中或地面觀測，均是在地球表面之外，然而大地起伏的 計算所需的重力或重力異常值卻必須是在大地水準面或擬大地水準面上，因此必 須將重力約化，如 Δg = Δg fa − Δg ref − Δg te. （3-23）. 式中 Δg fa ：自由空間重力異常 Δg ref ：參考之長波長重力異常. Δgte ：地形效應. 重力約化的方法主要有：布格（Bouguer） 、均衡（isostatic） 、Rudzki、Helmert 壓縮(consendation)、剩餘地形模型。由於一合適的約化方式需符合，1、使重力異 常小且值很平滑，2、約化的方法須符合一具有地球物理意義的模型，3、間接效 應應很小（Heiskanen and Moritz,1993） ，其中只有均衡、Helmert 壓縮、RTM 方法. 34.
(48) 能滿足上述要件。 本文使用 RTM 方法，其參考於 Molodensky 定義之擬大地水準面，因此需採 用擬大地水準面與大地水準面間之關係式進行改正，此法稱為 quasigeoid minus geoid separation。依 Bruns 公式得高程異常. ζ =. T (φ , λ , H ) γ (φ , H ). （3-28）. 式中 T ：擾動位. γ ：正常重力 φ , λ ：地理緯度、地理經度 H ：正高. 真實的大地起伏，依 Bruns 近似公式可得 N=. T (φ , λ , 0) γ (φ , 0). （3-29）. ζ 與 N 關係為（Omang and Forsberg, 2000）. ζ −N ≈−. Δg B. γ. （3-30）. H. 式中. Δg B ：布格重力異常， 可得 N ≈ζ −. 2π G ρ. γ. H2. （3-32）. 式中 2π G ρ ：布格項，當地球質量密度為 2.67 g / cm3 時，其值為 0.1119 mgal / m. 由本章 3-2~3-5 節理論計算得之區域重力擬大地水準面，經此參考面的改正， 可獲得一區域重力大地起伏模型。. 35.
(49) 3-7. 台灣區域性大地起伏模型 經由 RCR 方法獲得的為一區域性重力大地起伏模型，此模型參考於全球基. 準；但是一般民生使用的正高所參考的面為區域垂直基準面，此面一般以長年觀 測之驗潮站作為起算點。 幾何法大地起伏值在每一點的成果精度都很高，然而其需要進行水準測量才 能獲得，而水準測量所需時間、人力、金錢成本高，尤其在山區進行精密水準實 行困難度高且精度低，因此不太可能取代重力法大地起伏計算（Kuroishi et al., 2002） ；另外幾何高差法無法考量地球真實密度等物理性質，因此無法預測大地起. 伏變化趨勢。. 重力法大地起伏值採用重力的變化，實質上包括了量測位置、地球密度即地 球內部質量分佈的資訊，且重力法區域大地起伏計算成果，能獲得高解析度及優 良的區域精度，但是卻可能受制於球諧係數的誤差累積，而產生之長波長大地起 伏的系統誤差（Smith and Milbert, 1999） 。 因此重力法及幾何法計算成果具有互補作用（Smith and Milbert, 1999） ，結合 GPS/水準資料獲得之幾何大地起伏值與區域性重力大地起伏模型，所獲得之大地. 起伏模型，實質上即為 GPS 測量之橢球高、及水準測量之正高之合適的轉換中介 模型（Kuroishi et al., 2002） ，如美國 GEOID99 模型，結合美國 G99SSS 重力大地 起伏資訊及 GPS 與水準資訊，就是為了做為美國 NAD83 橢球高與 NAVD88 Helmert 正高之直接轉換面（Smith and Roman, 2000）。 為此，可採 GPS/水準資料控制網中修正重力大地起伏模型，當然第一要件為 確保水準及 GPS 高程儘可能要無誤差，以免造成重力大地起伏模型的誤差。GPS/ 水準獲得的大地起伏資訊 N GPS = hGPS − H levelling. （3-33）. 36.
(50) 兩高程系統差為. ε = N GPS - N gravimetric. （3-34）. 模型化此差值訊號 ε ，再加回區域性重力大地起伏模型，就可把區域性重力 大地起伏模型調整到 GPS/水準基準中，以獲得台灣區域性大地起伏模型。. 37.
(51)
數據
+7
Outline
相關文件
(a) The magnitude of the gravitational force exerted by the planet on an object of mass m at its surface is given by F = GmM / R 2 , where M is the mass of the planet and R is
According to a team at Baycrest’s Rotman Research Institute in Canada, there is a clear link between bilingualism and a delayed onset of the symptoms of Alzheimer ’s and other
Elsewhere the difference between and this plain wave is, in virtue of equation (A13), of order of .Generally the best choice for x 1 ,x 2 are the points where V(x) has
The empirical results indicate that there are four results of causality relationship between Investor Sentiment and Stock Returns, such as (1) Investor
知名的台南花園夜市，歷史雖不如其他台灣 夜市悠久，但在短時間內卻發展成近400個
The angle descriptor is proposed as the exterior feature of 3D model by using the angle information on the surface of the 3D model.. First, a 3D model is represented
運用 Zuvio IRS 與台日比較文化觀點於日本文化相關課程之教學研究 Applying Zuvio IRS and Perspective on Cultural comparison between Taiwan and Japan to Teaching
Meanwhile, the customer satisfaction index (SII and DDI) that were developed by Kuo (2004) are used to provide enterprises with valuable information for making decisions regarding
