以時間序列預測麥寮地區地層下陷之研究

國 立 交 通 大 學

土木工程學系

碩

士

論

文

以 時 間 序 列 預 測 麥 寮 地 區 地 層 下 陷 之 研 究

A Study of the forecasting Subsidence in

Mailiao Township Using Time Series Model

研 究 生：康家桂

指導教授：陳春盛

以時間序列預測麥寮地區地層下陷之研究

A Study of the forecasting Subsidence in Mailiao Township

Using Time Series Model

研 究 生：康家桂 Student：Chia-Kuei Kang

指導教授：陳春盛 Advisor：Chun-Sung Chen

國 立 交 通 大 學

土 木 工 程 學 系

碩 士 論 文

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 Master

in

Civil Engineering August 2005

Hsinchu, Taiwan, Republic of China

以時間序列預測麥寮地區地層下陷之研究

研 究 生：康家桂 指導教授：陳春盛

國立交通大學土木工程學系碩士班

摘 要

當地下水的抽取量大於補注量，地層因為水平衡壓力的關係，導致土壤的孔 隙中沒有地下水，產生地層往下沉降。台灣因超抽地下水，引起地層下陷和海水 倒灌，付出的社會成本很高，如本文研究之雲林麥寮地區，早期因養殖業，大量 抽取地下水導致地下水位急劇下降 20 至 25 公尺。近年來下陷雖已減緩但並未停 止，因此分析其近況將有助於進一步掌握下陷動態。 本研究透過地表及地下之監測方法掌握目前麥寮地區地層下陷之變化與機 制，並藉由時間序列分析，依照變數本身過去的資料所存在的變異型態來建立模 型之方法，以地表監測方法之日變化量、週變化量及月變化量對過去歷史資料進 行評估，就目前地層下陷之變化選定適合之預測方法建立預測模式。於資料量之 不足之情形可加入地下之監測方法以及雨量、地下水位資料，藉由複迴歸分析模 式提升模式之精度。日、週、月變化量之預測成果，以週變化量所得之 MAPE=18% 為最佳；但以月變化量加上季節性的影響，其 MAPE=12%，若是再加上地陷井 地層壓縮量、地下水位以及雨量資料進行複迴歸分析，其 MAPE=7%於 Lewis 所 劃分四等級中乃可歸為高精確度成果。

A Study of the forecasting Subsidence in Mailiao Township

Using Time Series Model

Student：Chia-Kuei Kang Advisor：Chun-Sung Chen

Department﹙Institute﹚of Civil Engineering

National Chiao Tung University

ABSTRACT

Land subsidences triggered by either natural or anthropogenic factors are considered one of the most serious environmental problems for Taiwan. The stratum starts to sink when the soil initially lied on top of the ground water has depressed after the ground water pressure failed to resort its balance.

Such disasters tremendously impact the environment of Taiwan. It is showed by this investigation that the Mailiao Township soil level has sharply decreased by 20-25 meters over the years, which can impose blame on the fishery industries around this area in the early days. [As a result, the groundwater had been overused, and caused the land above it to begin to sink into the ocean.] The land in south of Taiwan still keep descending at present.

To further understand the situation, this investigation analyzed the current status of subsidence in Mailiao by the application of the time series model, which involves the monitoring data collected on the surface and under the ground. The established time series model was used to evaluate the historical data by their daily, weekly and monthly variations. Furthermore, based on the current subsidence model, the subsidence trend can be forecast. The accuracy of this model has been improved by the reference of the precipitation and ground water level data. This research concludes that, the weekly data set yielded the best result among 2 other historical data sets where the MAPE=18%. Moreover, it is also the conclusion that the variation of the historical data of monthly had influenced by seasonal status would have the MAPE equal to12%.

誌謝

本論文承蒙恩師陳春盛教授在研究上的悉心指導，史天元教授、黃金維教授 給予我的鼓勵與幫助，不吝撥冗提供意見與指正，讓我能順利完成本論文並在學 業以及研究有所得，在此向三位教授致上最誠摯的謝意。同時也感謝口試委員楊 潔豪校長、李振燾主任與黃金維教授提供寶貴建議讓本論文更完善。 研究生兩年生活很快，一轉眼就步入社會為工作而忙碌，在交大必需感謝學 長偉嘉、小蜜、BOSS，在工研院提供我這麼多資源，以及Bernese超強的達哥和 大綱，沒想到曾經的這一步，成就了我的論文。俊寰學長的貼心小禮物，曾子榜 的愛的鼓勵，還有許多人在我對我的論文絕望失意時，給我的安慰，點點滴滴都 在心頭說不出。 最後，僅以本論文獻給我敬愛的家人，感謝他們也感謝大家。 藉由簡短的文字訴說我無盡的感激。

目錄

中文摘要……….i 英文摘要………ii 誌謝………...iii 目錄………...iV 圖目錄………...vi 表目錄……….viii 一、 緒論………..1 1-1 研究區域介紹與動機…………..………..1 1-2 研究目的…………..………..1 1-3 研究方法與架構…..………..2 二、 文獻探討………..4 2-1 地層下陷的成因及影響……….4 2-2 監測方法之分析……….7 2-3 預測模式之相關研究………...………..9 三、 資料分析………11 3-1 GPS 資料分析………11 3-1.1 GPS 固定站建置與運作……….………11 3-1.2 基線解算……….………17 3-1.3 座標轉換……….………24 3-2 監測井地層壓縮量之分析….………...………...26 3-3 其他相關資料……….………..………35 3-3.1 降雨量資料……….………..………..………...…35 3-3.2 地下水文特性分析.………..………..………...…37

四、 時間序列分析……….….….……41 4-1 時間序列組成成分………...……….………..……41 4-2 時間序列預測模式………...………...……44 4-2.1 AR(p)模式………..…45 4-2.2 MA(p)模式…..………...…46 4-2.3 ARIMA(p,d,q)模式………...……….………...47 4-2.4 ARIMA (p,d,q)×(P,D,Q)S 季節性自我迴歸整合移動平均模式...48 4-3 時間序列分析之建立…...……….………...50 4-4 預測評估分析…...……….………...53 五、 預估模式精度分析………55 5-1 時間序列成份………...…56 5-2 時間序列模式………...…58 5-2.1 GPS 高程差日變化量………..…...…58 5-2.2 GPS 高程差週變化量……….….…...…60 5-2.3 GPS 高程差月變化量……….….…...…61 5-3 擬合後之比較分析………...…67 六、 結論與建議………68 6-1 研究結論………68 6-2 研究限制與建議………69 參考文獻………..70 附錄一………..72 附錄二……….…….77

圖目錄

圖 1-1 研究流程圖……….…...2 圖 2-1 1995~2000 各年 GPS 連續觀測站基線分量座標高程向之平均精度...8 圖 3-1 GPS 固定站分佈圖……….…....14 圖 3-2 新興國小 GPS 固定站 MP1 訊號品質分析圖……….….15 圖 3-3 新興國小 GPS 固定站 MP2 訊號品質分析圖……….….15 圖 3-4 墾丁 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖………19 圖 3-5 陽明山 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖…………20 圖 3-6 北港 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖………21 圖 3-7 新興國小 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖……....22 圖 3-8 濁水溪沖積扇監測井之位置圖………..………28 圖 3-9 無線電波磁感應環分層式監測井構造示意圖………..29 圖 3-10 雲林新興國小監測井各感應環至地表之壓縮歷線……….31 圖 3-11 雲林新興國小地層柱狀圖及相對於井底之壓縮量……….32 圖 3-12 GPS 高程差和監測井地層壓縮量之相關性………...…33 圖 3-13 新興國小 GPS 固定站與地層下陷監測井下陷量比較表………34 圖 3-14 雨量和 GPS 高程變化之關係….………..36 圖 3-15 地層壓縮量與地下水位比較….……….…..40 圖 4-1 時間序列成分……….………..43 圖 4-2 預測流程圖（林靖等，2005）……….………..52 圖 5-1 不規則變動以及長期趨勢加週期變動時間序列成份分解…………...56 圖 5-2 季節性時間序列成份分解……….………..57 圖 5-3 GPS 高程差之日變化量原始時間序列圖（日）……….….58 圖 5-4 GPS 高程差自身相關圖（日）……….………..59 圖 5-5 GPS 高程差偏自身相關圖（日）……….………..59

圖 5-6 GPS 高程差日變化量一次差分時間序列圖（日）……….……59 圖 5-7 一次差分自身相關圖（日）……….………..……..60 圖 5-8 一次差分偏自身相關圖（日）……….………..…..60 圖 5-9 GPS 高程差之週變化量原始時間序列圖（週）………..….…..60 圖 5-10 GPS 高程差自身相關圖（週）……….………..60 圖 5-11 GPS 高程差偏自身相關圖（週）……….………..60 圖 5-12 GPS 高程差週變化量一次差分時間序列圖（週）………...61 圖 5-13 GPS 高程差自身相關圖（週）……….…………...…...61 圖 5-14 GPS 高程差偏自身相關圖（週）……….………...…...61 圖 5-15 GPS 高程差之月變化量原始時間序列圖………...62 圖 5-16 GPS 高程差自身相關圖……….………...…...62 圖 5-17 GPS 高程差偏自身相關圖……….………...…...62 圖 5-18 GPS 高程差月變化量一次差分時間序列圖………....…...62 圖 5-19 一次差分自身相關圖……….……….…...……..63 圖 5-20 一次差分偏自身相關圖……….………..63 圖 5-21 經一次差分與季節差分時間序列圖……….…………..63 圖 5-22 一次差分季節差分自身相關圖……….………...……...63 圖 5-23 一次差分季節差分偏自身相關圖……….……..………63 圖 5-24 GPS 日高程變化量預測分析圖……….……….………….64 圖 5-25 GPS 週高程變化量預測分析圖……….……….………….64 圖 5-26 GPS 月高程變化量預測分析圖……….……….………….65 圖 5-27 GPS 月高程變化量季節性影響預測分析圖……….…….…….65 圖 5-28 GPS 月高程變化量配合相關資料預測分析圖………….……….….66 圖 5-29 GPS 月高程變化量配合相關資料季節性影響預測分析圖……..….66 圖 1 Bernese 4.2 bpe 安裝概要流程圖……….…………..…73

表目錄

表 2-1 地層下陷之成因………..….………...…4 表 2-2 台灣地層下陷主要發生原因……….4 表 3-1 GPS 固定站平面與高程座標分佈統計表………...………22 表 3-2 地下監測井之特性比較……….………..………....26 表 3-3 地下水位站況……….………..………...39 表 4-1 七種 ARMA 模式的特徵……….…………..………...…51 表 4-2 MAPE 預測能力之等級……….………..…………....…54 表 5-1 日、週、月高程變化量精度評析……….…………...…67 表 1 BPE 作業程序表……….………...…...75

一、

緒論

1-1 研究區域介紹與動機 經濟部中央地質調查所之『濁水溪沖積扇之水文地質與地下水系統概念模 型』報告中指出，雲林縣地表下二百公尺深之範圍內，由上而下主要可劃分為富 水層及阻水層，而各富水層在上游段岩性均屬礫石，並且合而為一，缺乏顯著之 分隔，愈往下游分層愈明確，粒徑愈細，富水層厚度漸小，亦即各富水層末端應 是封閉的，並未與海連通，形同自然地下水庫，大部分地下水均經由此區域內之 河床和地表入滲而得，此地下水補注區之土地如不當利用，將降低入滲補注量或 增加地下水污染。而雲林縣地下水屬濁水溪沖積扇地下水區，由濁水溪、新虎尾 溪、北港溪等河川之沖積扇堆積而成，本區主要含水層豐厚，岩層孔隙率大，原 本為台灣地區的地下水資源區，惟由於本區的地質結構為尚未固結之沖積層，地 下水的補充量並不充沛，且養殖業者大量抽取地下水，造成地下水下降日益嚴 重。濱海地區雖然可獲得降雨與地表水之補注，但因為近年來地層下陷、海水入 侵與倒灌的緣故造成水質較差，為了獲得品質較穩定之地下水，當地之抽水深度 已超過地表下二百公尺，深層受壓含水層已受到相當程度的超抽。 早期雲林地區下陷中心皆集中在沿海地區，自民國85 年下陷中心逐漸往內陸 移動，民國88 年下陷中心集中在元長、土庫與褒忠，民國 90 年以後，最大年下 陷速率皆超過10cm/年，加上台灣高速鐵路目前正經過雲林地區最嚴重之下陷中 心，使得雲林地區地層下陷情況更應加以重視（水利署，2005）。地層下陷所在 區域易肇致國土流失，地利降低，生活環境品質下降，社會成本增加等之負面效 果，為防患沉陷問題的擴張與惡化，掌握地層下陷動態實為重要課題。 1-2 研究目的 抽水深度已有多深？地下水位的高低和雨量的多寡對地層下陷的影響量有 多大? 地層不斷下陷未來趨勢如何估計? 鑒於台灣西南地區地層下陷嚴重，目前 預測模式又多考慮土壤壓密及黏土之含量等問題，但這些因素卻是難以準確估 計，本文擬針對GPS固定站24小時長期接收資料之監測方法，以時間序列對其解 算求得之高程變化分別以日變化量、週變化量及月變化量數據加以探討，評估其 預估模式準確性，並加入監測井之地層壓縮量、地下水位、雨量資料三種變數， 藉由複迴歸模式配合時間序列分析期能提高預測精度。

1-3 研究方法與架構 時間序列分析和其他傳統分析不同的是，它不需藉助預測變數，僅依照變數 本身過去的資料所存在的變異型態來建立模型。欲以時間序列預測未來，須對以 往若干時日連續不斷所產生的時間序列，加以詳細的分析，以明瞭其變動的趨勢。 於本研究中主要探討方向為下列三項，研究流程圖參圖1-1： （1）欲藉由 GPS 的連續紀錄監測資料，藉由 Bernese 解算得到固定站之幾何高 （h），利用不同時間段的幾何高高差△h 與地層下陷監測井之相對於井底 300 公尺深之壓縮量資料比對，評估於地層於 300 公尺下是否有壓縮的現 象。 （2）欲藉由民國 90 年 7 月到 94 年 6 月底之雲林縣新興國小 GPS 固定站之觀測 資料之數據，選取適當之時間序列分析（Time Series Analysis）模式，評估 以日變化量、週變化量及月變化量三種型態預測之精度。 （3）月變化量將 4 年資料壓縮，數據減少，為提升預測模式準確性，加入監測 井之地層壓縮量、地下水位、雨量資料三種變數，使用時間序列分析（Time Series Analysis）配合複迴歸分析模式，由過去之資料為建立模型，預測雲 林地層下陷未來趨勢。 模式建立 資料蒐集與處理 確定研究範圍 選擇研究方法 文獻探討 實證分析 結論與建議 研究動機與目的 圖 1-1 研究流程圖

1-4 研究架構 本論文內容包括 5 個章節，茲分述如下： 第1章 緒論 本章主要針對本研究的區域、動機、目的、方法與架構做概略性介紹。 第2章 文獻探討 針對地層下陷的成因及影響、監測方法以及預測模式之相關文獻進行探討。 第3章 資料分析 採用工研院能資所提供GPS 固定站監測資料以及地層下陷監測井之監測數 據，以GPS 固定站之下陷量與站底的壓縮量作一相關性比對，並進一步分析中 央氣象局之雨量資料，水文水資源資料管理供應系統之地下水位資料。 第4章 時間序列分析 以自相關函數（ACF）及偏自相關函數（PACF）建構配適時間序列模式來 預測未來；本章主要介紹時間序列模式之數學公式，以及評估方法。 第5章 預估模式精度分析 將日變化量、週變化量及月變化量與加入複迴歸分析之預測結果，利用 MSE、MAPE 等方法評估模式優劣。 第6章 結論與建議。

二、 文獻探討

台灣地區由於經濟的發展，水資源的需求大幅增加，由於地下水取用方便， 平原區超抽地下水的情形非常普遍，以致西南沿海以及台北和宜蘭地區，相繼發 生嚴重的地層下陷問題。這些下陷地區所衍生的排洪問題與水土資源永久性損害 等問題，使社會付出巨大成本。 2-1 地層下陷的成因及影響 廣義之地層下陷泛指地層向下移動產生地面及地層沈陷之現象，其發生的原 因可概分為天然因素及人為因素兩種，說明如表2-1；而在台灣地區地層下陷區 域已近全島面積之十分之一，主要發生原因可分類如表2-2。 各項原因肇致地層下陷面積範圍各不相同，而造成之環境衝擊大致可歸類下 列數項： （1）地下水資源容易耗竭 地下水資源涵蓄能力降低地層發生壓密作用，其厚度減少、孔隙率減 低，所能涵養之地下水水量亦相應減少，造成下陷地區之地下水資源容易耗 竭。 （2）地勢低窪易生洪氾 地層下陷區因為地勢低漥，排水系統功能不彰，逢雨易積水不退。沿海 地區，因地勢低陷地面低於海平面，時有暴潮溢淹、海水倒灌之災損。 （3）建築物、工程結構物與維生系統設施損壞 地下水位下降導致地層有效應力增加，無論淺基礎、深基礎或橋墩與橋 台等，因沈陷率不均而損壞，自來水、瓦斯管陷與下水道等維生管線則因彎 曲損壞而洩漏 （4）海水入侵導致地下水鹹化 地層下陷區地下水水位面低於海水面，導致海水侵入地下含水層，造成 地下水鹹化，使得地下水可利用度降低，間接使得土壤鹽化，影響農作物種 植。

（5）濕地生態之轉變 河流水路受地層下陷區地勢低窪影響而改道，除易造成淹水外，部份土 地長期積水不退逐漸發展形成濕地型態。 （6）海水倒灌或海堤潰決 堤岸保護程度減小，維護費用增加，易引起海水倒灌或海堤潰決。 整體而言，地層下陷所在區域易肇致國土流失，地利降低，生活環境品質下 降，社會成本增加等之負面效果。 表 2-1 地層下陷之成因 天然因素 人為因素 1.地殼變動 1.抽取地下水或石油 2.地層自然壓密 2.深基礎開挖及採礦 3.地震 3.工程排水 4.火山爆發 4.結構物重力壓密 5.地盤本身之上升或下沉 5.淺覆蓋開挖 6.其他 6.其他 表 2-2 台灣地層下陷主要發生原因 原因 說明 地下水水位下降 地下水水位下降的原因主要為地下水天然補助量的減 少及地下水的超量使用，導致土體結構的有效應力增 加，發生壓密作用造成沉陷。 地面構造物興建 當地面承載負荷增加時，土體結構的有效應力也隨之增 加，進而發生壓密作用造成沉陷。 深基礎開挖或採礦 當開挖面支撐強度不足時，會引起一定面積之土體移動 產生地面陷落。 新生填土 包括新生填土自然壓密作用及其填築增加原地面荷重 發生之壓密作用。 其他原因 如地殼變動、車輛載重滾壓或地震等。

2-1.1 地層之抽水壓縮行為 地下水屬於天然循環之資源，當人為不當抽用，超出”安全出水量”(safety yield) ， 則 含 水 層 之 水 位 將 產 生 下 陷 的 現 象 ， 此 時 稱 為 地 下 水 超 抽 （Terzaghi,1948）。 抽水引起之地層下陷是複雜的地質災害，其形成之必要條件包括”超量抽用 地下水”及”地層中具有可壓縮之地質材料”，因此，地層下陷大部分發生在含水 豐富而地質軟弱之現代沖積地層，該類地層組織結構複雜，不易評估其壓縮性（柳 志錫，2004）。Terzaghi (1996)指出超抽地下水引起之地層下陷係來自地質材料間 之壓縮或壓密，並顯現於地表高程之降低，其原理為當地層中之地下水位降低 時，土體承受之有效壓力(effective stress)增大，並導致土壤產生壓縮現象。因此 對於地層下陷之偵測，除了調查地面之變動資訊外，探討地層之變形行為也很重 要，方可釐清確實的地層下陷機制，以輔助修正地層下陷之防治方向。地層下陷 監測井因為直接觀測地層的變形，尤其是分層式監測井之觀測分析資料，更可協 助區分不同地層間之變形。單信瑜（1998）指出西部沿海地區的地層下陷與含水 層砂質土壤的高壓縮性有密切的關係。 超抽地下水引起之地層下陷是水土資源開發失調之後遺症，如何正確評估地 層下陷與地下水開發間之互制行為及模式，是地下水資源開發過程中最重要之基 礎工作（柳志錫，2004）。而隨著全球定位系統 GPS 定位技術成熟及精度提升， 整合無線通訊設備元件之遠端監測系統變為可行(朱伯勛，2004)。

2-2 監測方法之分析 地層下陷監測主要可分為地表下陷監測及地下地層變形監測，其中地表下陷 監測主要包括水準測量及GPS 定位測量，其功能為實際偵測地表之陷落量，此 為瞭解地表沉陷最直接之方法；另外地下地層變形監測則是透過地層下陷監測井 之設置，實際監測地下土層之變形行為，所以透過地表及地下之監測方法將可確 實掌握目前地層下陷之變化與機制（水利署，2006）。 在過去進行地層下陷之監測工作時，大多採用直接水準測量之方法進行， 其過程費時費力。GPS定位系統經過二十多年的發展，其技術日臻成熟，並已成 為當前最具便利性、精密性與實用性之衛星定位測量技術。對於地層下陷的監測 工作也已廣泛採用（蔡明達，1998）。洪偉嘉(2002)指出GPS不僅具備測量快速 且利用GPS固定站監測地層下陷更可獲得主要下陷區的日變化量、週變化量及月 變化量，因此在台灣地區利用GPS監測地層下陷為確實可行的方法。 彭淼祥等（2001）施測範圍800平方公里，水準高程的測量精度達1 cm以內； 由同一時期兩個完全重複觀測的GPS網比較，內部精度高程達15 mm、平面精度 達5 mm、基線精度達1.6 mm。由兩個時期GPS網所得到之下陷量與兩個時期水 準所獲得之下陷量驗證精度，由8個驗證數據顯示，GPS與水準作業最大差值 11 mm，平均 4 mm。 郭隆晨（2001）的研究中，台灣 GPS 連續觀測站高程的精度依年份分析的 結果，顯示出精度有逐年提高的趨勢，參見圖2-4。這與 GPS 衛星分布、固定站 數增多及國際參考框架趨於完整與穩定有關。以在台灣基線長平均為200 公里估 算， 1995-2000 年的水平向平均精度約為± 2.5~5.7 ㎜，高程向平均精度約為 7 ~15 ㎜；而 2000 年的每週解精度水平向更可達約 ± ± 2 ㎜，高程向約 7~9 ㎜。 ±

圖2-1 1995~2000 各年 GPS 連續觀測站基線分量座標高程向之平均精度

下標1 表示每日解，下標 7 表示每週解；

2-3 預測模式之相關研究 預測技術在各領域皆扮演相當重要的角色，而精確的預測成果可以提供決策 者更多的資訊（王成財等，2002）。在高程方面，王德盛（1997）以卡門濾波法 （Kalman Filter）用於水準點位垂直變動情形之預測。 歷年來有不少學者專家就地層下陷總量與速率加以預估，多就地下水抽出後 孔隙水減壓，使有效應力增大導致地盤經力學上之彈性變形與壓密考量（單信 瑜，1997）。 柳志錫（2004）將超抽地下水引起之地層下陷的預測方法大致可分為經驗法 （empirical method）、半理論法（semi-theoretical method）以及理論法（theoretical method）等三類。 (1) 經驗法即是直接推求地層下陷與時間之迴歸關係，並利用外插方式預測未 來可能之下陷量（曹以松，1969）（Vega，1984）（鐘明松等，1997），一般 常用的迴歸關係有二次函數、指數函數及對數函數等。 (2) 半理論法利用地層下陷與相關現象之關係進行預測，例如 Kumai（1969） 與Yamamoto（1984）以抽水量多寡來預測地層下陷量；Lofgren（1969）、

Huang and Wu（1969）及廖日昇等（1990）則直接利用下陷量對地下水位 變化量之比值進行預測；Gabrysch（1969）認為下陷量對地下水位變化量之 比值與黏土之含量有關，因此近一步推導下陷量與黏土含量百分比之關係。 (3) 理論法則根據土壤的力學行為進行解析及預測，因超抽地下水產生之地層 下陷過程涵蓋地下水流動與土壤受壓變形兩步驟；耿慶志等(1994)曾對雲嘉 平原及離島工業區地下水與地層下陷，利用MODFLOW 地下水程式建立水 流模式，並利用Terzaghi 一維壓密理論解沉陷。杜富麗和歐陽湘(1995)利用 COMPAC 模式模擬預測雲林地區地層下陷，可能由於輸入參數部份再壓縮 係數有誤，導致部份模擬結果回彈量過高有違常理。洪明瑞等(1997)認為地 層中土壤分佈複雜，且黏土並非為主要壓縮層，遂以歷時為參數，用雙曲 線法預測下陷量。劉志純等(1997)以 MODFLOW 模擬水流，以 USGS 發展 出之INTERBED 模擬地層下陷量。劉等對砂土層下陷量採線性彈性計算， 認為此模式適合用在如雲林地區黏土層不連續之狀況。徐享崑等(1997)依據

改良之Terzaghi 單向度壓密公式，發展多層土壤壓密沈陷程式。其模擬預 測結果與觀測值比對沈陷量的急速增大與回彈反映模式可能有誤。賈儀平 和陳鴻泉(1997)利用 ABAQUS 程式分析台西地層下陷，土壤部份利用 Biot 三維壓密理論，利用試誤法調整模型致滿意後，進行預測。蘇苗彬等(1997) 利用Terzaghi 壓密理論，將土層以黏土層當量觀念分層分析，但用試誤法 求得參數。 台灣沿海地層下陷問題的嚴重後果正漸漸浮現。近年來許多研究者均致力於 了解地層下陷的現況、成因與補救方法，但往往無法確認引起地層下陷的因素， 而使的預測成效有限。就以往沿海地層下陷區之狀況判斷，在過去二十年中，地 下水位下降高達 10 ~ 20 m，對淺層土壤而言就可能已造成了約 2 m 以上大幅的 沉陷；然而這部份的下陷隨地下水位的停止下降而停止。目前沿海地區抽水深度 漸漸往下發展之中，就淺層土壤來說，其大規模的下陷應已停止。但目前沿海下 陷區地表仍持續以每年 15 cm 左右的速度逐漸下陷。受壓縮的土層應是深層黏 性土壤，即使抽水立即終止，其最終的下陷量目前仍無法準確地預估；亦即此時 抽水將造成的災害將在未來才會真正地降臨。然而，繼續大量抽水將引致大規模 的海水入侵，則除了地層解壓變形之外，黏土質地質材料因孔隙水鹽化的滲透壓 密也將使下陷更加嚴重。

三、

資料分析

3-1 GPS 資料分析

GPS 全球定位系統（Global Positioning System）在經歷過去二十多年的發展 過程，目前在各項理論及應用技術上已趨成熟，並已成為當前重要的空間大地測 量（Space Geodesy）定位技術。傳統的水準測量作業可提供高精度之高程成果， 但是高精度水準測量的作業過程費時費力，且不易自動化，相較於GPS 測量作 業的快速、方便，對於需高精度三維點位座標資料之科學應用而言，GPS 定位可 提供一個獲取、分析及建立此類大地觀測量之快速有效方法。 影響GPS 相對定位的誤差來源，主要有衛星軌道誤差、衛星及接收儀時錶 誤差、測站誤差、電離層延遲、對流層延遲、週波脫落、週波未定値（integer cycle ambiguity）、雜訊、多路徑效應及相位中心偏差等項，在處理資料、計算基線向 量時，須藉由差分方式或用有效模式改正，將上述誤差僅可能的消除。 本研究採用工業技術研究院及內政部於台灣地區所建置之GPS連續觀測站 的資料，以工業技術研究院量測中心（TNML）為座標參考站，同時解算內政部 七個衛星追蹤站陽明山站（YMSM）、北港站（PKNM）、墾丁站（KDNM）、太 麻里站（TMAM）、鳳林站（FLNM）、馬祖站（MZUM）、金門站（KMNM）、 以及工業技術研究院能資所設置之雲林縣新興國小固定站（YLSS），共9個站， 其地理分布位置如圖3-1。計算時間由民國90年7月至民國94年6月底止。 3-1.1 GPS 固定站建置與運作 GPS 固定站的設置與運作，主要可分為下列幾項：（1）站址的選定（2）固 定站之設計（3）訊號品質分析。 （1）站址的選定 地層下陷GPS 固定站站址的選定主要可分為主觀因素與客觀因素兩大項來 說明。客觀因素方面，固定站應具備有位置永久性、交通方便、附近有人看管。 主觀因素方面，因地層下陷GPS 固定站主要目的為監測地層下陷資訊，因此其 位址必須在地層下陷中心且附近最好應具備有地層下陷監測井、地下水位觀測井 及水準點。利用GPS 觀測搭配水準測量數據、地層下陷監測井與地下水位量測，

以綜合分析地層下陷的動態。根據上述條件，本計畫於雲林縣新興國小設置GPS 固定站。 （2）固定站之設計 關於GPS 固定站的設計主要包括 A、基樁設計 B、GPS 接收器與天線選擇 C、 環境設計（包含電力供應、不斷電系統、自動溫控系統、防突波裝置）D、遠端 傳輸及監控系統設計。（詳參洪偉嘉等，2000）。 （3）訊號品質分析 GPS 訊號會因電離層閃爍（Ionospheric Scintillation）、多路徑效應 （Multipath）、自然與人為雜訊對訊號的干擾、信號通過障礙物所產生的繞射現 象、天線增益與接收器本身之雜訊等，未模式化之誤差因素影響接收訊號的品質 （李振燾，2000）。因此對於 GPS 接收器所接收的訊號應採取品質分析與管制， 以確保GPS 後續解算成果的品質。 史天元（2001）GPS觀測所得之數據，有多種可能的原因致使觀測之數據品 質不符合需求。諸如，儀器故障、老化，測站之多路徑效應過高，或測站附近電 磁干擾嚴重，或有暫時性之遮蔽，致使訊號中斷等，均有可能影響數據之品質。 因此，在大量數據同步進行嚴密計算之前，檢查個別數據以確定其品質，均有實 質上對單站GPS觀測量進行品質檢核之需求。TEQC程式(Estey & Meertens, 1999)，該程式包含格式轉換(Translation)，編修RINEX檔案(Editing)，及品質控 制(Quality Control)，其為針對GPS單站之觀測量進行初步之檢核。由其mp1及mp2 均值，可以研判是否該測站有多路徑效應之問題。同時，由訊雜比亦可提供觀測 量品質之量度。 GPS 衛星定位測量之多路徑效應誤差，是指 GPS 接收儀接收之衛星訊號包 含直接路徑及非直接路徑兩部分。非直接路徑的訊號是指為衛星訊號傳輸過程 中，經對流層、電離層折射所造成的時間延遲效應（衛星多路徑效應）或接收天 線周圍環境經折射干擾造成測站多路徑效應，又可分成L1 載波的多路徑效應 （mp1）、L2 載波的多路徑效應（mp2）兩種。 不同路徑長，會使得天線受到不同路徑訊號之干擾，兩種信號的重疊將會引 起相位中心位置的變化，由於直接訊號與反射訊號存在著時間差與相位偏移量，

將會導致接收儀無法辨識真正的相位觀測量，嚴重時甚至訊號失鎖產生週波脫落 （黃德昌，2003）。 2 1 2 2 1 1 1 1 2 1 2 2 1 2 1 2 1 1 Φ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ₋ − + Φ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − − = Φ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + Φ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + − = α α α α P mp P mp （3-1） 雜訊量在此以公尺為單位 其中 P 為虛擬距離觀測量 Φ 為載波相位觀測量 2 2 1 _⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = f f α f1=L1之頻率=1575.42 MHz f2=L2之頻率=1227.60 MHz 本研究利用TEQC 軟體，對雲林縣新興國小 GPS 固定站測試訊號品質，測 試時間由民國90 年 7 月 1 日至民國 94 年 6 月 31 日止，測量的觀測數據如圖 3-2 及圖3-3。94 年度新興國小 GPS 固定站的 MP1、MP2 值約介於 0.6～0.8 之間， 除91 年中至 92 年初，因固定站附近興建体育館，影響固定站之透空度，造成 MP2 值異常，其它時段 GPS 訊號接收品質皆無太大變化。 有鑑於新興國小GPS 固定站該儀器年限已長達五年以上，加上該 GPS 監測 站位置靠海，長期受海風及季節氣候影響，導致儀器發生鏽蝕故障；該固定站於 93 年 1 月至 5 月間，因 GPS 接收儀器故障導致資料空缺，於同年 5 月進行儀器 維護；93 年 8 月至 94 年 4 月間故障頻率密集導致監測活動中斷無法順利取得監 測資料；建議更新儀器以避免地層下陷監測活動中斷以及提升監測資料的穩定 性。

0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 MP1 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 MP2 圖 3-2 新興國小 GPS 固定站 MP1 訊號品質分析圖 圖 3-3 新興國小 GPS 固定站 MP2 訊號品質分析圖

3-1.2 基線解算 GPS 定位測量主要分為兩個部分，一為外業資料的收集（硬體部分），二為 資料的處理（軟體部分）；硬體部分，多以外購專業的工具為主，而在軟體部分， 則可以分為自行開發及購買套裝軟體兩類（李振燾,2001）。目前 Bernese 軟體是 國際上學術界普遍使用的軟體之ㄧ，並附有原始程式，可依實際需求進行修改， 亦可加入新的模式。

Bernese 軟體是由瑞士伯恩大學天文研究所（Institute of Astronomy,

University of Berne）所研發的，以 Fortran 的語言攥寫而成，並附有原始程式，

可依實際需求進行修改，亦可加入新的模式。該軟體除了一般性GPS 資料處理

功能外，另有BPE（Bernese Processing Engine）自動化處理 GPS 資料、計算軌 道參數、極運動（Polar Motion）參數、地球轉動（Earth Rotation）參數、推求

區域性及全球性之電離層模式、對流層折射附加參數與模擬GPS 資料之功能。

本研究將固定工業技術研究院量測中心（TNML）為座標參考站，並蒐集內

政部七個衛星追蹤-站陽明山站（YMSM）、北港站（PKNM）、墾丁站（KDNM）、

太麻里站（TMAM）、鳳林站（FLNM）、馬祖站（MZUM）、金門站（KMNM），

以及工業技術研究院能資所自行設置之雲林新興國小固定站（YLSS），共9 個追

蹤站，採用學術軟體Bernese 4.2 版 bpe，配合 IGS 精密星曆進行精密解算以求得

GPS 測站三維座標，並且 畫 出 三 維 座 標 變 化 圖，確 認 GPS 追 縱 站 運 作 正 常 。基線解算軟體Bernese 4.2 版 bpe 必須建置在 Linux 系統上，其安裝流程與 執行步驟請參照附錄一。 圖3-4~圖 3-7 為使用 Bernese 4.2 bpe 求解 GPS 資料後，選定墾丁（KDNM）、 陽明山（YMSM）、北港（PKGM）三個內政部 GPS 衛星固定站於民國 90 年 7 月至民國94 年 6 月底每日之座標解，與新興國小（YLSS）解算成果做比較，其 中北港和新興國小同樣位於雲林下陷區域，而陽明山、墾丁分別位於台北及屏東 地區，由圖可知，四個GPS 固定站之水平位移除墾丁之 E 座標方向（E 座標最 大值與最小值差異約22 cm，往西偏移）並無明顯偏向，而新興國小與北港兩站 其垂直向高程卻有顯著變化（參照表3-1）。就高程部分，墾丁（KDNM）、陽明 山（YMSM）高程並無明顯的下陷趨勢；而北港（PKGM）與新興國小同樣位於 雲林縣，於圖3-6 中可查知北港（PKGM）站於每年 1~6 月有明顯的下陷趨勢， 17

7~12 月下陷趨勢減緩;而在本研究觀測時段民國 90 年 7 月至 94 年 6 月底期間內 北港累積下陷量約為10 cm，新興國小 GPS 固定站累積下陷量約 17 cm。新興國 小（YLSS）GPS 固定站於 93 年 1 月至 5 月間，因 GPS 接收儀器故障導致資料 空缺，同年5 月進行儀器維護；93 年 8 月至 94 年 4 月間故障頻率密集導致監測 活動中斷無法順利取得監測資料；目前該固定站運作恢復正常迄今且觀察其下陷 趨勢無異常現象。

墾丁(KDNM) N座標分布圖 2427973.44 2427973.45 2427973.46 2427973.47 2427973.48 2427973.49 2427973.50 2427973.51 2427973.52 2427973.53 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 N 座 標 ( 公 尺 ) 墾丁(KDNM) E座標分布圖 227485.20 227485.25 227485.30 227485.35 227485.40 227485.45 227485.50 227485.55 227485.60 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 E 座 標 ( 公 尺 ) 墾丁(KDNM) H座標分布圖 58.20 58.22 58.24 58.26 58.28 58.30 58.32 58.34 58.36 58.38 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 高 程 ( 公 尺 ) 圖 3-4 墾丁 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖 19

陽明山(YMSM) N座標分布圖 2784255.78 2784255.79 2784255.80 2784255.81 2784255.82 2784255.83 2784255.84 2784255.85 2784255.86 90年1月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 N 座 標 ( 公 尺 ) 陽明山(YMSM) E座標分布圖 307867.76 307867.78 307867.80 307867.82 307867.84 307867.86 307867.88 307867.90 90年1月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 E 座 標 ( 公 尺 ) 陽明山(YMSM) H座標分布圖 783.98 784.00 784.02 784.04 784.06 784.08 784.10 784.12 784.14 90年1月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 高 程 ( 公 尺 ) 圖 3-5 陽明山 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖

北港(PKGM) N座標分布圖 2608670.64 2608670.65 2608670.66 2608670.67 2608670.68 2608670.69 2608670.70 2608670.71 2608670.72 90年7 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 N 座 標 ( 公 尺 ) 北港(PKGM) E座標分布圖 179111.82 179111.83 179111.84 179111.85 179111.86 179111.87 179111.88 179111.89 90年7 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 E 座 標 ( 公 尺 ) 北港(PKGM) H座標分布圖 42.80 42.85 42.90 42.95 43.00 43.05 90年7 91年1 7月 92年1 7月 93年1 7月 94年1 6月 高 程 ( 公 尺 ) 圖 3-6 北港 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖 21

新興國小(YLSS) N座標分布圖 2626354.22 2626354.23 2626354.24 2626354.25 2626354.26 2626354.27 2626354.28 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 N 座 標 ( 公 尺 ) 新興國小(YLSS) E座標分布圖 170718.72 170718.73 170718.74 170718.75 170718.76 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 E 座 標 ( 公 尺 ) 新興國小(YLSS) H座標分布圖 25.55 25.60 25.65 25.70 25.75 25.80 25.85 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 高 程 ( 公 尺 ) 圖 3-7 新興國小 GPS 水平變位(N 及 E 方向)與高程變化(H)結果圖

表 3-1 GPS 固定站平面與高程座標分佈統計表 墾丁 KDNM N E H 標準差( cm) _{0.8 4.9 1.5} 平均值(公尺) _{2427973.50 227485.43 58.32} 最大值(公尺) _{2427973.52 227485.54 58.36} 最小值(公尺) _{2427973.47 227485.32 58.24} 陽明山 YMSM N E H 標準差( cm) _{0.8 0.7 1.0} 平均值(公尺) _{2784255.82 307867.84 784.08} 最大值(公尺) _{2784255.85 307867.89 784.11} 最小值(公尺) _{2784255.80 307867.81 784.00} 北港 PKGM N E H 標準差( cm) _{0.4 0.4 3.7} 平均值(公尺) _{2608670.69 179111.86 42.90} 最大值(公尺) _{2608670.71 179111.88 43.01} 最小值(公尺) _{2608670.67 179111.84 42.81} 新興 YLSS N E H 標準差( cm) _{0.9 0.4 5.5} 平均值(公尺) _{2626354.25 170718.74 25.71} 最大值(公尺) _{2626354.27 170718.75 25.83} 最小值(公尺) _{2626354.22 170718.73 25.59} 23

3-1.3 座標轉換 解算測站座標時，是採每日資料計算成一組座標，此座標是以地球幾何中心 為原點的三維空間卡式座標。為TWD97（X, Y, Z）須經程式轉換成 TWD97

(

φ, λ, h

)

；

( )

( )

( )

N P h a e P b e Z X Y − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ × × − × × + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − − φ θ θ φ λ cos cos sin tan tan 3 2 3 2 ' 1 1 （3-2） 其中 a 為參考橢球長半軸 ， b 為參考橢球短半軸 第一離心率 ₂ 2 2 2 a b a e = − ， 第二離心率

( )

₂ 2 2 1 ' e e e − = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ × × = + = − = − b P a Z Y X P e a N 2 2 1 2 2_sin tan 1 φ ， ， θ φ ＝緯度 λ＝經度 h ＝橢球面高程（Ellipsoid height） 而GPS 測量所測得之高程為一數學橢球面（如 WGS84 橢球體）為起算基準， 稱為幾何高（橢球高），幾何高是純粹的幾何量，代表從參考之橢球面沿法線量 至點位之高度距離量。在求得橢球面高程h 後，須轉換為以物理量之大地水準面 為起算基準之正高，一般通用以H 表示，代表從參考之大地水準面沿垂線量至 點位之高程。此大地水準面即定義了高程的參考面，通常乃以近似大地水準面之 平均海水面，作為高程參考基準台灣區即以基隆驗潮站之平均海水面為基準。 可用下式（3-3）求得大地水準面高程（Geoid height）： 大地水準面高程＝橢球面高程－大地水準面距橢球面高程 H ＝ h － N （3-3）

監測地層下陷的方法可依據不同時刻的高程值是否變動來做判斷，如（3-4） ∆d = H₁−H₂ （3-4） d ∆ 為下陷量值； H1為T1時刻高程值； H2為T2時刻高程值。 將（3-3）式代入（3-4）式

(

h H

) (

h H

)

h N H H d = − = − − − =∆ −∆ ∆ _{1 2 1 1 2 2} （3-5） 由（3-5）式可知正高和橢球高之間最重要的差異為大地起伏值的決定，然而大 地起伏理論上會隨時間而改變，但實際上改變量甚小，而小範圍內的測區中地形 變化不大，其相應之大地水準面的趨勢變化微小，因此可以假設 ，則橢 球高差△h與正高差△H應該相等。 0 = ∆N ∆d =∆H =∆h （3-6） 蔡明達（1998）在探討新測 GPS 高程與舊有水準高程之間進行基準化算所 遭遇之問題。由研究成果可知，大地起伏值之處理是影響其沉陷量分析成果之主 要因素之一，當分別採用「單點法」或「相對法」之沉陷處理模式時，「相對法」 可以減少大地起伏值的影響。在本研究中使用相對法進行後續應用。 25

3-2 監測井地層壓縮量之分析 濁水溪沖積扇自民國40、50 年代開始大規模抽用地下水作為灌溉用途，當 時地層已略為受壓變形，但緊接著在60 年代沿海地區亦密集抽用地下水作為養 殖用途，使地層下陷問題漸漸浮現，並在70 年代達到高峰期，80 年代以後因沿 海養殖型態改變，及海水養殖技術之推廣，地層下陷集中在仍大量使用淡水養蜆 之地區(成功大學水工所，2003)。因地層下陷嚴重，乃陸續開始設置地層下陷觀 測井。地層下陷監測井是探討地中土層壓縮機制的最有效方法，其原理為鑽鑿一 深度足以涵蓋土層受壓範圍之井體，並在井體內設置量測設施，對該井定期實施 觀測後即可瞭解井底至地表間土層之變形行為（柳志錫，2004）。 目前台灣地區主要下陷監測井分布範圍涵蓋宜蘭、彰化、雲林、嘉義、屏東 等地，監測井型式可分為磁感應環分層式、鋼索式、鋼管式及伸縮儀式（表3-1）； 上述這些不同型式之監測井之設置方法、年限及量測方式均不同，其中以磁感應 環分層式監測井之監測效果最佳，其定位精度可達1~5 mm，且可相當多層次及 彈性調整監測層次等，相當適用於複雜多變之沖積扇地質區之地層下陷監測（工 研院，2003）。 對於地層下陷監測井之設置深度，基本上應參考地質概況、地下水集中抽用 之含水層深度範圍以及監測井之基本功能限制等，以台灣西南部之濱海地區而 言，根據相關調查資料，濁水溪沖積層之沖積層厚度就可能超過1000 公尺(林朝 棨，民國52 年)，另外沿海地區之地下水取用深度大部約介於 200~300 公尺（陳 鴻泉，民國86 年），所以西南部沿海地區之地層下陷監測井設置深度一般建議 為300 公尺（柳志錫等，2000）。 圖3-4為濁水溪沖積扇之監測井分佈之位置圖。本研究所選定雲林地區新興 國小之監測井，其型式為無線電波磁感應分層式監測井（圖3-5），設置深度為 300公尺，其原理係利用無線電波感應的方式，將磁感應環預先利用鑽井方式錨 定於地下不同深度之地層，再利用無線電波偵測設備，以「月」為時間單位量測 每一感應環距地表之深度，進而分析各層次之壓縮量，以瞭解各層次之沈陷情 形。這種監測井之特點為可依據分佈選擇適當之監測深度並且可多層監測，其偵 測感應環深度之設備可讀至公釐（工研院能資所，2001）。

表 3-2 地下監測井之特性比較 監測井型式 優點 缺點 精度 無線電波磁 感應環分層式 (Multilayer magnetic Sensing subsidence monitoring wells) ¾ 測讀精度高 ¾ 耐久性佳 ¾ 適合大深度量測 ¾ 監測層次具彈性 ¾ 需以人工量測 ¾ 成本較高 1 mm 鋼索式 ¾ 簡單且容易設置 ¾ 成本低廉 ¾ 誤差較大，精度較低 ¾ 設備易受鹽化侵蝕 10 mm 伸縮儀 ¾ 可多層量測 ¾ 量測容易 ¾ 可自動化監測 ¾ 設置過程需非常謹慎，否則 極易影響監測之準確性 ¾ 設置成本較高 0.1 mm 鋼管式 ¾ 可自動觀測 ¾ 設置較容易 ¾ 易受地層摩擦力之影響 ¾ 僅單層監測 0.1 mm 27

資料來源:工研院能資所 圖 3-8 濁水溪沖積扇監測井之位置圖

資料來源:工研院能資所 圖 3-9 無線電波磁感應環分層式監測井構造示意圖

新興國小監測井各感應環至地表之壓縮歷線（圖3-10），自85/9至94/9之總 壓縮量(0~292 m)為34.1 cm，其平均壓縮速率約0.32 cm/月(約3.8 cm/年)；而94年 至9月為止壓縮1.2 cm，其壓縮速率約0.13 cm/月(約1.6 cm/年)。94年至9月之壓縮 速率甚低，某些地區甚至有回彈現象，顯示地層壓縮有趨緩之趨勢。圖3-11為雲 林新興國小地層柱狀圖及相對於井底之壓縮深度變化圖，圖上每一條曲線代表一 次觀測量，曲線斜率變化大即為壓縮較嚴重之處。新興國小監測井呈現比較平滑 的下陷曲線，像一條開口向上的拋物線，表示越接近井底處其地層壓縮率越大。 由(圖 3-11)雲林新興國小地層柱狀圖及相對於井底之壓縮量之資料，可知 越接近井底處其地層壓縮率越大，但因繪圖之參考點為深度277 m 的地層，故不 確定地層壓縮範圍是否超過300 m。而 GPS 的連續紀錄監測資料，藉由 Bernese 解算得到地殼至地表之高程值，若求出高程差，用以和監測井監測方法所得相對 於井底之壓縮量比對，藉以評估於地層於300 公尺下是否有壓縮的現象。 雲林新興國小GPS 固定站之下陷量與地層下陷監測井之地層壓縮量監測時 間是由民國90 年 7 月至民國 94 年 6 月底(圖 3-13)。地陷井之資料為以月為單 位，而GPS 固定站資料為每日解，但因每日之高程值變化大，較不容易看出其 趨勢性，故以一週為單位求其平均值，短期無資料部份藉由內插方式求取概略 值，長期無資料部份則予以空缺。該成果顯示新興GPS 固定站於本觀測時段內 累積下陷量約17.8 cm，監測井累積下陷量約 16.7 cm。取監測井每月觀測之地層 累積壓縮量值與GPS 對應之幾何高程差異量值進行迴歸分析，用以分析兩資料 之相關性（參圖 3-12），計算結果其相關係數達96.7%；地層下陷監測井的監測 數據乃是相對於井底之壓縮量，GPS 固定站則為地表至地殼的高程變化，由數據 顯示兩者的下陷趨勢與下陷累積量相當一致，故可推知絕大部份之下陷發生在地 表至深度300 公尺之間之地層，地層 300 公尺下並無明顯的地層壓縮現象，也印 證柳志錫（2004）研究地層總壓縮量與同時期之水準測量資料做比較，指出雲彰 地區絕大部份之地層下陷發生在地表至深度300 公尺間地層。

-40 -35 -30 -25 -20 -15 -10 -5 0 85/9/19 86/9/14 87/9/9 88/9/4 89/8/29 90/8/24 91/8/19 92/8/14 93/8/8 94/8/3 95/7/29 雲林縣新興國小監測井 各感應 環 至地表 之 壓縮量 （ 公分） 0~12m 0~22m 0~32m 0~42m 0~57m 0~72m 0~87m 0~97m 0~107m 0~117m 0~137m 0~147m 0~152m 0~162m 0~177m 0~187m 0~202m 0~212m 0~222m 0~237m 0~252m 0~277m 0~292m

雲林縣新興國小監測井

圖 3-10 雲林新興國小監測井各感應環至地表之壓縮歷線

31

0 50 100 150 200 250 300 0 5 10 15 20 25 30 35 40

相對於井底之壓縮量（公分）

深度 （ 公 尺 ） 85/9/19 85/12/3 86/12/11 87/12/2 88/12/8 89/12/6 90/12/12 91/12/23 92/12/12 93/12/30 94/9/29 0 50 100 150 200 250 300 0 5 10 15 20 25 30 35 40

相對於井底之壓縮量（公分）

深度 （ 公 尺 ） 85/9/19 85/12/3 86/12/11 87/12/2 88/12/8 89/12/6 90/12/12 91/12/23 92/12/12 93/12/30 94/9/29

圖

3-11 雲林新興國小地層柱狀圖及相對於井底之壓縮量

33

GPS和監測井地層變化相關圖

R2 = 0.9671 -20 -15 -10 -5 0 -20.00 -18.00 -16.00 -14.00 -12.00 -10.00 -8.00 -6.00 -4.00 -2.00 0.00 GPS高程差(公分) 監測井地層壓縮量( 公分)

圖 3-12 GPS 高程差和監測井地層壓縮量之相關性

GPS 固定站和監測井地層變化相關圖

3-3 其他相關資料 3-3.1 降雨量資料 降雨量資料來源為「中央氣象局」，為了配合地陷井資料分析以及後續應用， 選取的資料型態為統計每月之降雨量，以迴歸分析探討GPS 連續觀測站資料與 雨量參數間變化的關係。因為新興國小連續觀測站附近並無連續的雨量觀測紀 錄，所以在本研究中配合中央氣象局氣象觀測站的分布，以台南農改場雲林分場 逐日雨量資料，代表GPS 連續觀測站址周圍環境的變化。在進行迴歸分析時， 假設地層之下陷為持續下陷之線性關係，將GPS 連續觀測站各參數資料求其線 性之趨勢線，減去該參數之值，利用趨勢線和觀測值之間的差值，可以得到地層 於此一時間明顯之變化量，藉此探討雨量資料與GPS 資料之間的相關性。 圖3-14 可察知，降雨量於某一時期增加，其地層會有緩慢回升現象，雨量多寡 與地層下陷量的關係為落後指標，不能及時反應其變化，但仍具有影響力。 35

GPS高程變化之趨勢線和觀測值之間的差值 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 季 節 性 影 響 25.50 25.55 25.60 25.65 25.70 25.75 25.80 25.85 高 程_︵ 公 尺_︶ 季節影響 下陷量 線性 (下陷量) （cm） 0.00 100.00 200.00 300.00 400.00 500.00 600.00 700.00 800.00 900.00 1000.00 90年7月 91年1月 7月 92年1月 7月 93年1月 7月 94年1月 6月 季 節 性 影 響 -2.50 -2.00 -1.50 -1.00 -0.50 0.00 0.50 1.00 1.50 2.00 2.50 雨 量_︵ 毫 米_︶ 雨量(毫米) 季節性影響 （cm） 圖 3-14 雨量和 GPS 高程變化之關係

3-3.2 地下水文特性分析 地面降雨部分會經由岩層、土壤下滲，貯存於地下含水層內，形成地下水， 其蘊藏量與降雨多寡、土壤滲透性、地層狀況等有關，而地下水位在自然補注下 會維持一穩定的高度。一般而言，地下水位之變化常隨自然補注之時序呈季節性 之變化，就長期而論，並無明顯之升降趨勢。但若有人為持續超抽地下水，易使 其變化頻率增加，則地下水位將呈長期下降之趨勢，表示該地區之抽水量已大於 補注量。 台灣地區的雨量雖然非常豐沛，但由於降雨在時間及空間上的分布極不均 勻，可利用的地面水資源相當有限，因此除了大量興建水庫蓄水，以供枯水期間 之用水外，並於各地普遍鑿井抽取地下水，以補充地面水源之不足；台灣在平原、 台地等沖積層間有廣大的地下水區，蘊藏著相當豐富的地下水資源，但因長期大 量抽取地下水的結果，使得各地區的地下水位長期呈明顯的下降趨勢，而部分地 區的地下水位已降至海平面以下，尤其在沿海地區有逐漸往內陸延伸之趨勢，部 分地區的地下水位甚至已降至海平面 20 公尺以下。由於地下水位的下降，使得 地下水水壓減小，加上地面的荷重和土壤、岩層的壓密作用，便逐漸造成了地層 下陷。因此地下水位的長期下降，可以視為即將發生地層下陷的先期跡象，所以 瞭解地下水的升降對地層的影響是必要的。 水利署於雲林地區已陸續建立水位觀測井，進行地下水位觀測，由於自然的 地下水位係隨著豐水或枯水季節而變化，若有人為抽水現象則易使其變化頻率增 加。本研究地下水井水位變化的資料取自「水文水資源資料管理供應系統」，以 日為時間單位，記錄每日平均水位，最靠近新興國小之地下水位觀測站為「和 豐」，其有井深 126 公尺以及 226 公尺之資料（表 3-3），可配合地層下陷監測 井各地層之壓縮量，藉以分析不同地層受地下水位的變動影響之壓縮行為。 近幾年來的水利處水準點檢測資料與工研院分層觀測資料配合水文地質資 料顯示，台灣西部沿海的地層下陷主要壓縮土層並非黏土層，而是砂質與沈泥質 土層；尤以含水層的夾沈泥砂層為主要壓縮層(柳志錫和歐陽湘，1995)。彰化、

雲林沿海地區黏土層並不豐富，不僅在各深度層次都出現且厚度薄、連結性欠 缺；不構成引發大量壓密的條件。反倒是最大壓縮量都在抽水深度附近出現；從 早期較淺的50 m以上地層，到目前100 m附近及以下的地層。亦即，地層下陷的 主要機制是抽水深度附近的含水層材料之壓縮導致。 「水文水資源資料管理供應系統」提供深度126公尺及226公尺之地下水位 量，對照(圖3-11)可知此兩層皆屬於細砂土壤層，而含水層主要由砂性土壤組成， 因為排水性高，地下水下降所造成之壓縮屬於即時沉陷。在地層壓縮量與地下水 位變動的資料比對上(圖3-15)，地下水位深度126公尺與226公尺兩者間變化量相 當，而在地層壓縮量方面，深度222~237公尺處之變化較大。地下水位有較快速 的回升時，地層有些許回彈現象，但回彈量並不高（約在1 mm之間），但於民 國94年地下水位量有明顯的升高趨勢，但是對地層壓縮量並沒有回升現象。根據 這一個現象，可以做出兩點推論，第一是主要壓縮層為透水性高的砂土層，所以 回彈與地下水位回升的時間差極短；第二是，砂土層的行為非線性彈性，因回脹 量遠不及壓縮量大。

井名 和豐(1) 設置日期 19930312 距地面高(公尺) 1.4 井深(公尺) 125.8 水位高(公尺) 1.7 完鑿深度(公尺) 126.0 井名 和豐(1) 設置日期 19930330 距地面高(公尺) 1.4 井深(公尺) 226.1 水位高(公尺) 1.7 完鑿深度(公尺) 226.0 表 3-1 地下水位站況

四、 時間序列分析

時間序列是指一組隨著時間而記錄下來的觀測值，對時間序列之研究，稱為 時間序列分析。時間序列的資料往往不能以迴歸分析的方法來建立模型加以分 析，因為迴歸分析主要強調不同變數之間的關係模式，想要建立的是因果模型； 而時間序列中之各觀測值間通常都存在相關性，時間相隔越短之兩觀測值，其相 關性越大，時間序列並不滿足所謂『各觀測值為獨立』的必要假設。 GPS 高程差資料型態為隨時間連續觀察所累積之數值資料，欲經由過去時間 的數值變化趨勢來預測未來數值可能的變化，藉由時間序列單變量以及多變量分 析來進行預測模式。多變量分析是將地層下陷監測井所得之地層縮量、雨量資料 及地下水位量三種變數數列納入考量，使用 ME、MSE、MAE、MAPE、MPE 四項指標來評估模式。 4-1 時間序列組成成分 為分析時間序列資料的模式或趨勢，通常須先了解時間序列資料的組合成份 (component)。通常將時間序列寫成 ，其中 t=0，1，2…為下標代表時間；並 以時間為橫軸，將各時間點的觀測值描繪出，如此或可大略瞭解該變數隨著時間 而變動的趨勢。一般視時間序列由長期趨勢(long-term trend)、季節變動(seasonal variation)、週期變動(cyclical fluctuation)及不規則變動(irregular fluctuation)，四成 份所構成。 } {y_t （1）不規則變動（irregular fluctuation） 不規則變動（也稱 white noise）是在時間序列中將長期趨勢，季節變動以及 循環變動等成份隔離後，所剩下隨機狀況的部份。一般而言，長期趨勢，季節變 動以及循環變動皆受到規則性因素的影響，而只有不規則因素是屬於隨機性的， 其發生原因為﹕自然災害、人為的意外因素、天氣突然改變以及巨大變化等。

（2）季節性變動（seasonal variation）

季節變動是一種週期性變動，對 GPS 高程差而言週期為半年，每週期皆會 有循環性的變動產生，周而復始，呈現重複性之行為的序列，通常與自然的氣候 季節有關。

（3）長期趨勢（long-term trend）加週期變動（cyclical fluctuation）

時間序列依時間進行而逐漸增加或減少的長期變化之趨勢，在一較長的時間 內，呈現出遞減的趨向，相當穩定會隨時間呈現一個趨勢，在此趨勢為一線性 (linear)。 循環變動是一種圍繞趨勢線上下波動的情形如鐘擺般地循環變動，主要包含 四個階段:上升或擴張(expansion)、高峰(peak)、下降或衰退(recession)、谷底 (trough)。循環變動的週期大約二至十五年，其變動的原因甚多，而且週期的長 短與幅度亦不一致。 分析時間序列的初步工作，係將時間序列繪製歷史資料曲線圖（圖 4-1），根 據圖形觀察出單一時間序列的變動情況，亦可知多種時間序列變動的相互關係。 除了根據歷史資料圖可概略了解時間序列之變動情況外，尚可根據時間序列之四 種成份的不同結合方式，而提出了所謂的相加模型與相乘模型來進一步地分析時 間序列。茲將此二種模型簡略說明如下﹕ （1） 相加模型 假定時間序列係基於四種成份相加而成的。相加模型中，各成份彼此間互相 獨立，無交互影響﹔亦即長期趨勢並不影響季節變動。若以 Y 表示時間序列， 則其方程式為﹕ Y=T+S+C+R （4-1） T：長期趨勢(long-term trend)； S：季節變動(seasonal variation)； C：週期變動(cyclical fluctuation)； R：不規則變動(irregular fluctuation)。

（2） 相乘模型 假定時間序列係基於四種成份相乘之結果。相乘模型中，各成份之間明顯地 存在相互依賴的關係，即假定季節變動與週期變動為長期趨勢的函數。如以方程 式來表示此模型，即為﹕ R S C T Y = × × × （4-2）

平穩時間序列

未來值 未來值 未來值

線性趨勢之時間序列

線性趨勢與季節性之時間序列

時間序列數值

Time

圖 4-1 時間序列成分

4-2 時間序列預測模式

自我迴歸整合移動平均模式(AutoRegression Integrated Moving Average, ARIMA）乃為 Box&Jenkins（1976）所提出的，主要方法是對歷史資料分析，檢 視其自相關與偏自相關等特性，應用三階段模式建構過程，配適一個最佳的模 式，以進行資料分析與預測。自我迴歸整合移動平均模式（ARIMA）在各類型 的時間數列走勢當中，以穩定型的時間數列(stationary time series)較常見也較為

容易分析。一般對穩定型時間數列的定義為：若一時間數列 ，其 、 、 三者皆為不受時間 t 影響的常數，則稱 為一穩定型時 間數列。所謂 Box&Jenkins 預測方法，就是根據穩定型時間數列。 } {Xt E(Xt) ) (Xt Var Cov(Xt,Xt+s) {Xt} 時間數列的資料特性是『每個觀察值之間並非是互相獨立』，前後期的觀測 值Xt與Xt+1之間往往會有相關性存在。這與迴歸模式不同，在迴歸模式中，觀測 值之間是互相獨立的，所以不同觀測值Ｘt之間，假設無相關。但是在時間數列 分析中，觀測值間的自相關性質探討卻相當重要。因為前一期觀察值與這一期(或 下幾期)觀察值有某種程度的關係存在。若能計算出它們之間的相關程度， 則對 未來的模式建構與預測有很大幫助。 自相關函數(ACF)與偏自相關函數(PACF)的分析便是最常用的計算方法。在 對時間數列建構模式之前，需藉樣本來估計自相關函數與偏自相關函數之型式， 再由此資訊來建構配適時間數列的模式。 a.自相關函數(ACF) 給定一組時間數列 ，觀測值 與 之間的自相關係數(在時間上相差 k 期)為：

{ }

Xt Xt Xt+k

(

)

(

)

k t t X X k t t k t t X X Cov X X + + + = _{σ σ} ρ , , （4-11） 稱做時間數列

{ }

Xt 的自相關函數(ACF)，以符號ρk表示。對於穩定的時間數列而 言，式（4-11）時差 k 期的自相關係數可寫成

(

)

2 2 x k t t k X X E σ µ ρ = + − _（4-12）

因為Cov

(

Xt,Xt+k

)

=Cov

(

Xt+k,Xt

)

=Cov

(

Xt,Xt−k

)

，故對任意整數 k，ρk =ρ−k。

b.偏自相關函數(PACF) 給定一組時間數列

{

。對任意時差 k，當欲將 與 之間所有變數 的影響除去，再考慮 與 之間的自相關程度時，（4-13） 為條件相關係數

}

t X X_t X_t+k 1 2 1, + , , + − + t t k t X X X _L Xt Xt+k

(

)

(

_{( ) (}

₎

)

k t t k t t t k t t k t t t k t t X X X X X X X Cov X X X X X + − + + + + − + + + + = _{σ σ} ρ 1 2 1 1 2 1 , , , , ,... , , L （4-13） 稱做偏自相關函數(PACF)，可用符號ρkk來表示之，且稱ρkk為相差 k 期的偏自相 關函數。 ARIMA 模式可用 p、d、q 三個引子(argument)描述，p 為 AR 階數，d 為差 分階數，q 為 MA 階數。例如： ARIMA(1,0,0)=AR(1) ARIMA(0,0,1)=MA(1) ARIMA(1,1,1)=ARMA(1,1) ARIMA 模式經由 p、d、q 三個引子來描述，藉由 p、d、q 的變化，可產生 無數多的模式，每個 p、d、q 的組合就是一個模式。觀察一串時間序列之後，研 究者必須找出最佳模式來描述此一時間序列。 4.2.1 AR(p)模式 AR 模式假定外生變數僅受變數本身之影響，而不受其他變數的影響，故模 式由外生變數本身(Y)及誤差項(e)所構成。AR 之一般模式，稱為 p 階自我迴歸 模式，以 AR(p)表示。 AR(p)之中的 p 表示最大落差階數，故 AR(p)模式可描述當期數值受到前面 p 期影響之關係，統計模式如下： Yt = φ1Yt-1 + φ2Yt-2 + ...+ φPYt-p + et （4-14）

其中

Yt =第t期的觀察值

Yt−i=落後i期之數值，i=0,1,2,…,p

φi=i階自我迴歸係數，i=0,1,2,…,p

et=對應之誤差項，遵循NID(0,σ2)

B= backward shift operator

m t t m t t

Y

Y

B

Y

BY

− −

=

=

₁ 將此一模式以後退運算子表示，則可改寫如下： (1−φ1B −φ2 B2− ...−φPBp)Yt = et （4-15） 可再改寫如下： φ(B)Yt = et （4-16） 其中 φ(B)=1−φ1B−φ2B2−...−φPBp， （4-17） 即為 B 之多項式 4.2.2 MA(p)模式 時間序列之當期數值(Yt)，除了會受前數期(Yt−p)之影響外，也可能受當期誤 差項(et)或前數期誤差項(et−q)之影響。 移動平均模式之一般式，以 MA(q)表示之，q 為最大落差階數。MA(q)模式 在描述當期數值受當期及前 q 期誤差項影響之關係，統計模式如下： Yt = et−θ1et−1−θ2et−2−...−θqet−q （4-18） 其中 Yt=第t期之觀察值； -θj=j階移動平均權重係數，j=1,2,…,q； et−j=落後j期之誤差項，遵循NID(0,σ2)，j=1,2,…,q。

MA(q)亦可用後退運算子(B)表示之： Yt = (1−θ1B−θ2B2−...−θqBq)et （4-19） 若將(1−θ1B−θ2B2−...−θqBq)以隱函數θ(B)表示，（4-19）式可改寫為： Yt = θ(B)et （4-20） 4.2.3 ARIMA(p,d,q)模式 單純的自我迴歸(AR)模式及單純的移動平均(MA)模式，二者所描述的時間 序列各具特色，但是真實世界中，研究者面臨時間序列往往無法僅用單純的 AR 或單純的 MA 模式加以描述，必須混合二者模式，才有較高的配適度。此種混合 模式稱為自我迴歸混合移動平均模式(ARMA，Auto Regression mixed Moving Average)。如同 AR 與 MA 模式，ARMA 之一般模式，係以 ARMA(p,q)表示。理 論上，p、q 二值可設定為任意非負整數，但實務上以 0、1、2 最為常見。

ARMA(p,q)模式是由 AR(p)混合 MA(q)而成，統計模式如下：

Yt = φ1Yt−1+φ2Yt−2…+φpYt−p+et−θ1et−1−θ2et−2−…−θqet−q （4-21） 運用後退運算子(B)表示 ARMA(p,q)，則上式可改寫如下：

(1−φ1B−φ2B2−…−φpBp)Yt =(1−θ1B−θ2B2−…−θqBq)et （4-22）

上式中，左式為 AR 之 B 多項式，右式為 MA 之 B 多項式。這兩個多項式 除參數不一樣外(φ與θ)，其他型式完全一樣，參數之符號亦皆為負號。

ARIMA 模式與 ARMA 模式間最大不同，在於 ARIMA 將差分引入模式，以 處理非平穩序列。差分乃藉由差分運算子(∇, Differencing Operator)進行計算，原 理是將當期數值與前期數值相減。差分之階數(d)愈高，計算式將愈形複雜。 ∇d_{=(∇∇…∇)，即差分d階，d=1,2,…；} Yt−i=落差i期之數值。 差分運算對 ARIMA 模式非常重要，可以使非平穩序列轉為平穩序列，以適 用 ARMA 模式繼續分析。 因此，ARIMA 模式之建立，經常需要先將原始序列進行差分，使之成為平 穩序列。差分 d 階後再配上 ARMA(p,q)模式，即形成 ARIMA (p,d,q)模式。 綜上所述，ARIMA (p,d,q)模式可開列如下：

(

)

[

]

(

)

t q q t d p pB Y B B B e B B φ φ µ θ θ θ φ − − − ∇ − = − − − − − _L 2 _L 2 1 2 2 1 1 1 （4-23） 其中 Yt =第t期數值； µ =差分 d 階後之均值； p= AR 之階數； d=差分之階數； q= MA 之階數； φ1,φ2,...,φp=AR之參數； θ1,θ2,...,θq = MA之參數； et =誤差項，遵循NID(0,σ2)。 不同(p,d,q)之組合，可構成無限多種 ARIMA 模式。因此，研究者在分析一 個時間序列時，必須先要辨認此一序列適合何者 ARIMA 模式，方能進一步估計 參數及產生預測（林靖等，2005）。 4-2.4 ARIMA (p,d,q)×(P,D,Q)S 季節性自我迴歸整合移動平均模式 (Seasonal AutoRegression Integrated Moving Average, SARIMA)

季節性差分與連續性差分之目的均在消除趨勢。季節性差分乃對季節性跨距 s(Seasonal Span)，差分 D 階，常見之跨距有 4、7、12、24，分別適用之資料 為季資料、日資料、月資料、小時資料，遇季節性資料，可將季節性差分與連續 性差分共同運用。亦即： ∇d∆D_sY_t 其中 d=連續性差分階數； D=季節性差分階數； s=季節性跨距 。