陸地上大地起伏模型的計算主要有三種方法：重力法（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
及 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）。
圖3-1、重力法大地起伏計算理論發展
3-2 去除-計算-回復法（remove-compute-restore, RCR）
目前全球最廣為使用於決定區域重力法大地起伏模型的方法，即為去除計算 回復法（Featherstone et al., 2004），即利用觀測之重力異常，去除利用全球大地位 模式理論所計算得之長波長重力異常，再去除利用剩餘地形效應理論所計算得之 短波長重力異常，其值稱為殘餘重力異常（中波長重力異常）。將殘餘重力異常經 最小二乘配置法計算後，獲得殘餘大地起伏值（中波長大地起伏值），最後再回復 長波長大地起伏及短波長大地起伏，若以擾動重力位表示為
GGM RTM res
T =T +T +T （3-1）
式中
TGGM：全球大地位模式的擾動重力位
TRTM：剩餘地形模型理論的擾動重力位
T ：殘餘的擾動重力位，即未模型化的殘餘重力場 res
本文使用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 節），
即可得台灣區域大地起伏模型。
圖3-2、本文大地起伏模型計算流程
3-3 全球大地位模式（Global Geopotential Model, GGM）
S v
, , ,
利用Bruns 公式（Heiskanen and Moritz, 1985），可得參考之長波長大地起伏 值為
3-4 剩餘地形模型理論（Residual Terrain Model, RTM）
由於陸測、船測、衛星測高重力分佈約2′～3′ 一點，為獲得高於重力資料解 析度的區域大地起伏模型，我們就需要結合高程資料（Forsberg, 1984）。
大地起伏的決定當中，針對地形效應的處理主要有三種計算方法，一、地形 均衡改正，其利用觀測得之重力場參數描述地球深處質量如何均衡補償之主要趨 勢，但實際上，地球上很多區域，尤其在深海溝與中海脊處，在深層並不是一定 處於均衡狀態（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、剩餘地形模型
RTM 之大地起伏與重力異常計算公式如下：
(
,)
2 1( ) (
,)
3-5 最小二乘配置法（Least Square Collocation, LSC）
最小二乘配置法可以結合多種資料來源，為此法優點。本文中重力資料的使
S：預估值（於此為大地起伏值）
gL
式中
P ：n 階勒建德多項式 n
C ：異常階數變方（anomaly degree variance） n
Cn
δ ：誤差異常階數變方（error anomaly degree variance）
（3-19）式即為（3-17）式中之第一及第三項，亦為（3-15）式之Ctt項。
C ：大地起伏、水準梯度的協變方矩陣，意義近似（3-15）式之Ne C st
3-6 擬大地水準面（Quasigeoid）
Molodensky 提出擬大地水準面觀念，如圖 3-4 將 P 點依據 Helmert 定理投影
將每一點的高程異常ζ 值加於橢球面之外，所得之面 Molodensky 稱為擬大地 水準面，只是此面並非水準面，亦無物理意義。此面在海面上與大地水準面是一 致 的 ， 在 陸 地 上 則 否 ， 然 而 擬 大 地 水 準 面 與 大 地 水 準 面 仍 舊 是 非 常 接 近 的
（Heiskanen and Moritz, 1979）。
ζ
H*
Qo
圖3-4、telluroid 面（Heiskanen and Moritz, 1979）
重力資料不論是在空中或地面觀測，均是在地球表面之外，然而大地起伏的 應應很小（Heiskanen and Moritz,1993），其中只有均衡、Helmert 壓縮、RTM 方法
能滿足上述要件。
本文使用RTM 方法，其參考於 Molodensky 定義之擬大地水準面，因此需採 用擬大地水準面與大地水準面間之關係式進行改正，此法稱為 quasigeoid minus geoid separation。依 Bruns 公式得高程異常
( , , )
3-7 台灣區域性大地起伏模型 行困難度高且精度低，因此不太可能取代重力法大地起伏計算（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/
水準獲得的大地起伏資訊
GPS GPS levelling
N =h −H （3-33）
兩高程系統差為
GPS gravimetric
= N -N
ε （3-34）
模型化此差值訊號 ε ，再加回區域性重力大地起伏模型，就可把區域性重力 大地起伏模型調整到GPS/水準基準中，以獲得台灣區域性大地起伏模型。