國 立 交 通 大 學
土 木 工 程 學 系
博 士 論 文
臺灣地區重力變化：觀測及分析
Gravity Changes in Taiwan：
Observations and Analysis
研 究 生：程自強
指導教授：黃金維
II
臺灣地區重力變化：觀測及分析
Gravity Changes in Taiwan：
Observations and Analysis
研 究 生：程自強 Student： Tze-Chiang Cheng
指導教授：黃金維 Advisor：Dr. Cheinway Hwang
國立交通大學
土木工程學系
博士論文
A Dissertation
Submitted to Department of Civil Engineering College of Engineering
National Chiao Tung University in partial Fulfillment to the Requirements
for the Degree of Doctor In
Civil Engineering July 2013
Hsinchu, Taiwan, Republic of China
臺灣地區重力變化：觀測及分析
研究生：程自強 指導教授：黃金維
國立交通大學土木工程學系
摘 要
本論文主要以臺灣地區各式地面重力觀測成果，檢視並分析各地區重力變化 情形，配合高程變化等相關資訊對重力變化機制作初步分析。2003 年及 2007 年， 共進行了兩次的全臺一等水準點上重力測量，所有重力觀測量皆須經過環境改正 後再行應用或做後續計算，其中本文特別針對海潮負載效應改正作較深入的探討 及分析。海潮負載效應對地處西太平洋的臺灣而言，重力觀測量的影響極為顯著， 本文即以適用臺灣地區的 OTL10 區域海潮模型，配合 SGOTL 程式計算求得海 潮負載效應，對臺灣各地的重力觀測量作改正，並同時以 g7 及 SPOTL 程式計算 所得的海潮負載效應對絕對重力觀測量做改正，其中以 SGOTL 所改正後的絕對 重力值，可以得到最佳的精度表現。相對重力觀測量的標準偏差表現，主要來自 於觀測時環境干擾的程度，而藉由增加重力的觀測數量，不但可重複檢核觀測量 的閉合差，更可有效提升重力成果的精度。就以 2003 至 2007 年的重力變化情形 來看，臺灣中央山脈北段部分重力值變化量可達-0.085 mgal/year，臺灣西部地區 則達到 0.061 mgal/year，經比較正高變化情形後，可以明顯研判中央山脈部分地 區的造山運動仍持續進行中，而臺灣西部則有地層下陷的情形發生。於 2004 年 至 2008 年期間，本文利用絕對重力及相對重力方法檢視雲林地層下陷情況，以 FG5 絕對重力儀於雲林縣同安國小觀測所得重力值變化率達到 22.72 μgal/year， 同時以相對重力進行六次重力網形觀測，經平差計算後求得雲林地層下陷區的重 力變化值，並與該區域各地的正高變化率作比較，可初步驗證重力方法對於地層 下陷監測的可行性。II
Gravity Changes in Taiwan：
Observations and Analysis
Student : Tze-Chiang Cheng Advisor : Dr. Cheinway Hwang
Department of Civil Engineering
National Chiao Tung University
Abstract
The temporal gravity changes in Taiwan are used to examine the crustal deformation in Taiwan. All gravity observations are corrected by the environmental effects, such as solid earth tide, ocean tide loading, polar motion, atmosphere and underground water.
This study will put a special emphasis on ocean tide loading. Ocean tide loading(OTL)
effect of gravity is large in Taiwan and around its nearby area. Three different OTL models are tested to improve the accuracy of absolute gravity around the coast of Taiwan and its islets. The model SGOTL with the latest regional ocean tide model of Taiwan is assimilated with local tide gauge records. In Taiwan, OTL10 and the high resolution DEM result in the optimum accuracy in the OTL corrections of absolute gravity measurements. It’s obviously to shown the standard deviations of adjusted gravity have reduced by using mixed method of ladder and star. Two major relative network campaigns around the whole Taiwan are made from 2003 to 2007. The relative gravity observations are adjusted using the weight-constraint least-squares method. With the 4-year gravity changes, we have found that orogeny occurs on the central ridges and
land subsidence occurs in southwestern Taiwan with average rates of -0.085 mgal/year
Taiwan is established to determine gravity variations caused by large land subsidence. At station TAES, the gravity change rate from absolute measurements by FG5 is 22.72 μgal/year. The gravity-height admittance factor in Yunlin subsidence area is -5.25 μgal/cm, instead of -1.97 μgal/cm, due to the complicated mechanism of subsidence. This result can be a prototype method of subsidence monitoring by gravimetry.
IV
致 謝
能夠完成這本博士論文，進一步地說能夠完成博士的學程，這些年來要感謝 的人很多，首先要感謝的就是我的指導教授 黃金維老師。黃老師在做學問及研 究的熱忱與執著，一直是我景仰和學習的對象，今生有幸能成為他的門生，更在 他嚴謹地指導下做研究，使我的學識能夠逐漸增長。他除了指導我的學業及研究 之外，也讓我有機會與國際權威級學者專家共同研習，更帶領我前往澳洲及美國 等地參加國際學術發表會議，讓我能放眼世界增廣見聞。由於在我完成學分修習 及通過博士資格考後，即回到工作崗位復職，又因為工作時間的不穩定，導致我 無法安排較完整的時間繼續進行研究，而黃老師在那段期間總是督促我要持續地 做研究，最後因為考量工作與學業無法兼顧，在獲得黃老師全力支持後，毅然離 開工作回到交大，於最後一學年全力衝刺完成博士的學程。 另外要感謝當時工作單位長官的提攜與支持，也感謝工作時的夥伴甯方璽、 姜震昇、高豫麒及劉榮寬四位學長的鼓勵及指點迷津，讓我不後悔地選擇進入交 大；並感謝工研院量測中心的李瓊武博士及文祺，讓我與 FG5 絕對重力儀結下 了不解之緣，並且能以此題材完成研究學業，更要特別感謝鄭景中博士對我在做 學問及研究上的指導與幫助。 在做研究的枯燥乏味生活中，多虧了在研究室一同努力的學長姐及學弟妹們 的鼓勵，而你們努力的精神與態度正是我學習的榜樣，特別是我在最後一學期如 無頭蒼蠅般埋首研究期間，小光(煥欽)的精神鼓勵與長儀的經驗指導，讓身陷研 究泥沼中的我能夠再度邁步前進。 當然，在此衷心感激我的愛妻珮瑜，當我拋妻棄女衝刺學業的這段時間，妳 毫無怨尤地獨撐家計，讓我無後顧之憂地專心完成學業及研究，近十年的博士班 學業，當初兩個小女娃的宇恬及宇柔，如今已亭亭玉立並且聽話懂事，都要感謝 愛妻的辛勤教養，在此我只想說~珮瑜，這輩子有妳真好！最後，我要感謝辛苦 養育我的父母親，我所有努力的成果以及這一切的榮耀，都是屬於你們倆老！目 錄
摘 要... I Abstract ... II 致 謝... IV 目 錄... V 表目錄... VII 圖目錄... VIII 一、緒論... 1 1-1 研究目的 ... 1 1-2 文獻回顧 ... 2 1-3 論文架構 ... 4 二、重力測量及平差方法 ... 6 2-1 絕對重力測量 ... 6 2-1-1 FG5 絕對重力儀 ... 7 2-1-2 重力梯度及重力化算 ... 9 2-1-3 絕對重力值計算 ... 10 2-2 相對重力測量 ... 132-2-1 Lacoste & Romberg G 型重力儀 ... 15
2-2-2 Lacoste & Romberg Graviton-EG 型重力儀 ... 17
2-2-3 Scintrex CG-5 型重力儀 ... 21
2-3 重力環境改正 ... 23
2-4 重力觀測量平差計算 ... 24
2-4-1 自由基準平差（datum-free adjustment） ... 25
2-4-2 加權約制平差（weighted constraint adjustment） ... 27
2-4-3 顯著性測試 ... 29 三、海潮負載改正 ... 31 3-1 概述 ... 31 3-2 牛頓引力效應 ... 32 3-3 彈性效應 ... 34 3-4 海潮負載格林函數 ... 36 3-5 海潮負載振幅及相位 ... 38 3-6 海潮模型 ... 40 3-7 三種海潮負載模式於近岸及離島的比較與重力值改正 ... 41
VI 3-7-1 海潮負載重力效應計算軟體 ... 42 3-7-2 絕對重力值改正成果比較分析 ... 44 3-8 全臺相對重力觀測量海潮負載效應改正 ... 51 3-9 雲林地區相對重力觀測量海潮負載效應改正 ... 58 四、全臺 2003 年及 2007 年重力變化及分析 ... 60 4-1 概述 ... 60 4-2 觀測資料 ... 60 4-3 全臺相對重力網平差成果 ... 62 4-3-1 自由基準平差 ... 62 4-3-2 加權約制平差 ... 63 4-4 重力變化分析 ... 66 4-4-1 重力變化及其原因 ... 66 4-4-2 重力值精度分析 ... 78 4-5 本章小結 ... 83 五、以重力方法檢視地層下陷 ... 84 5-1 概述 ... 84 5-2 重力檢視地層下陷基本原理 ... 86 5-3 絕對重力測量 ... 87 5-4 相對重力測量 ... 88 5-5 重力變化分析 ... 90 5-6 本章小結 ... 101 六、結論與建議 ... 103 參考文獻... 105 附錄 A FG5 絕對重力儀操作說明 ... 111 附錄 B Scintrex CG-5 型重力儀操作說明 ... 121 作者簡歷... 131 學術著作... 132
表目錄
表 2-1：絕對重力值觀測成果 ... 11 表 2-2：本文臺灣西部中海拔山麓地區的相對重力觀測量 ... 20 表 2-3：本文各式相對重力儀的諸元比較表 ... 22 表 3-1：Farrell 的各階洛夫常數... 36 表 3-2：G-B 地球彈性體模型海潮負載格林函數表 ... 37 表 3-3：NAO 海潮模型... 40 表 3-4：本文綠島絕對重力觀測量使用三組海潮負載效應改正後之成果 ... 47 表 3-5：本文東河國小絕對重力觀測量使用三組海潮負載效應改正後之成果 ... 48 表 3-6：本文澎湖絕對重力觀測量使用三組海潮負載效應改正後之成果 ... 49 表 3-7：本文墾丁絕對重力觀測量使用三組海潮負載效應改正後之成果 ... 50 表 3-8：本文 2003 年全臺重力網自由基準平差重力值成果概述 ... 53 表 3-9：本文 2007 年全臺重力網自由基準平差重力值成果概述 ... 55 表 3-10：本文雲林重力網自由基準平差重力值各時段成果概述 ... 59 表 4-1：本文共同觀測之絕對重力點歷次紀錄 ... 61 表 4-2：本文 2003 年及 2007 年全臺相對重力自由基準平差成果統計概述 ... 62 表 4-3：本文 2003 年全臺相對重力加權約制兩處平差成果統計概述 ... 63 表 4-4：本文 2007 年全臺相對重力加權約制兩處平差成果統計概述 ... 64 表 4-5：本文 2003 年全臺相對重力加權約制三處平差成果統計概述 ... 64 表 4-6：本文 2007 年全臺相對重力加權約制三處平差成果統計概述 ... 64 表 4-7：本文全臺相對重力加權約制四處平差成果統計概述 ... 65 表 4-8：2000 至 2008 年全臺一等水準測量作業內容略表 ... 65 表 4-9：本文 Scintrex CG-5 重力儀於各地質層重力觀測量標準偏差表現情形 .. 79 表 5-1：雲林地區 1992 年至 2002 年下陷資料統計表 ... 85 表 5-2：本文各次相對重力網形平差重力值精度表現情形 ... 94 表 5-3：地表沉積層各種物質孔隙率分布略表 ... 97 表 5-4：本文各重力觀測位置地下水位變化率及其改正後重力變化率 ... 98 表 5-5：本文以重力方法及水準測量推得之地層變化率 ... 101VIII
圖目錄
圖 2-1：自由落體計算重力加速度 ... 6 圖 2-2：FG5 絕對重力儀觀測現況 ... 8 圖 2-3：簡易重力梯度值測量方法 ... 10 圖 2-4：嚴密重力梯度值測量方法 ... 10 圖 2-5：相對重力觀測進行方式 ... 13 圖 2-6：全臺重力觀測網形圖 ... 14圖 2-7：Lacoste & Romberg G 型重力儀 ... 15
圖 2-8：Lacoste & Romberg G 型重力儀內部組件圖 ... 16
圖 2-9：Lacoste & Romberg Graviton-EG 型重力儀 ... 17
圖 2-10：Lacoste & Romberg Graviton-EG 型重力儀內部組件配置圖 ... 18
圖 2-11：地震震波傳遞波形 ... 19 圖 2-12：本文臺灣西部沿海地區 Gravition-EG 型重力儀相對重力觀測量 RMSError 分布直方圖(2004/12) ... 19 圖 2-13：本文日本京都地區 Gravition-EG 型重力儀重力觀測量 RMSError 分布直 方圖(2004/11/23) ... 21 圖 2-14：Scintrex CG-5 重力儀 ... 22 圖 2-15： t 分布之機率密度函數 ... 30 圖 3-1：本文中海潮負載效應於觀測位置 p 與水體 q 之間關係 ... 32 圖 3-2：本文中臺灣離島及本島近岸絕對重力站觀測位置 ... 44 圖 3-3：本文中綠島絕對重力站(AG1)海潮負載效應表現 ... 45 圖 3-4：本文中東河國小絕對重力站(AG2a)海潮負載效應表現 ... 45 圖 3-5：本文中澎湖絕對重力站(AG9)海潮負載效應表現 ... 46 圖 3-6：本文中墾丁絕對重力站(KDNG)海潮負載效應表現 ... 46 圖 3-7：本文 2003 年全臺海潮負載改正前與後的重力值標準偏差分布情形 ... 53 圖 3-8：本文 2003 年全臺海潮負載改正後的重力值標準偏差變化情形 ... 54 圖 3-9：本文 2003 年全臺海潮負載改正後的重力值變化情形 ... 55 圖 3-10：本文 2007 年全臺海潮負載改正前與後的重力值標準偏差分布情形 ... 56 圖 3-11：本文 2007 年全臺海潮負載改正後的重力值標準偏差變化情形 ... 57
圖 3-12：本文 2007 年全臺海潮負載改正後的重力值變化情形 ... 57 圖 4-1：本文 2003 年及 2007 年全臺相對重力測量共同觀測絕對重力點分布圖61 圖 4-2：全臺正高變化率及本文重力變化率略圖 ... 66 圖 4-3：本文 2003 至 2007 年蘭陽平原地區重力變化及高程變化初步比較 ... 67 圖 4-4：本文 2003 至 2007 年蘭陽平原地區重力及高程變化初步比較圖 ... 67 圖 4-5：1992 至 2012 年宜蘭地區累積下陷量分布圖 ... 68 圖 4-6：本文 2003 至 2007 年雲嘉平原地區重力變化及高程變化初步比較圖 ... 69 圖 4-7：本文 2003 至 2007 年雲嘉平原地區重力及高程變化初步比較圖 ... 69 圖 4-8：本文 2003 至 2007 年中北部山脈地區省道 7 號、7 甲、8 號、21 號及 18 號公路沿線重力變化及高程變化初步比較... 70 圖 4-9：本文 2003 至 2007 年中北部山脈地區省道 7 號、7 甲、8 號、21 號及 18 號公路沿線一等水準點重力及高程變化初步比較圖... 71 圖 4-10：本文通過顯著性測試之 2003 至 2007 年一等水準點重力變化與相對應位 置正高變化分布略圖... 72 圖 4-11：本文通過顯著性測試之 2003 至 2007 年中央山脈北部地區省道 7 甲及 8 號部分公路沿線一等水準點重力與正高變化情形... 73 圖 4-12：本文通過顯著性測試之 2003 至 2007 年於中央山脈北部地區省道 7 甲及 8 號部分公路沿線一等水準點重力變化與相對應位置正高變化表現略圖... 73 圖 4-13：青藏高原地殼增厚示意圖 ... 75 圖 4-14：本文試以布格平板概念推算台灣中央山脈北部地區省道 7 甲及 8 號部分 公路沿線地殼厚度增厚之可能情形... 75 圖 4-15：台灣地區的地殼厚度分布略圖 ... 76 圖 4-16：本文 2003 至 2007 年雲林地區一等水準點重力與正高變化情形 ... 77 圖 4-17：本文 2003 至 2007 年於雲林地區一等水準點重力變化與相對應正高變化 表現略圖... 77 圖 4-18：本文 2007 年全臺相對重力觀測量(Scintrex CG-5)標準偏差及平差重力值 精度分布情形... 78 圖 4-19：本文重力觀測位置於地質層分布略圖 ... 80 圖 4-20：本文 AGTO 相對重力觀測點位示意圖 ... 81 圖 4-21：本文 AGTO 相對重力觀測量(Scintrex CG-5)標準偏差分布情形 ... 82 圖 4-22：本文南橫公路相對重力測量平差重力值殘差及標準偏差分布情形 ... 82
X 圖 5-1：雲林地區 2002 至 2003 年年平均下陷速率圖 ... 85 圖 5-2：本文中布格平板示意圖 ... 86 圖 5-3：本文雲林相對重力觀測混合進行方式示意圖 ... 89 圖 5-4：本文雲林相對重力觀測基本網形 ... 89 圖 5-5：本文 2004 年 10 月至 2013 年 5 月同安國小絕對重力值變化情形 ... 90 圖 5-6：本文同安國小及北港衛星追蹤站絕對重力值變化情形 ... 91 圖 5-7：本文 EG 與 CG-5 相對重力儀觀測量標準偏差分布圖 ... 92 圖 5-8：本文第一次(2004 年 12 月)相對重力網形殘差與標準偏差分布情形... 92 圖 5-9：本文第二次(2005 年 4 月)相對重力網形殘差與標準偏差分布情形... 92 圖 5-10：本文第三次(2005 年 12 月)相對重力網形殘差與標準偏差分布情形.... 93 圖 5-11：本文第四次(2006 年 7 月)相對重力網形殘差與標準偏差分布情形 ... 93 圖 5-12：本文第五次(2006 年 11 月)相對重力網形殘差與標準偏差分布情形 .... 93 圖 5-13：本文第六次(2007 年 4 月)相對重力網形殘差與標準偏差分布情形... 94 圖 5-14：本文各相對重力觀測位置重力值變化情形 ... 95 圖 5-15：本文濁水溪流域地質分布略圖 ... 97 圖 5-16： 雲林地區水準測量環線分布圖 ... 99 圖 5-17： 本文雲林地層下陷區各重力觀測位置重力與正高變化比較 ... 100 圖 5-18：濁水溪沖積扇之水文地質概念分層圖 ... 102
一、緒論
1-1 研究目的
本文首先針對全臺一等水準點上重力測量成果進行計算。自 2001 年起，由 內政部委由國立交通大學進行一等一級與一等二級水準點上重力測量，即開啟全 面及系統性重力測量工作，於 2003 年 6 月完成全臺一等水準點重力測量(內政部， 2003)，另於 2006 年起，內政部委由寰宇測量公司進行並再次完成全臺一等水準 點重力測量(內政部，2006)。 在所有重力觀測量中，均須經過環境改正，本文特別針對海潮負載效應改正 做計算及分析比較，由於臺灣位於西太平洋邊緣，海潮負載效應所產生的重力值 對重力觀測成果的影響極為顯著，此影響部分需藉由海潮負載理論模式進行觀測 量改正，以期獲得較佳的重力成果。本文引用 SGOTL(Superconducting Gravity for Ocean Tide Loading; Hwang and Huang, 2012)海潮負載改正模式，對臺灣全島相對 重力測量觀測量做改正，並比較海潮負載改正前後之重力值精度變化情形。海潮 負載效應對於離島及近岸地區的重力觀測有明顯的影響量，由於絕對重力是以連 續觀測的方式進行，因此每個觀測時間的海潮負載效應改正量均應列入考量，除 了使用 FG5 絕對重力儀內建軟體 g7 計算海潮負載效應之外，本文另使用 SGOTL 與 SPOTL(Some Programs for Ocean Tide Loading; Agnew, 2013)兩種海潮負載改 正模式，並對臺灣近岸及離島絕對重力觀測量做改正，最後比較三種海潮負載改 正模式於重力計算後成果的差異，選取最適合的改正模式，藉此使各絕對重力值 達到最佳的觀測精度。 內政部於 2003 年底陸續購進 FG5 絕對重力儀及超導重力儀，使我國可自主 觀測臺灣地區絕對重力值及精確訂定重力基準，並於新竹市十八尖山設置國家重 力基準站，為我國重力測量開啟重要里程碑。以 FG5 絕對重力儀所測量的觀測 精度，在連續觀測 12 小時以上均可達到 1 至 2 μgal 等級的精度，而相對重力測2
量的觀測標準偏差多半取決於觀測地點的安靜程度，若增加觀測量至相對重力網 形中，不但可重複檢核觀測量閉合差，更可有效提升重力成果的精度達 5.4μgal 等級。其中本文所分析比較的 Lacoste & Romberg G 型、Lacoste & Romberg Graviton-EG 型及 Scintrex CG-5 型相對重力儀，在觀測量的標準偏差與重力成果 的精度表現上，均屬 Scintrex CG-5 型相對重力儀為最優。
由於臺灣地處歐亞板塊與菲律賓海板塊交界處，地殼變動與地震發生頻繁， 菲律賓海板塊仍以 7 至 8 cm/year 速率向歐亞板塊聚合(Seno, 1977； Yu et al., 1986)，除了地表可觀測得到變化量之外，針對地表下所造成的重力變化，可由 多次的重力觀測值，依照其時變結果，分析及研判臺灣各地的地層變動表現情形。 臺灣地區重力變化表現，經由與正高變化情形比較後，中央山脈部分地區呈現重 力值下降且正高呈現上升的趨勢，而西部平原地區重力值普遍增加且正高值呈現 下降的趨勢，可明顯看出中央山脈隆起現象與西部地層下陷情形。 臺灣西部平原地區多為人口聚居及工商業廠所，另外則是農業種植及灌溉區 域，雲林地區西部因為農業種植區域廣大，境內所需灌溉水源缺乏，轉而以抽取 地下水方式進行灌溉或養殖，多年來已造成嚴重地層下陷，又由於臺灣高速鐵路 路線行經該地層下陷區，因此對於行車安全也是項隱憂。本文利用重力測量方法 檢視雲林地區地層下陷狀況，自 2004 年 12 月至 2013 年 5 月的長期觀測資料中， 雲林縣同安國小(TAES)所觀測得到的絕對重力值累積變化已達到 124.4±2.7 μgal， 此重力變化量所顯示的物理意義，即為是地層下陷以及其他的各種影響量的表 現。
1-2 文獻回顧
臺灣地區重力測量的約制點位來自於絕對重力點，在國內未引進FG5絕對重 力儀之前，絕對重力測量曾經由美國國家大地測量局(National Geodetic Survey) 來臺做觀測，於1991年起，先後利用JILLA-4與FG5絕對重力儀在大溪、新竹、 臺中、鳳山、綠水及太魯閣等六處進行選點及觀測(張瑞剛，1997)，2003年內政部購進FG5絕對重力儀後，即由國立交通大學與工研院量測技術中心合作，展開 對臺灣地區的絕對重力測量，於2005年完成全臺15處絕對重力站設置與測量(內 政部，2005)，並於2006年完成國家重力基準站的建置(內政部，2008)，2006年11 月起，與法國合作進行AGTO(Absolute Gravity of the Taiwanese Orogeny)計畫， 對臺灣中央山脈南部東西向一帶，進行逐年的絕對重力測量，以分析該區域的地 殼變動情形(內政部，2009)，Mouyen et al.( 2012）在研究臺灣南部山區侵蝕作用 時，即結合絕對與相對重力測量進行地殼變動監測；Kao(2011)更對國內絕對重 力與超導重力應用在重力基準、地體動力及環境變遷做相關的研究及探討。 內政部自2001年起即展開全國性的一等水準點上重力測量，並檢測國內各絕 對重力點變動情形，以做為重力網形平差的約制條件，計算出全國各地重力值， 另分析各式平差模式對重力成果精度的影響(Hwang et al., 2002），另於2006年起， 又再次進行全國一等水準點上重力測量作業；在水準測量部分，Yang et al.（2003） 經由分析2001年全國1010個一等水準點測量成果後，重新訂定全新的臺灣垂直基 準系統(TWVD2001)，後續並逐年實施一等水準測量，Chen et al.（2011）研究並 分析2001年至2008年的正高成果後發現，位於臺灣西南部地表的最大下陷率為 -109.4 mm/year；另有Ching et al.（2011）結合GPS與水準資料，並運用各種模式 分析全臺灣地區地殼變動；Sun et al.(2009)於中國青藏高原及雲南等地以重力及 GPS觀測資料研究並驗證板塊增厚。 陳南松(2003)曾探討固體潮及海潮負載效應對臺灣地區衛星追蹤站位置及 重力影響的變化；在本論文中各式重力觀測資料，海潮負載重力效應對觀測量的 影響不可輕忽，黃鉅富(2012)對此研究開發一套適用臺灣地區的海潮模型 OTL10， 並以 SGOTL 程式計算海潮負載效應；國際上另有 Agnew(2013)所開發的 SPOTL 相關程式用以計算海潮負載效應，Zhou et al.(2013) 的研究中，使用 SPOTL 計算 海潮負載效應，進一步分析中國大陸東南沿岸地區重力觀測位置的絕對重力值改 正情形，藉以改善絕對重力觀測精度。
4
使用 GPS 連續觀測或以地陷監測井做為監測方法(經濟部水利署，2004；陳鶴欽， 2008)，洪偉嘉(2009)對於大面積的監測更增加應用 InSAR(Interferometry Synthetic
Aperture Radar)技術，Hung et al.(2011)並以 InSAR 應用於臺灣濁水溪沖積扇的地
表變形監測，Jacob et al.(2010)於 2006 年至 2008 年期間，在法國南部以絕對重 力及相對重力測量來檢視該地區水蝕石灰岩地形地下水位的變化，藉此推測出的 地下水的儲量。Hwang et al.(2010)則以重力測量方法對雲林地層下陷區進行初步 監測實驗，主要是以 FG5 絕對重力儀及相對重力儀觀測下陷區與週邊相關位置 重力變化做比較分析。
1-3 論文架構
本論文主要分成六大章節及相關附錄，其重點內容架構如下： 一、「緒論」：闡明本研究之目的、文獻回顧及論文架構。 二、「重力測量及平差方法」：主要介紹地面重力測量理論及方法，分別闡述絕對 重力測量觀測理論及方法、相對重力測量觀測理論及方法、重力環境改正與重力 觀測量平差計算。 三、「海潮負載改正」：介紹海潮負載的原理及相關理論公式，並利用 SGTOL 程 式計算海潮負載改正，計算並比較相對重力網平差成果精度。另針對近岸與離島 的絕對重力觀測量，以 g7、SGOTL 與 SPOTL 三種程式進行海潮負載效應計算， 並加入觀測量中進行改正，以求出最佳化的絕對重力值及觀測精度，並對各相對 重力網於改正海潮負載效應前後的成果精度做比較。 四、「全臺 2003 年及 2007 年重力變化及分析」：闡述 2003 年與 2007 年兩次於臺 灣全島一等水準點上相對重力觀測與成果，藉由比較正高變化與重力變化做初步 分析，檢視全臺各地的地表變動情形。 五、「以重力方法檢視地層下陷」：介紹以重力測量方法觀測雲林地層下陷情形， 探究其重力變化所代表的物理意義，並分析重力與地層下陷的初步相關性。 六、「結論與建議」：綜整本論文之研究成果及結論，並對研究中不足之處提出說明與建議。
6
二、重力測量及平差方法
2-1 絕對重力測量
絕對重力測量即是於某一點位置量取其重力加速度的大小。本文中所討論的 絕對重力測量原理，主要是以自由落體法則來測定。假設有某一物體由一高度 h 自由落下至h ，則落下距離0 h h 0與落下時間 t 的關係式為： 2 0 0 1 2 hh v t gt 即為 2 0 0 1 2 hh v t gt (2-1) 上式(2-1)中 g 為重力加速度值，h 表示於時間 t=0 物體所在位置，h 則是時間為0 t 時刻物體所在位置，v 表示 t=0 時物體的速度。根據上式，只要量測不同落下0 時間t 及與其對應的落下距離i hih0，則可求解出重力加速度值 g。由於式(2-1) 中有三個未知數(h 、0 v 與 g)，因此必須測定至少三組的0 h 及i t 值以組成聯立方i 程式求解出 g 值。本文中所討論的 FG5 絕對重力儀，即是採用如上觀測原理。 圖 2-1：自由落體計算重力加速度如圖 2-1，自由落體在三個位置上的落下時間及距離分別為( , )t h 、1 1 ( ,t h 及2 2) 3 3 ( ,t h ，再假設 ) 1 1 0 d h h 2 2 0 d h h 3 3 0 d h h 依據(2-1)式，可得： 2 1 0 1 1 1 2 d v t gt 2 2 0 2 2 1 2 d v t gt 2 3 0 3 3 1 2 d v t gt 以上三式中，將第二式及第三式均減去第一式，可得： 2 1 2 1 0 2 1 1 ( ) ( )[ ( )] 2 d d t t v g t t 3 1 3 1 0 3 1 1 ( ) ( )[ ( )] 2 d d t t v g t t 令D1d2d1，D2d3d1，T1 t2 t1，T2 t3 t1，代入上式後，消去 v0，簡 化後則可得 2 1 2 1 2 1 2 ( ) ( ) D D T T g T T (2-2) 上式中的 g 為重力加速度值，亦即是我們要測定的絕對重力值。以 FG5 絕 對重力儀為例，該儀器即是以碘穩頻雷射來測定(2-2)式中的落下距離 D，並以銣 原子鐘測定落下時間 T，觀測精確度可達 10-8 ms-2，亦即μgal 等級(Kao, 2011)。 2-1-1 FG5 絕對重力儀 FG5 絕對重力儀主要部份是一部極精準的碘穩頻雷射干涉儀，用來監測物體
8 自由落下的直角稜角運動情形，干涉儀所收集到的光學資料即提供了非常精確的 距離觀測量系統。簡單來說，其原理就是以精密的測距系統搭配精密計時的原子 鐘，精準地測定物體落下的重力加速度值，經由大量地收集落下的數據，見圖 2-2，最後以電腦軟體程式計算得出一平均加速度值，即為絕對重力值。 圖 2-2：FG5 絕對重力儀觀測現況 國內自 2003 年底即陸續引進兩部 Micro-g Lacoste FG5 絕對重力儀，使我國 重力測量的技術可與世界同步，除可自主測定臺灣地區的絕對重力值之外，再搭 配設置超導重力儀，並於新竹市十八尖山完成建立國家重力基準站(Laboratory of Geodesy and Geodynamics)。
絕對重力的觀測作業方法，主要是架設絕對重力儀於待測位置上，收集該位 置的各次重力加速度值，經由統計分析原理以求定平均重力加速度值，即為所測 定的絕對重力值。然而 FG5 絕對重力儀所測定重力加速度的位置，是在該儀器 上半部的落體裝置，實際與待測地面仍有高度落差，FG5 儀器所訂定的此高度落 差為 100 公分，亦即表示所測得的絕對重力值必須再經由化算至地面後，才是實 際待測位置的絕對重力值。作者於國內剛引進 FG5 絕對重力儀時，即參與內政
部委託工研院量測技術中心接收儀器事宜，並隨該中心工程師於 Micro-g Lacoste 公司技師指導下，學習 FG5 絕對重力儀架設與測量等事務，之後除了於工研院 量測技術中心實驗室內進行 FG5 絕對重力儀的測量工作外，更前往臺灣最南部 及最北部的內政部墾丁(KDNM)及陽明山(YMSM)衛星追蹤站進行絕對重力觀測， 此舉為國內首次自主性的絕對重力觀測實驗任務，其後又於 2005 年完成全臺 15 個絕對重力點的設置與測量(內政部，2005)，另又以絕對重力觀測應用於雲林地 層下陷地區的重力監測(Hwang et al., 2010)，就重力監測地層下陷部分於第五章 做深入討論。本論文並將 FG5 絕對重力儀相關作業模式與操作經驗整理後，收 錄於附錄 A 做詳細說明。 2-1-2 重力梯度及重力化算 重力值所根據的高度化算，首要測定該位置的重力梯度值，以自由空間理論 的重力梯度值而言為 -3.086 μgal/cm，然而考量各觀測重力位置的經、緯度與高 度皆不同，必須於待測位置實際測定其重力梯度值。簡易觀測方法是以相對重力 儀於待測位置上零高度先測定一重力值，再於待測位置上方某單位高度再測定一 重力值，重複上述觀測步驟數次後，取其重力差值的平均值，再除以單位高度後， 即得到該位置的重力梯度值，見圖 2-3。另外較嚴密的重力梯度值測定方法，是 利用某些固定高度的支架，見圖 2-4(Kao, 2011)，依照各種不同的高度 h2、h1， 以相對重力儀測量出上、下高度的相對重力值 g2、g1，重力梯度值則可如下式(2-3) 表示： 2 1 2 1 grad g g g h h (2-3) 重複以上做法，並改以各不同高度組合施測，最後取其平均值，即得到該位置的 重力梯度值。
10 圖 2-3：簡易重力梯度值測量方法(左：零高度測定重力值，右：單位高度測定重 力值) 圖 2-4：嚴密重力梯度值測量方法(Kao, 2011) 2-1-3 絕對重力值計算 絕對重力的觀測可依照需求，藉由 Micro-g Lacoste 公司所開發的 g7 程式軟 體(g7 User’s Manual, 2006)，設定欲進行的觀測組數(set)及各組內的落下數(drop) 與間隔時間，控制 FG5 絕對重力儀進行觀測作業，完成觀測後即可得到平均之 絕對重力值，而其重力值的觀測精度來自於總落下數以及儀器本身的系統誤差 (不確定度)。以多組數量連續觀測同一位置的平均重力值標準偏差如下：
2 1 ( ) 1 n i i g g n
(2-4) 上式(2-4)中， 為組離散度(Set scatter)，n
為重觀測組數量(Set)，gi為組重力平均值，g為平均重力值，平均重力值的精度即為： g n (2-5) FG5 在完成絕對重力觀測作業後，會由軟體 g7 計算並產生觀測成果記錄檔，觀 測成果紀錄檔範例見表 2-1。 表 2-1：絕對重力值觀測成果(2005 年同安國小觀測範例) Micro-g Solutions g Processing Report
File Created: 12/22/05, 08:23:43 Project Name: 1214C
g Acquisition Version: 3.1021 g Processing Version: 4.0405 Company/Institution:
Operator: Cheng, Tze-Chiang Station Data
Name: Tong Ann Elementary School Site Code: TAES
Lat: 23.71164 Long: 120.26461 Elev: 5.00 m Reference Height: 13.10 cm
Datum Height: 100.00 cm Gradient: -3.900 uGal/cm
Nominal Air Pressure: 1012.65 mBar Barometric Admittance Factor: 0.30 Polar Motion Coord: 0.0633 " 0.3874 " Earth Tide (ETGTAB) Selected
Potential Filename: C:\Program Files\Micro-g Solutions Inc\gWavefiles\ETCPOT.DAT Delta Factor Filename: C:\gData\OceanLoad.dff
Delta Factors
Start Stop Amplitude Phase Term 0.000000 0.002427 1.000000 0.0000 DC 0.002428 0.249951 1.160000 0.0000 Long 0.721500 0.906315 1.154250 0.0000 Q1 0.921941 0.974188 1.154240 0.0000 O1 0.989049 0.998028 1.149150 0.0000 P1 0.999853 1.216397 1.134890 0.0000 K1 1.719381 1.906462 1.161720 0.0000 N2 1.923766 1.976926 1.161720 0.0000 M2 1.991787 2.002885 1.161720 0.0000 S2 2.003032 2.182843 1.161720 0.0000 K2 2.753244 3.081254 1.07338 0.0000 M3 3.791964 3.937897 1.03900 0.0000 M4 Ocean Load ON, Filename: C:\gData\OceanLoad.olf
Waves: M2 S2 K1 O1 N2 P1 K2 Q1 Mf Mm Ssa Amplitude (uGal): 2.000 0.650 2.352 2.088 0.683 0.778 0.251 0.468 0.000 0.000 0.000 Phase (deg): 116.5 25.5 69.0 37.9 45.4 66.4 15.8 31.1 0.0 0.0 0.0
12 Instrument Data Meter Type: FG5 Meter S/N: 224 Factory Height: 116.40 cm Rubidium Frequency: 10000000.00000 Hz Laser: WEO (200) ID: 632.99117754 nm ( 0.36 V) IE: 632.99119473 nm ( -0.09 V) IF: 632.99121259 nm ( -0.46 V) IG: 632.99123023 nm ( -0.81 V) IH: 632.99136890 nm ( -1.45 V) II: 632.99139822 nm ( -1.20 V) IJ: 632.99142704 nm ( -0.90 V) Modulation Frequency: 8333.350 Hz Processing Results
Date: 12/14/05 Time: 14:55:16 DOY: 348 Year: 2005 Gravity: 978867771.63 uGal
Set Scatter: 16.84 uGal
Measurement Precision: 2.54 uGal Total Uncertainty: 37.94 uGal Number of Sets Collected: 44 Number of Sets Processed: 44 Set #s Processed:
1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44 Number of Sets NOT Processed: 0
Set #s NOT Processed: Number of Drops/Set: 200 Total Drops Accepted: 8673 Total Drops Rejected: 127 Total Fringes Acquired: 700 Fringe Start: 30
Processed Fringes: 650 GuideCard Multiplex: 4 GuideCard Scale Factor: 250 Gravity Corrections
Earth Tide (ETGTAB): 8.08 uGal Ocean Load: 0.29 uGal Polar Motion: 5.17 uGal Barometric Pressure: 3.22 uGal Datum Height: 115.05 uGal Reference Xo: 0.01 uGal Uncertainties
Earth Tide Factor: 0.500
Average Earth Tide Uncertainty: 4.04 uGal Ocean Load Factor: 0.20
Average Ocean Load Uncertainty: 0.06 uGal Barometric: 1.00 uGal
Polar Motion: 0.05 uGal Laser: 0.05 uGal Clock: 0.50 uGal System Type: 1.10 uGal Tidal Swell: 0.00 uGal Water Table: 0.00 uGal Unmodeled: 0.00 uGal System Setup: 1.00 uGal
Gradient: 0.89 uGal ( 0.03 uGal/cm) Comments
2-2 相對重力測量
相對重力測量是測定兩個位置的重力差值，也就是由一已知重力點位做為基 準，利用其與未知點位的重力較差求得該未知點位的重力值，此觀測原理與水準 測量的做法相似，相對重力測量的作業主要採往、返測，進行方式主要以階梯式 及星狀式為主，見圖 2-5。 圖 2-5：相對重力觀測進行方式(聯合勤務總司令部，1983；Torge, 1989) 本論文中所討論的相對重力測量作業，主要有兩個部分，一為範圍遍及全臺 灣的一等水準點上相對重力測量，觀測路線主要是沿各省、縣道公路施測，觀測 位置主要為內政部一等一級與二級水準點，觀測進行方式為階梯式(Ladder)，見 圖 2-5，以每一日可進行往返觀測的路程為一個觀測測段，多個觀測測段則組成 測線及環線，連結全臺灣所有觀測測線及環線形成全島的相對重力觀測網，見圖 2-6；另一為雲林地區相對重力網，觀測範圍僅於雲林縣境內，觀測行進主要是 以相對重力網中央位置所觀測的同安國小(TAES)絕對重力點為起迄位置，對週邊 各個選定的待測重力點做相對重力觀測，進行方式是混合以星狀式(Star)及階梯 式(Ladder)並連接引測網形內各重力觀測點，再對網形外圍各重力點做相對重力 觀測，見圖 5-4，如此可增加觀測數，並提升網形平差後的重力值精度。14 圖 2-6：全臺重力觀測網形圖(內政部，2003) 內政部於 2003 年 6 月委由國立交通大學完成全臺相對重力測量點位共計 1183 點，另又於 2007 年 6 月委由寰宇測量工程顧問有限公司完成全臺相對重力 測量點位共計 2198 點，雲林地層下陷監測部分是以小區域相對重力網做觀測， 來研究該區域的重力變化情形，本文將於第五章做討論。 相對重力儀使用的測量原理分為動力法（dynamic method）及靜力法（static method）兩種，現代的相對重力儀幾乎均採用靜力法原理。國內用來觀測相對重 力所使的儀器，主要有 Lacoste & Romberg G 型、Graviton-EG 型及 Scintrex CG-5 重力儀等三種。本文中所收集並計算的 2003 年全臺相對重力觀測資料，即是以 Lacoste & Romberg G 型重力儀做觀測，而 2007 年全臺相對重力觀測資料，是以
Scintrex CG-5 重力儀做觀測，而雲林地區相對重力觀測網，則是使用 Lacoste & Romberg Graviton-EG 型重力儀及 Scintrex CG-5 重力儀觀測收集數據。
2-2-1 Lacoste & Romberg G 型重力儀
Lacoste & Romberg G 型(LCR G)重力儀是近四十年來，世界使用最廣泛獲取
區域重力資料的標準儀器。重力儀外觀尺寸為 20
15
25 cm，重量約 4.2 kg， 另外包含外箱、電池、定平底盤和充電器，總重量約為 12.3 kg，見圖 2-7(LCR G Manual, 1997)。1
4
3
6
5
2
圖 2-7：Lacoste & Romberg G 型重力儀。1 為重力儀本體、2 為 12V 充電式電池、 3 為電池充電器、4 為儀器外箱、5 為定平底盤、6 為安培計 (聯合勤務總司令部， 1983) 該儀器在設計時，主要為了讓重力儀能適應於全球重力值觀測(約 977 gal～ 984 gal)，LCR G 型重力儀觀測相對重力的讀數範圍為 7000mgal，而讀數精度為 0.001 mgal，儀器的每月漂移量小於 0.5 mgal。為了減低大氣壓力改變對儀器造 成的影響，重力儀出廠時已完全密封，其內部並增加大氣補償以提高精度，其感 應器在消磁之後密閉於防磁盒中，見圖 2-8。
16
圖 2-8：Lacoste & Romberg G 型重力儀內部組件圖(聯合勤務總司令部，1983)
LCR G 型重力儀的操作原理，主要為儀器的重力感應系統及量測系統。重 力感應系統的部分，包含水平桿一端的秤錘，水平桿由一防震彈簧組成浮動式支 軸，用來消除觀測運動過程中的任何摩擦。由於重力感應系統完全懸掛彈簧上， 必須嚴密確保儀器外殼不受到劇烈撞擊及震動，以免導致系統損壞。量測系統部 分，是由槓桿系統與量測螺旋所組成，該系統須於已知的重力值之間精確校正， 校正參數(重力單位換算常數)則與量測螺旋及水平桿系統的品質相關，儀器的零 點漂移率為每月可小於 0.5 mgal。為了確保重力儀的保持最佳的工作狀態測量作 業前及作業期間，特別是每三週或長距離搬運後，必須對重力儀進行檢定及校 正。 重力儀使用時，須接上 12 伏特之電池，保持重力儀於攝氏 49 度至 51.5 度 間的最佳工作溫度狀態，電池充電器與電池連接充電時，須將電池移出儀器外箱， 以防電池過熱。未使用儀器時，則利用 110 伏特之交流電充電器保持儀器暖機， 以維持重力儀之工作溫度(魏祥鴻，2005)。由於此型重力儀是以傳統光學目鏡讀 取重力觀測數據，並記錄於觀測手簿中，另需詳細記錄觀測地點的經度、緯度、
高程及大氣壓力值，以便內業計算重力值之前，先行對重力觀測量做各式環境改 正，期能較理想的觀測數據，重力相關環境改正部分，留待第 2-3 節與第三章再 做討論。
2-2-2 Lacoste & Romberg Graviton-EG 型重力儀
Lacoste & Romberg Graviton-EG 型(EG)重力儀，見圖 2-9，為全功能、全自 動化且自動定平的相對重力儀。該儀器藉由內建伺服馬達驅動嵌腳即可達真正水 平，因此觀測的讀數已不必如其他電子重力儀需透過數值推算來補償傾斜，更不 必像 LCR G 型用旋轉刻度盤做觀測及校正，並免除了換算各觀測站重力值的繁 雜程序。此型儀器的優點計有：容易操作、自動定平、重量輕、全功能、高靈敏 度、整合資料讀取及可擴充性等。
圖 2-9：Lacoste & Romberg Graviton-EG 型重力儀(LCR Graviton-EG Manual, 2002)
EG 型重力儀的感應器是以零長彈簧為主，此感應器與 LCR G 型相同，見圖 2-10，觀測解析度出廠前可達到 0.001 mgal 的範圍，其重複性在可控制情況下可 達到 0.01 mgal，而在野外環境下則約 0.003 至 0.020 mgal，儀器漂移量每月可小
18
於 0.001 mgal，內建 32 MB 快閃記憶體，可以觀測記錄 100,000 站以上(LCR Graviton-EG Manual, 2002；魏祥鴻，2005)。
圖 2-10：Lacoste & Romberg Graviton-EG 型重力儀內部組件配置圖(LCR Graviton-EG Manual, 2002) EG 型重力儀於觀測時，可外接 GPS 裝置或直接輸入觀測位置的坐標經、緯 度及高程值，藉此儀器可計算出地球固體潮效應，直接對觀測重力值做改正，也 可設定觀測時間長度，在記錄觀測資料時，更可以濾波視窗(長度設定為 1、10、 30、60 或 120 秒)方式獲得最佳的觀測品質，並可顯示觀測量的標準偏差(LCR Graviton-EG Manual, 2002)。最後可將觀測數據下載至電腦，直接做為相關計算 參考，可省去紀錄手簿的繁雜並可避免人為錯誤。由於此型儀器可長時間連續觀 測，在每一筆設定觀測時間長度的數據中，可以檢視出該觀測時間長度所測得重 力平均值及標準偏差(RMSError)，並可於觀測進行時，直接檢視觀測時間長度內 各筆觀測重力的解析度(Error)。 根據實際觀測經驗， EG 型重力儀觀測所得的重力值標準偏差，可以明顯地 反應觀測環境的噪聲(Noise)，當觀測重力值的 RMSError 值越大時，表示環境的 震動干擾越明顯，其中不乏人為因素如交通或工業所產生的震動影響，若以自然
產生的地震來看，Yeh et al.(1982)研究中指出，當臺灣以外地區所產生地震後， 其震波傳至臺灣地區時，根據相對重力儀的觀測，可產生約 50~100 µgal 之重力 變化量。 圖 2-11：地震震波傳遞波形(Yeh et al., 1982) 本文利用EG型重力儀(S/N 1184)在彰化縣大城鄉、雲林縣麥寮鄉、崙背鄉與 臺西鄉等之沿海平原地區做觀測時，該地區為海拔約20至60公尺之沖積層地質帶 (經濟部中央地質調查所)，重力儀觀測量RMSError介於 0.434 ~ 2.686 mgal之間 不等，見圖2-12。 RMSError ( mgal ) 圖 2-12：本文臺灣西部沿海地區 Gravition-EG 型重力儀相對重力觀測量 RMSError 分布直方圖(2004/12)
20
另於嘉義縣梅山鄉與竹崎鄉進行觀測，而該地區位於中央山脈西麓，海拔約 在1100至1500公尺，地質為三峽群及其相當地層(經濟部中央地質調查所)，所得 到的觀測量RMSError介於 0.037 ~ 0.061 mgal之間，見表2-2。
表 2-2：本文臺灣西部中海拔山麓地區的相對重力觀測量 (單位：mgal) DATE TIME RAW GRAV TC GRAV RMSERR TIDE LATITUDE LONGITUDE ELEV
S330 1-Jul-05 05:53:05 2225.248 2225.228 0.037 0.0207 23.56218 120.7101 S330A3 S690 1-Jul-05 07:11:48 2134.537 2134.499 0.04 0.0378 23.43699 120.6579 S690A2 S690 1-Jul-05 07:17:31 2134.55 2134.512 0.061 0.0381 23.43699 120.6579 S690B1 S330 1-Jul-05 08:22:28 2225.362 2225.328 0.04 0.0339 23.56218 120.7101 S330B1 以上兩組觀測數據正好代表兩種不同地質條件之重力觀測量，震波發生在沿 海地區沖積層之軟質沉積地層時，會使得震波於行進時因為地質鬆軟而有共振現 象發生，進而使得震波有加大的趨勢(黃有志，2003)，沿海地區人口稠密農工商 業交通活動頻繁，因此易有各種微震發生，又由於靠近海岸邊，海浪衝擊岸線所 造成的震波也容易在此區域發生，綜合以上所產生的微震震波，乃是造成該區域 重力觀測量RMSError較大之主要因素。 相對於近中央山脈西麓地帶而言，地質為較堅硬之三峽群及其相當地層，震 波於此地質區域發生時，較易被地層所吸收，進而使得震波減小，又由於此觀測 區域離海岸較遠且人口較疏散，農工商業及交通活動不若西部沿海地區頻繁，因 此造成環境微震的因素較少，使得該區域重力觀測量RMSError遠小於西部沿海 地區。 另檢視此重力儀(S/N 1184)前往日本京都所做的觀測量表現情形，見圖 2-13， 更可以明顯的看出其 RMSError 遠小於在臺灣地區約十至數百倍，根據經濟部中 央地質調查所地質圖圖資顯示，臺灣西部沿海平原地質構造皆為沖積層，多為礫 石、砂與黏土所構成，屬於軟質地帶，又位於主要省道與線道公路旁，交通頻繁 並且常有載重車輛經過，因而使得觀測之環境噪聲過大而導致觀測量 RMSError 過高，導致 EG 型重力儀無法達到所預期觀測的 1 至 10 µgal 等級解析度。
RMSError ( mgal ) 圖 2-13：本文日本京都地區 Gravition-EG 型重力儀重力觀測量 RMSError 分布直 方圖(2004/11/23) 2-2-3 Scintrex CG-5 型重力儀 Scintrex CG-5 型(CG-5)重力儀相較於其他重力測量儀器，優點是更為簡易的 操作模式，且量測迅速、即時測量與統計分析等功能，可提高觀測精度與可靠度。 該儀器的讀數範圍為 8000 mgal，觀測解析度於出廠前檢定可達 0.001 mgal 的範 圍，其重複性在可控制情況下可達 0.005 mgal，而在野外環境下則可達約 0.005 mgal，儀器漂移量每日可小於 0.02 mgal (CG-5 Operation Manual, 2002)，見 圖 2-14。該儀器於測量前，可先輸入觀測位置坐標經緯度與高程，儀器即可於觀測 時自動改正由地球固體潮所造成的誤差，並即時將資料儲存於內建微電腦處理器， 為一自動化的重力測量儀器。CG-5 型重力儀精度高、量測速度快，更可自動記 錄讀數、自動補償傾斜、溫度、固體潮等改正，另外於觀測者在輸入測站之經緯 度、高程值後，更可進一步施行地形改正以供其他相關應用。本論文將 CG-5 型 重力儀相關作業模式與操作經驗整理後，收錄於附錄 B 做詳細說明。
22 圖 2-14：Scintrex CG-5 重力儀 表 2-3 為本節三型相對重力儀的各式諸元比較表，依作者使用 CG-5 型重力 儀的經驗來看，利用此型重力儀觀測所得到的重力值標準偏差，普遍較 EG 型重 力儀來得低，見圖 5-7。 表 2-3：本文各式相對重力儀的諸元比較表 儀 器 項 目
Lacoste & Romberg G
Lacoste & Romberg
Graviton-EG Scintrex CG-5 重複觀測精度 0.01~0.02 mgal 0.001~0.01 mgal 0.005 mgal 觀測讀數範圍 全球 7000 mgal 全球 7000 mgal 全球 8000 mgal
定平功能 手動 自動 手動 觀測值型態 乘以換算常數 後即為重力值 重力值 重力值 改正功能 無 傾斜改正 固體潮改正 傾斜改正 固體潮改正 漂移率 每月可小於 0.5 mgal 每月可小於 0.001 mgal 每日可小於 0.02 mgal
2-3 重力環境改正
由於在地球表面上任何一處的重力值會隨時間而改變，對於這些引起重力值 變化的環境因素大致可分為日、月引潮力所產生的地球固體潮、海潮變化造成的 海潮負載、氣壓改變、地球極移、地下水位及土壤濕度變化等。 本論文中的重力觀測量，其中就已包含了上述環境因素的重力影響量，而所 謂的點位重力值，是定義在無日、月引潮力、無海潮變化（即海水面為平均海水 面）、標準大氣壓、無地球極移、地下水位及土壤濕度為長時間平均值的“正常” 地球狀態下。上列任何一項若非屬正常狀態，則會產生重力的變化，因此必須由 重力觀測值中移除。地下水位及土壤濕度變化部分，由於觀測時間不長，各絕對 重力觀測時間僅一至三日不等，而各相對重力觀測時間更只有數分鐘至數小時不 等，在短時間內地下水位與土壤濕度均無明顯變化，對於各重力觀測量的改正而 言，此兩項均可忽略不計。本論文中所考慮之環境引起之重力變化僅為前五項， 其所使用的模式大部分可見 Melchoir（1983），Vanicek and Krakiwsky （1986）， Moritz and Mueller（1987），Torge（1989）等相關文獻。本文另僅就海潮負載效 應所產生的重力變化做介紹，留待第三章進行討論。24
2-4 重力觀測量平差計算
開始平差計算的首先步驟，是組成相對重力測量觀測方程式，這個方程式包 含環境改正後的重力值g、基準未知數N 與漂移參數0 d t( )。相對重力測量觀測 方程式如下： l v g N0d t( ) (2-3) 由相對重力的讀數，最後是以重力差值(la b, lb la)做為觀測量，重力讀 數差值所組成的相對重力觀測方程式如下： , , [ ( ) ( )] a b a b b a b a l v g g d t d t (2-4) 上式(2-4)中，va b, 為la b, 的殘差值，t 與a t 為觀測時間，原本(2-3)式中的基準未b 知數因為相減而相消除。因為觀測時間差很小，可用低階多項式來模擬漂移參數。 當有 n 個相對重力觀測量時，就可以計算出重力值及漂移參數，並且採用最小二 乘法來求解，並可由如下的矩陣來表式： 1 1 2 2 1 1 1 ( ) ( ) 1 1 ( ) ( ) 1 1 ( ) ( ) b n a b a b a n b a b a n n n b a b a n g g t t t t t t t t L V BX d t t t t d (2-5) 上式(2-5)中，L 矩陣共有 n 組觀測量，X 矩陣是待解的 u 個未知數，V 矩陣則是 觀測誤差。 觀測量的權為： 2 0 2 i i P (2-6) 2 0 單位權變方(P=1) 相對重力網的基本觀測量是相對重力值，即為兩個重力點的重力差值，並加上環境變化及系統誤差影響量，由於在網形中， 若是基準不足或圖形缺失 (configuration defect)造成秩虧現象(Datum defect)，會讓法方程式產生奇異而致無 法求解逆矩陣，解決的方法就是加入約制條件。 在網形中，若是任意給定一個已知點，平差後所得到的結果，即是相對於該 已知點的成果；相反的，若是沒有給定已知點進行平差，法方程式矩陣 N 會是 秩虧的奇異矩陣，將會造成無窮多組解。魯林成(1982)稱這類的平差問題是秩虧 自由網，或簡稱自由網。為了讓自由網的平差中能有唯一解，主要能滿足改正數 V 加權平方和為最小的條件(VTPV)，Koch(1987)與 Caspary(1988)曾提出以下兩種 特性： (一) Xˆ為未知數估計值，Xˆ的 NORM 最小，即為： minimum ) (X TX X 2 (2-7) (二) 點位精度之後驗協變方矩陣的軌跡(trace)最小，即為： X trace(Q ) minimum (2-8) 其中 Xˆ 2 0Q ˆ 未知數之後驗協變方矩陣 2 0 ˆ 後驗單位權變方 本論文中介紹以下兩種平差模式：自由基準平差與加權約制平差，其中引用 李莉華(2001)、Hwang et al.(2004)與魏祥鴻(2005)等文獻做說明。 2-4-1 自由基準平差（datum-free adjustment） 若組成的重力觀測方程式並未加上任何約制條件的話，且無圖形缺失的狀況， 此時的秩虧度（rank defect）則為 1。假設重力點數有 n 個，則重力儀參數會有 u-n 個，重力觀測方程式的係數矩陣 B 中的每一列前 n 個元素是由 0、-1、1 所組 成，以下表示：
26 By 0 (2-9) 1 1 1 0 0 0 n u n c T y (2-10) 上式(2-10)中，c 表任一不為零的常數，係數矩陣 B 的 Null Space 有一非零元素， 即表示係數矩陣 B 的秩虧度為 1，若不以加權約制平差而讓法方程式有唯一解， 則可加入下列基準條件（Koch, 1987）： ˆ T S X = 0 (2-11) 且 BS = 0 (2-12) 其中
1 1 1 0 0 0
ST (2-13) 為了滿足(2-7)式解為 NORM 最小解或自由基準解（datum-free）：
ˆ+ T T -1 T b T -1 X = B PB + SS B PL = N + SS U (2-14) 經由誤差傳播，可得Xˆ後驗協變方矩陣為：
0 2 2 0 ˆ ˆ ˆ T 1 T 1ˆ x Σ N SS N N SS N (2-15) 後驗單位權變方為： r u n T V PV 2 0 ˆ r = 1 (2-16) 求解得出Xˆ，滿足(2-11)，故： 0 ˆ 1 i k i g (2-17) 上式(2-17)中， ˆg 為第i i站重力估值（包含在向量Xˆ中），(2-16)式可解釋為所有 點的重力值之平均為零。因此，由自由基準解得到的 ˆg 值並非真正之重力值，而i 自由基準的改正數 V 是唯一解且有意義的（Koch, 1987）。假設： P X X X g ˆ ˆ ˆ g P B = B B (2-18) 其中 ˆX 是重力值， ˆg X 是儀器參數，則兩點位間改正後之相對重力p g為唯一解， 可表示為： ˆ ˆ g g P P Δg = B X = V + L - B X (2-19) 因此，可經由一個已知重力值來推求重力網中所有點之重力值，此即為自由 基準平差（datum-free adjustment）解法，相關重力成果詳見第 4-3-1 節。
2-4-2 加權約制平差（weighted constraint adjustment）
另一種排除重力網平差奇異的方法，就是加入未知數觀測方程式，或稱為加 權約制方程式，並以此來做為約制條件。利用已知重力點的重力值及其先驗標準 偏差來進行加權約制平差，其未知數觀測方程式如下： X X X V + L = B X ，L 權矩陣=X P X (2-20) 以 L 和L 組成之觀測方程式如下： X ˆ a X X X V B L V = = X - = BX - L V B L ，權矩陣= X P 0 P = 0 P (2-21) 令VTPV最小值，利用最小二乘法求解出未知數之估值 a Xˆ ： ˆa T T -1 T T X X X X X X X = (B PB + B P B ) (B PL + B P L ) (2-22) 觀測量改正數的加權平方和為： T T T X X X V PV = V PV + V P V (2-23) 後驗單位權變方為： 2 0 ˆ n u r T V PV (2-24)
28 其中 n 相對重力觀測之數目 u 未知數的個數 r 約制點之數目(未知數觀測方程式個數) 令 T X X X X B P B = P ，則(2-22)式可寫為： ˆa T -1 T T X X X X X = (B PB + P ) (B PL + B P L ) (2-25) 經誤差傳播，可得未知數 ˆa X 之變方-協變方矩陣如下： ˆa ˆ 2 T -1 0 X X Σ = σ (B PB + P ) (2-26) 在加權約制平差的過程中，約制用控制點設計矩陣BX為一特殊矩陣，其行 數（columns）和列數（rows）皆等於未知數個數 u，除了相對於約制用控制點未 知參數對角線元素值是+1 外，其餘元素都是零的 u u 矩陣。 PX是約制用已知點權值所組成的 u u 對角線矩陣，且PX與ATPA皆為 u u 矩陣。而PX是除了相對於約制用控制點未知參數的對角線上元素值為控制點權 值之外，其餘元素都為零的 u u 對角線矩陣。 假設只以第 2、3 點為約制點，且其先驗權分別為P 和2 P ，則： 3 2 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X P P P (2-27) 由上式可以得知，約制點的權可以直接加入由相對重力觀測所組成的法方程 矩陣 T B PB，其中為相對於約制用控制點未知參數的對角線元素。當約制點的權 為無限大時，平差後的值會與平差前所給的控制點值相同，相應殘差值為零，即
為固定該點，這就是加權約制平差（weighted constraint adjustment）解法。 另外，當約制點的權介於零與無限大之間時，先驗值變動的大小與權成反比。 若沒有圖形缺失的問題，此平差系統中僅約制一點即可。這種只約制一點的平差 做法，稱為最小約制平差（minimum constraint adjustment）解法。
2-4-3 顯著性測試 在同一個點位上，於兩個不同時間所測得兩筆該點位的重力值，首先得判定 這兩筆重力值是否相異，即檢視這兩筆的差異值是否為零。若這兩筆差異值不為 零，則表示此點位的重力值產生變化，或稱該變化為顯著。 由於這兩筆重力值皆為由觀測及平差計算所得，因此均內含誤差值，在此假 設g 、1 g 為同一點位的前、後次觀測重力值，其標準偏差分別為2 1、2，df 、1 2 df 則是在計算g 、1 g 時所對應的平差自由度，則利用以下的假說來測試重力值2 變化是否為顯著（Koch, 1987）： 0: 2 1 0 H g g 1: 2 1 0 H g g 測試因子如下： 2 1 2 2 1 2 c g g T T （2-28） 則H 為拒絕，即為變化顯著。 0 上式（2-28）中，Tc值是自由度為m(df1df2)及某一信心水平(1
)（如 95％）下，t 分布之(1
)臨界值，見圖 2-13，例如，當(1
)=95%，df1df210， 則TC 2.228（Mikhail, 1976）。 若值g 的標準偏差為未知，而所對應的自由度與觀測儀器皆與1 g 相同，則130 可假設1 2，df1df2 2df1，若 T 及 Tc滿足下式，則變化為顯著： 2 1 2 c g g T T (2-29) 假設ˆe， ˆe分別為平差所得的參數估值及標準偏差，以下假說可測試參數 估值是否為顯著： 0: 0 H e 1: 0 H e 若測試因子 T 滿足下式： ) ( ˆ ˆ m T e T c e
(2-30) 則稱參數估值為顯著，其中 ( )T m 的定義與(2-28)同，而(2-30)式中 m 為平差自由c 度。 圖 2-15： t 分布之機率密度函數及(1
)臨界值 ( )T m ，m 為自由度 c三、海潮負載改正
3-1 概述
重力測量觀測量之中，海潮負載(Ocean Tide Loading, OTL)效應是僅次於地 球固體潮(Solid Earth Tide)效應的週期性重力變化影響量。由於潮汐在近岸或海 島觀測位置的影響較顯著，一方面由於海洋在靠近岸邊的深度開始改變，使得近 岸的地殼受到海洋載重的不同而產生變位，進而影響重力的變化，另一方面則是 因為海洋的質量變化造成產生對觀測位置的不同引力。因此，海潮負載效應對於 重力的影響包括兩個部分：一為海潮質量吸引所產生引力影響的直接效應，稱為 牛頓引力效應 (Newtonian effect)；另一為地殼受到海潮負載導致變形的間接效應， 稱為彈性效應 (elastic effect) (黃鉅富，2012)。 海潮負載效應是利用平衡位概念，以格林函數形式推導出全球積分模式，來 計算地表各位置的海潮負載重力效應改正量。本章所介紹計算海潮負載效應的理 論公式，其中以 Farrell（1972）、Melchoir（1983）及 Agnew（1997）等著作論 述為主，並以 SGOTL 程式做為計算海潮負載效應理論的工具，對本論文中的絕 對與相對重力改正後結果做初步比較與分析，以下各公式推導與相關理論，摘錄 自黃鉅富(2012)論文做說明。
32
3-2 牛頓引力效應
圖 3-1：本文中海潮負載效應於觀測位置 p 與水體 q 之間關係(參照 Hwang and Huang, 2012 繪製) 由於牛頓引力效應是來自海水質量的引力，根據牛頓萬有引力定律，可將觀 測位置 p 以外的全球海水面上海潮質量進行引力積分，見圖 3-1，各水體 q 對觀 測位置 p 產生的直接引力位 W，公式如下：
2 2 2 2 , , 2 cos q w q e p p p q C C e p e p pq Gdm h R W r G d S R r R r
(3-1) 其中 pq
水體q
q, q
與觀測位置 p 間的弧距角度，
q, q
為 q 點所在緯度及經度 q h 水體q
於平均海水面上之潮高w 海水平均密度（1030 kgm-3）
G
牛頓萬有引力常數 e R 觀測位置 p 點之球面半徑C
以整個球體為範圍之球積分範圍 q d 積分單位，相當於dq cos qqd dq q，於球面上相當於 2 e q dS R d
p r 觀測位置 p 點至地心間的距離，近似於R
e+ H
p； H
p 為p
點正高 由於潮高是定義在平均海水面上，上式（3-1）採用球面法近似平均海水面 進行積分，其中觀測位置 p 點所在球面半徑 Re是參考高斯公式進行計算(Torge, 1989)，而水體 q 至地心距離近似於其球面半徑 Re，觀測位置 p 點至地心距離同 樣假設近似於其球面半徑及正高之和（Re + Hp）。 根據引力位與加速度的關係，對引力位 W 進行反向垂直梯度計算，牛頓引 力效應於觀測位置 p 點的重力影響量gNewtonian，可見下式表示：
Newtonian p, p, p w q Newtonian p, pq q p p C W W g H G h K H d r H
（3-2） 式(3-2)中的核函數KNewtonian
Hp,pq
是以高程相依並不考慮方向(黃鉅富， 2012)，整理後可得核函數：
Newtonian 3 2 2 , 1 2 p pq a b K H a ab (3-3) e p e R H b R 根據黃鉅富(2012)的分析結果，近岸的觀測位置重力影響量若要達到 1 μgal34 等級的精度，就應當進行此牛頓引力效應改正。
3-3 彈性效應
彈性效應（elastic effect）在重力的影響量上，乃是海水質量重新分布而改變 引力位（potential）並造成地殼微幅變形，進而影響垂直方向的重力加速度，可 以藉由高程相依的格林函數（Green’s function）與海水面上海潮質量重新分布來 進行積分計算，相關理論及計算公式，在 Farrell（1972）、Melchior（1983）及 Moritz and Muller（1987）等文獻中有介紹。本節根據 Moritz and Muller（1987）、Guo et al.（2004）及圖 2-1 呈現觀測位 置 p 與水體 q 的相關性，介紹並推導高程相依的地球彈性重力效應理論；格林 函數可由高程相依函數構成，以垂直位移V H
