重力資料及計算大地起伏之原理

3-1 重力資料

台灣地區陸上重力資料、船測重力資料與衛星測高資料間會有系統誤差和座 標系統不一致的問題存在，因此在使用前必須經過處理改正將所有資料化算至同 一座標系統，並將系統誤差減至最小。衛星測高資料取自交通大學，點位分佈如 圖 3-1。本文共蒐集六組台灣地區重力資料，各組重力資料詳細內容見黃（2003），

重力資料處理詳細內容見郭（1998）、黃和王（2002）。

第一組重力資料是中央研究院地球科學研究所為研究台灣地區地殼及上部地 涵構造所測得之重力資料（Yen et al., 1990; Yen et al., 1995）。施測點位包含既有 的水準點、三角點及五千或一萬分之一航照圖的標高點，共 603 個分佈均勻的重 力觀測點，其中 308 個測點分佈於資料較為缺乏山區。此組資料內容包含施測點 位之經度、緯度及自由空間重力異常值，點位分佈如圖 3-2。此組重力點之坐標系 統為 TWD67，其坐標及高程取自內政部公佈之三角點及水準點高程。若點位無平 面坐標者，則自地形圖推算之。

第二組重力資料是中國測量學會以 LCR- G 型重力儀於 1986 年 10 月至 1988 年 6 月間測得（黃金維等, 1998），資料內容包含點位之經度、緯度、高程、地表 重力值和自由空間重力異常值，共有 276 點，其中有兩點資料或有平面坐標位置 之錯誤，因此將其刪除而剩下 274 點，圖 3-3 即為此組重力點位分佈圖。

第三組重力資料為中國測量學會於 1997 年至 1999 年間，在台灣地區沿水準 路線施測水準點重力而得，一共有 747 點。此組資料內容包含施測點位之經度、

緯度、正高、地表重力值和自由空間重力異常值，圖 3-4 即為此組重力點位分佈圖。

第四組重力資料為內政部在一等一級水準點上實施重力測量而得，共含 6 個 絕對重力點、10 個一等重力點及 1010 個一等一級水準點之重力值。施測時使用 LCR- G 型重力儀，點位重力值平均標準偏差為 0.07mgal。圖 3-5 為本組重力資料 點位分佈圖。詳細資訊見黃（2001）。

第五組重力資料為內政部提供在在一等二級水準點上實施重力測量而得，施 測時使用 LCR- G 型重力儀，經篩選後剩 1092 個點(分佈見圖 3-6)，新設點位水準 點上的重力值之平均精度為 0.071mgal，最大之重力標準偏差則為 0.14mgal。詳細 資訊見陳(2003)。

第六組重力組資料為台灣環海船測重力資料。本文分析的船測重力資料主要 是取自英國牛津大學之重力資料庫，而較靠近岸邊之船測重力資料主要取自 Hsu et al.（1998）。圖 3-7 之海上位置資料即為此組重力點位分佈圖，本組資料內容包含 所有點位之經度、緯度、重力異常及施測時間。

圖 3-1 衛星測高資料分佈圖

圖 3-2 第一組重力資料分佈圖，係中研院於 1980 至 1987 年施測之 603 個重力觀 測點

圖 3-3 第二組重力資料分佈圖，係中國測量學會於 1986 至 1988 年施測之 276 個 重力觀測點

圖 3-4 第三組重力資料分佈圖，中國測量學會於 1997 至 1999 年施測之 747 個重 力觀測點

圖 3-5 第四組重力資料分佈圖，內政部於 2000 至 2001 年施測之 1010 個重力觀測 點

圖 3-6 第五組重力資料分佈圖，內政部於 2002 至 2003 年施測之一等二級重力觀 測點（篩選後 1092 點）

圖 3-7 第六組重力資料船測重力點及所有陸上重力點位(第一至第五組資料)分佈圖

3-2 去除-回復技術（remove-restore technique）

在大地起伏計算的方法中有一種去除-回復技術（remove-restore technique），

此類型的方法結合積分法、球諧展開式與配置法的優點，克服各自的缺點，如 Forsberg et al.（1981）利用位模型計算長波長分量，利用一個平均高程面的剩餘地 形模型計算地形效應的短波長分量，再以 Stokes 公式對兩者的殘差計算中波長分 量。

重力異常值 g∆ （觀測值）可被分解為三個部分分別計算，即由全球擾動位模 式計算得到的長波長重力異常值，地形效應（短波長）重力異常值與殘餘重力（中 波長）異常值。大地梯度亦同樣可分解為相對應的三個部分（Hwang, 1997；Hwang, 2003）。首先計算殘餘重力異常值與殘餘大地梯度（由衛星測高資料計算而得）：

RTM EGM

res N E

N

N = + ₉₆ + （3-3）

其中N_res、N_EGM96、N_RTM分別為殘餘大地（中波長）起伏、長波長大地起伏與地 形效應（短波長）之大地起伏。圖 3-8 為大地起伏計算流程圖。

圖 3-8 本文大地起伏計算流程圖

3-3 大地位模式

=0

∆V （3-5）

此即為拉普拉斯方程式（Laplace’s Equation），它的解為諧函數。因此地球外部引 力位為諧函數，內部則符合 Poisson 方程式。

若選取地球質心為座標原點，主慣性軸為座標軸，則（3-5）式之解（外部引 力位V ）為（Heiskanen and Moritz, 1985）：

( ) ( ) ^{( )}

Cnm、Snm、C′n0：分別為完全正規化之地球引力位球諧係數與完全正規化之 有一階球諧函數的各項。（Heiskanen and Moritz, 1985）

設地球重力位為W，參考橢球之正常重力位為U，擾動位之定義為同一點上

nm

EGM96 等，其分辨率相當於解算至半波長

美國國防製圖局（NIMA, the National Imagery and Mapping Agency）、美國太 空總署（GSFC, the NASA Goddard Space Flight Center）和俄亥俄州立大學（OSU, the Ohio State University）共同合作發展的一套地球重力位至 360 階的球諧係數模 型，這個新的球諧係數模型稱為 Earth Gravitation Model 1996（EGM96）。EGM96

（Lemoine et al., 1998）包含地形改正資料及 TOPEX/POSEIDON（T/P）、ERS-1 與 GEOSAT 的衛星測高資料，還有大量的衛星追蹤資料包含 SLR（Satellite laser ranging）、GPS、TDRSS（NASA’s Tracking and Data Relay Satellite System）、法國 的 DORIS 和美國海軍的 TRANET 都卜勒追蹤系統。EGM96 最後的解包含三部分： 10∼50 倍。GGM01S 可提供至 120 階的解，GGM01S 結合衛星軌道追蹤資料、地 球表面重力資料與海面測高資料解算得之 GGM01C 則可提供至 200 階的解。最新 的 GGM02S 單純使用 GRACE 任務的飛行資料可提供至 160 階的解，結合地面重 力資料與平均海面高程資料之 GGM02C 可提供至 200 階的解。由於現階段僅提供 至 200 階的解，故本文並未使用此一模式。

3-4 剩餘地形模型理論（Residual Terrain Model, RTM）

計算短波長地形質量效應之大地起伏與重力異常，最常用的有剩餘地形模型

（Residual Terrain Model, RTM）理論與 Helmert 第二壓縮理論（Helmert Second Method of Condensation）等。為配合第二章製作之 DEM 本文採用剩餘地形模型理 論計算短波長之大地起伏與重力異常，RTM 理論在考慮地形質量化算的效應時，

僅考慮地表至參考地形面數間之高頻地形質量效應，如圖 3-9，計算出大地起伏、

重力異常與大地梯度之 RTM 效應後，在於原始觀測資料中扣除 RTM 效應之重力 異常與大地梯度。本文使用之實際地形面為第二章所製作之 3 秒網格解析度之 DEM，參考地形面則為利用 GMT 軟體將 3 秒網格解析度之 DEM 以”filter”指令 製作成 1 分網格解析度之 DEM。

圖 3-9 剩餘地形模型改正圖，斜線部分即為剩餘地形

本文計算剩餘地形效應理論對大地起伏、重力異常、垂線偏差之效應所採用

（1） 大地起伏

(

) (

) (

)

大地座標間的轉換亦需垂線偏差（Heiskanen and Moritz, 1985）。

3-5 最小二乘配置法（least squares collocation, LSC）

最小二乘配置法起源於根據已知點位重力異常來推估（內插和外推）未知點 位重力異常的最小二乘推值法。Krarup（1969）將此法發展為用於不同類型的觀測

（如重力異常和垂線偏差）來推估擾動重力場中的元素，如大地起伏和擾動位等。

後來 Moritz（1980）提出了帶有非隨機參數的最小二乘配置法。最小二乘配置法 可結合各類型的大地測量數據以推算地球形狀及其外部的重力場，但實際運用 時，參與解算之數據的數量受到矩陣大小的限制，因為大矩陣的求逆是非常困難 的，因此必須選擇有代表性的觀測數據，以減少計算量。

根據 Moritz（1980）依帶有雜訊之 LSC 方法，則大地起伏可以表示為：

(

)

C

N = _{sl ll} + _nn ⁻¹ （3-21）

其中

s：訊號（預估量），此為大地起伏

l：觀測量，此為重力異常值與垂線偏差分量ξ 和η C ：訊號與觀測量之協變方矩陣（covariance matrix） sl

C ：觀測量與觀測量之協變方矩陣 ll

C ：隨機雜訊之協變方矩陣，為一對角線矩陣，其對角線元素為誤差變方 nn

最小二乘配置法應用的關鍵是求訊號和觀測量的協變方，由於重力異常值、

垂線偏差、大地起伏與擾動位等都是相互聯繫的，其中最基本的是擾動位，因此

只要知道擾動位的協變方，便可以據此推出其他量的協變方。相關理論見 Moritz

n 係推導而得，詳見 Moritz（1980）、Hwang and Parsons（1995）以及 Tscherning and Rapp（1974）。

考慮本文使用之資料，本文所採用的最小二乘配置法公式為（Hwang, 1997）：

( )

e：為水準梯度

使用最小二乘配置法必須採用去除-回復技術，因為最小二乘配置法理論上要

3-6 週邊填零測試（Zero-padding）

（3-14）、（3-16）與（3-20）式僅為 RTM 效應的第一階近似公式，實際計算 時本文採用 Forsberg（1984）所發展的程式”tcfour”，此程式計算 衝區（buffer zone）填入零值（Press et al., 1992）。本文在 x、y 方向皆填入 50％的 零值。表 3-1 與表 3-2 分別為不同解析度之 RTM 效應之重力異常值與大地起伏填 零與未填零之差值統計，圖 3-9 與圖 3-10 為30′′×30′′解析度之 RTM 效應之重力異 常值與大地起伏填零與未填零之差值等值圖，由圖中可明顯看出差值皆分佈於網 格邊緣，且有等值線穿越台灣北部地區濱海地區，再對照表 4-1、表 5-1 與表 5-2 北部濱海路線之大地起伏值統計結果，發現該測試路線之標準偏差值約 5∼6cm，

表 3-1 RTM 效應重力異常值之填零差異統計表（單位：mgal）

3

3′′× ′′ 30′′×30′′ 2′×2′

mean 2.81897 2.0919 2.53298

std 72.4275 21.9611 23.8853

rms 72.4824 22.0605 24.0181

max 6522.4 632.317 150.684

Min -1735.28 -169.165 -70.5844

表 3-2 RTM 效應大地起伏之填零差異統計表（單位：公尺）

3

3′′× ′′ 30′′×30′′ 2′×2′

mean 0.0196391 0.038494 0.206585 std 0.0604658 0.112731 0.79977 rms 0.0635752 0.119121 0.825984 max 0.571282 0.899476 4.03134 Min -0.361801 -0.339807 -1.53089

圖 3-10 RTM 效應之重力異常值填零與未填零之差值等值圖，等值線間隔 5mgal

圖 3-11 RTM 效應之大地起伏填零與未填零之差值等值圖，等值線間隔 0.05m