Volume19, No3, January 2015, pp. 225-237
1國立交通大學土木工程學系 兼任助理教授 收到日期:民國 102 年 10 月 10 日
2國立交通大學土木工程學系 碩士 修改日期:民國 103 年 07 月 18 日
3國立交通大學土木工程學系 教授 接受日期:民國 103 年 12 月 24 日
4香港理工大學土地測量及地理資訊學系 助理教授
＊通訊作者, 電話: 03-5910279, E-mail: [email protected]
時域相關點雷達干涉技術應用於雲林嚴重 地層下陷區監測
洪偉嘉
1*黃大任
2黃金維
3張磊
4摘 要
雲林地區在地質條件上富含有易壓縮土壤，加上不當過量抽用地下水，導致該區域為臺灣目前最嚴 重的地層下陷區。傳統監測方法是設置多元化監測系統，該系統包括：全球定位系統（Global Positioning System, GPS）、水準測量、磁環分層式地層下陷監測井（以下簡稱地陷監測井）與地下水位井，分別從空 中、地面與地下不同面向監測下陷區域的變化。上述方法會因經費的限制影響監測點位的密度，進而影 響整體監測的精度。本文以時域相關點雷達干涉技術（Temporarily Coherent Point SAR Interferometry , TCPInSAR）對地表變形進行監測，該方法不僅可獲大範圍地表變形資訊，同時因為不需經由相位解纏
（Phase Unwrapping），避免可能存在相位模糊（Phase Ambiguity）的錯誤解算。
本研究使用從 2007 年 3 月至 2011 年 3 月共 15 幅的 ALOS 衛星影像，組成 29 對的影像對，獲得 TCP 密度約為 196 像素/平方公里，遠高於水準點密度 0.22 點/平方公里。透過 268 個水準點的交叉驗證，
TCPInSAR 與水準測量之垂直變形量的均方根誤差（RMSE）達到 0.6 公分/年，本研究驗證時域相關點雷 達干涉技術搭配 ALOS 影像能夠高精度與高解析度監測雲林嚴重地層下陷區的下陷。
關鍵詞：地層下陷、高速鐵路、雷達干涉技術、相位模糊
1. 前言
由於世界人口與經濟的持續成長，人類對於 水資源的需求亦不斷提昇，當地面水資源不足時將 轉而使用地下水資源以彌補不足，但不當超量抽用 地下水將會導致地層下陷，進而產生其他衍生性災 害，如地表變形、建物傾斜或龜裂、積水不退、海 水入侵與土壤鹽化等，此一問題已普遍存在於各已 開發或開發中國家（Gabrysch and Ronald, 2005;
Tosi et al., 2007; Abidin et al., 2008），而臺灣亦是其 中之一。臺灣西部平原地區由於經濟快速發展，水 資源需求大幅增加，加上地下水取用方便，使得平 原區超抽地下水的情形普遍，導致西南沿海地區嚴 重的地層下陷問題，其中又以雲林地區最為嚴重
（工業技術研究院，2011）。
就區域地質而言，雲林地區位於濁水溪沖積 扇的南側，地理位置北起濁水溪，南至北港溪，東 起斗六丘陵，西止於臺灣海峽，面積約 1291 平方 公里，是臺灣西部最重要的農業地區（圖 1）。本 區地層為全新世未固結的礫石、砂、粉土及黏土等 所組成的現代沖積層，厚度超過 200 公尺以上（陳 文福與江崇榮，1999），圖 2 為雲林地區地表至深 度 300 公尺間之水文地質概念分層圖，可視為由四 個阻水層及四個含水層所交互組成之水文地質架 構（賴典章等，2003）。扇頂的含水層主要成分以 粗顆粒土壤為主，例如：礫石及粗砂。扇尾的含水 層主要成分以細顆粒土壤為主，例如：細砂及黏土，
扇央的含水層成分則介於扇頂與扇尾之間，主要成
分以粗顆粒土壤與細顆粒土壤的土層交錯。由於扇 頂的粗顆粒土壤具有高透水性與高強度，因此地下 水位高度不易大幅下降，亦不容易產生顯著地層下 陷。扇央與扇尾地區因富含有易壓縮土壤，因此當 大量抽用地下水，容易導致孔隙水壓下降，進而引 發地層下陷，尤其在扇央地區，因為土壤組成為粗 砂、細砂及黏土交錯的薄層，當抽用地下水時，其 壓縮速率將會比扇尾地區快速（洪偉嘉，2009）。
根據水利署統計資料顯示，目前分布於扇央與 扇尾地區的地下抽水井口數超過 10 萬口以上（雲 林科技大學，2009），因此當大量抽用地下水時，
將導致地下水位下降，引發嚴重之地層下陷。為監 測雲林地區地表變形情況，經濟部水利署自 1992 年起於該區域設置多元化監測系統，此系統包括水 準測量、GPS 連續觀測站網與地陷監測井網。但因 傳統地層下陷監測方法皆為點狀測量，後續需利用 點狀測量獲得之下陷量，再應用克立金內插法內插 成下陷趨勢面，因此地層下陷整體的監測準確度與 佈設的點位密度有關。高密度監測點雖可獲得較高
量測精度，但相對投入的監測經費也較高，導致監 測點位密度受限於測量經費的規模，因此利用 DInSAR（Differential SAR Interferometry）與 PSI
（Persistent Scatterer Interferometry）技術，可以有 效率且經濟的降低監測成本，並獲得大範圍地表變 形資訊的監測方式，已成為一個重要的趨勢 （例 如 Galloway et al., 1998； Hoffmann, 2003）。但因 大氣效應造成之干涉誤差不易於干涉計算中消除，
造 成 應 用 DInSAR 技 術 監 測 地 表 變 形 之 限 制
（Hoffmann, 2003），以及 PSI 技術會受衛星影像數 量的限制，影響 PS 點的品質。因此本研究利用 TCPInSAR 技術，克服傳統 DInSAR 及衛星影像數 量 等 限 制 ， 並 以 新 的 方 法 選 擇 時 域 相 關 點
（Temporarily Coherent Point, TCP），以雲林嚴重地 層下陷區為監測區域，解算該區 2007 至 2011 年地 表高程變形速率，上述成果並與水準測量進行比對，
分析可能的誤差，並以 TCPInSAR 成果分析雲林 嚴重下陷區的下陷趨勢，該成果可以作為未來地層 下陷防治的重要參考。
圖 1 雲林地區位置圖
圖 2 雲林地區水文地質概念模型
2. 傳統監測方法成果
圖 3 為雲林地區的水準網、GPS 固定站網與地 陷監測井網的分佈圖。雲林地區水準測量範圍北起 濁水溪南岸，南至北港溪北岸，整個水準檢測網係 由 67 條主要水準測線連結組合成 13 個環線閉合網，
實際水準總里程數約為 450 公里，平均每 1.5 公里 設置一個水準點。
分析 2007-2011 年期間水準測量所獲得的年下 陷速率圖如圖 4，由圖中顯示，雲林地區主要下陷 區域分佈在中山高速公路以西，對照圖 1 顯示，主 要下陷發生在扇央與扇尾，該現象也與地質條件相 符，4 年最大累積下陷達到約 27 公分，主要下陷 位置集中在虎尾鎮、土庫鎮、元長鄉與褒忠鄉，而 高速鐵路正通過下陷中心，因此地層下陷也將會對
高速鐵路產生影響。
如同第一節所敘述，雲林地區的水文地質條件 複雜，各個含水層也具有不同特性，為能監測地層 下不同含水層的壓縮量，因此透過磁環分層式地層 下陷監測井（參考圖 5）可釐清地表下不同地層之 壓縮量。地陷監測井的原理，主要是依據現場鑽探 的資料與地球物理探勘的方法，在地層下不同深度 土層的分界位置設置磁環（每口地陷監測井設置約 20-26 個磁環），透過量測磁環位置的變化，可獲得 不同深度的土層壓縮變化量，該監測井的量測精度 經由量測不確定統計分析，可達到優於 1mm 的監 測精度（Hung et. al., 2012），目前在濁水溪沖積扇 則已設置 22 口地陷監測井（監測井分佈位置參考 圖 3）。
圖 3 雲林地區監測點位與水準網分布圖
圖 4 雲林地區 2007-2011 年平均下陷速率圖
由於地陷監測井目前量測的深度只達到 300 公尺，超過 300 公尺以下的土層便無法量測，因此 在雲林地區高鐵沿線共設置三個 24 小時連續觀測 之 GPS 固定站，其功能不僅可以彌補監測井在深 度上的不足，同時因為 GPS 可以獲得每日監測位 置的三維變化量，不僅增加整體的監測頻率，同時 可以作為水準測量與地陷監測井比對的基礎。本研 究利用水利署設置之 GPS 固定站，並與國內其他 相關單位共 20 個固定站，組成 GPS 聯測網，座標 參考站以工業技術研究院量測中心（TNML）為座 標參考站（全球 IGS 站），並採用 Bernese 軟體每 日進行解算，以分析雲林地區主要下陷中心的每日 變化。
圖 6 為高鐵沿線下陷最嚴重的土庫國中監測 站的水準測量、GPS 固定站與地陷監測井的監測成
果比較圖，由圖中顯示水準測量與 GPS 固定站之 成果相當一致，但其累積下陷量都明顯大於地陷監 測井，分析土庫國中之地陷監測井之數據如圖 7，
由圖中顯示，主要的壓縮深度發生在地下 200 公尺 以下，該深度之土層主要為砂層與泥層之薄層相互 交疊而成，因此只要抽水，因為厚度較薄，會有快 速壓縮的效果。另外，因為地陷監測井之深度為 300 公尺，使得若地層壓縮量發生在地底 300 公尺 以下則無法量測該壓縮量，但 GPS 固定站所量測 到的是地表到岩盤的壓縮總量，因此透過 GPS 與 地陷監測井量測成果的時間序列比對，可以發現兩 者隨時間增長而有明顯的差異量，該差異量即是 300 公尺以下土壤所發生的壓縮量，該資料亦說明 該地區在土壤的深層（超過地下 300 公尺）亦有抽 水行為。
圖 5 地陷監測井示意圖
圖 6 土庫國中之 GPS 固定站、地陷監測井與水準測量成果比較圖
圖 7 土庫國中地陷監測井年平均壓縮速率圖
3. 研究方法
雷達干涉量測技術（Interferometry Synthetic Aperture Radar， InSAR）用於偵測大面積的地層 下陷，具有大面積涵蓋而非點位性分布的優勢，其 優點包括：
(1) 主動性：測量時無需太陽光輔助，僅透過衛星 主動發射雷達波即可進行觀測，因此在夜間亦 能順利施測。
(2) 高穿透性：雷達波相對於可見光具有較優越之 穿透特性，因此較不受限於雲、霧及水氣之影 響，易於穿透大氣層。
(3) 高觀測精度：對於衛載 SAR 系統而言，系統 量測之地表變形相位變化比地表起伏相位變 化靈敏。因此，利用衛載雷達干涉技術於製作 數值高程模型時，精度為公尺級；而應用於地 表變形監測時，精度則可達到公分級。
(4) 高取樣密度：以取樣密度而言，InSAR 技術的 取樣密度遠高於水準測量，因此更能獲得較精 細之面狀地表資訊，同時透過時間序列的雷達 影像解析，使得 InSAR 技術能提供不同時間尺 度之地表變化情形。
過去前人的研究包括 Massonnet et al. （1997）
利用合成孔徑雷達差分干涉技術來偵測加州 East
Mesa 的地層下陷情形，驗證 InSAR 技術應用在地 層下陷監測的可行性，而後陸續有相關的研究應用 在地層下陷的監測，亦獲得良好成果（如 Galloway et al., 1998 ； Hoffmann, 2003 ； Galloway and Hoffmann, 2007； Hung et al., 2010； Hung et al., 2011；Zhang et al., 2011；Zhang et al., 2012）。
永久散射體（Persistent Scatterer, PS）的概念 與 處 理 方 法 最 先 由 Ferretti 等 人 提 出 ， 稱 為 Permanent Scatters TechniqueTM （ Ferretti et al., 2001）。其概念是在多張 SAR 影像中，以像素之強 度值為基準，尋找各影像中高相關性的固定亮點，
此亮點即為 PS 點，再將 PS 點的相位於一時間序 列影像中萃取出來進行相位解纏，以獲得長時間衛 星視角方向之地表位移量；而後多位學者陸續提出 其他改良之 PSI 技術（Colesanti et al., 2003; Hooper et al., 2004）應用於各方研究之中。由於此種有效 率且經濟的概念不但能降低監測成本，更可獲得較 為詳細的地表變形資訊，因此，雷達干涉測量技術 已成為目前地表變形與地上結構物監測的主流趨 勢。
然而上述方法需根據相位的穩定性來辨識 PS 相關點，以得到穩定的變形量，因此必須要有足夠 的衛星影像，一般而言必須要有 25 幅雷達影像
（Colesanti et al., 2003），Hooper et al. （2007）提
0
50
100
150
200
250
300
0 2 4 6 8 10
深 度( 公 尺)
年平均壓縮量 (公分/年)
2004 2005 2006 2007
2008 2009 2010 2011
次壓縮層
主壓縮層
出的 StaMPS/MPI 的解算方法，則至少需要 13 幅 衛星影像。另外 PSInSAR 在進行處理時，相位解 纏是一項關鍵的工作，常用的方法為對每一幅干涉 圖解纏。但因時域相關點無法被精確選定，在某些 相關點上可能含有較大的相位噪聲，使得這些點的 相位模糊（phase ambiguity）解算失敗，並導致整 體的解纏誤差。為了解決上述問題，Zhang et al.
（2011）提出了 TCPInSAR （Temporarily Coherent Point SAR Interferometry）技術，應用強反射體對 於罩窗尺寸的改變，並不會影響像素偏移量標準差 的特性，找出強反射體作為 TCP 點，並使用最小 二乘法進行參數解算，使用離群檢測方式估算殘差 大小來區分含有相位模糊與否的相關點，取代了容 易產生錯誤或失敗的相位解纏步驟，進一步估算出 地表變形速率。
3.1 ALOS 衛星雷達資料
本研究採用日本宇宙航空研究開發機構
（Japan Aerospace Exploration Agency, JAXA）的 ALOS 衛星拍攝的 L 波段 SAR 系統衛星影像為資 料來源，軌道編號為 track 447、frame 450，軌道涵 蓋範圍示意圖如圖 8 所示，影像格式為升軌、Fine 拍攝模式之多觀點 HH、HV 偏極複數雷達影像，
其影像幅寬約 73.9 公里、空間解析度約 30 公尺、
雷達平均入射角度（Incident angle）為 38.7 度。本 研究共使用 15 幅衛星影像，影像期間分佈在 2007 年 3 月至 2011 年 3 月之間，參數選用空間垂直基 線的上限為 1100 公尺，時間基線的上限為 900 天，
衛星影像資料如表 1，共組成 29 對影像對如圖 9。
表 1 衛星影像基本資料表 垂直基線長上限（m） 1100 時間基線門檻（day） 45/900
影像時間 2007/03/04-2011/03/15 影像數（Images） 15 幅
影像對數
（Image pairs） 29 對
圖 8 台灣地區 ALOS 衛星軌道涵蓋範圍示意圖
圖 9 雷達干涉影像對空間與時間基線分佈圖
3.2 時域相關點雷達干涉技術
本 研 究 依 Zhang et al. （ 2011 ） 提 出 之 TCPInSAR 演算方式進行解算，其處理過程依序 為：
A. 挑選 TCP 候選點
TCPInSAR 利用 DInSAR 所解算的強度影像，
以移動式罩窗在衛星行進方向與垂直方向對各組 主副影像對的每個像素，利用公式（1）進行影像 偏移量標準差計算（Bamler, 2000）。初始使用 256
×256 罩窗，計算每個像素偏移量 O（i,j），O（i,j）
為影像中任一像素之偏移 量，其中 i=1,2,…l；
j=1,2,…m，並組成此影像對的偏移量矩陣 O（l×m）
後，利用直方圖統計可得到影像的偏移量峰值 Oc，
利用公式（2），設定 A 為初始門檻，刪除偏移量 過大的像素，後續改變罩窗尺寸進行迭代計算（由 4×4 至 64×64），則每個候選點可得到一組偏移量序
列，再計算此序列之標準差，當標準差小於套合精 度（如 0.1pixel）時，保留其為 TCP 候選點。
𝜎𝑥,𝑖 = √2𝑁3 ∙√1−𝛾𝜋𝛾2χ32 (1) 𝜎𝑥,𝑖表示第 x 組影像對中第 i 個像素的偏移量 標準差；N 為罩窗尺寸；χ為過取樣參數，於計算 中可視為常數；γ為影像對於此罩窗的相關係數。
𝑂𝑐 = 𝑝𝑒𝑎𝑘(ℎ𝑖𝑠𝑡2𝐷(𝑂𝑙×𝑚)) (2)
|𝑂𝑖,𝑗− 𝑂𝑐| < 𝐴 (3) ℎ𝑖𝑠𝑡2𝐷(𝑂𝑙×𝑚)表示沿衛星行進方向與垂直方向，利 用像素偏移量各取一直方圖統計。
B. TCP 點
使用一多項式擬合步驟 A 所獲得 TCP 候選點 之偏移量序列，若某一像素的偏移量偏離該多項式 時，則剔除該像素，計算得到所有副影像的 TCP 候選點後，取其交集，於每幅副影像中皆出現的點 才選定為最終的 TCP 點。
C. 最小二乘法求解
（1）觀測方程式
假設 N 幅影像生成了 M 幅滿足基線要求的差 分干涉圖，第 i 幅干涉中相關點 p 上的差分相位可 表示為公式（4）：
𝜙𝑝𝑖 = 𝑊{𝜙𝑡𝑜𝑝𝑜,𝑝𝑖 + 𝜙𝑑𝑒𝑓,𝑝𝑖 + 𝜙𝑎𝑡𝑚,𝑝𝑖 + 𝜙𝑜𝑟𝑏,𝑝𝑖 + 𝜙𝑁,𝑝𝑖 } (4) 𝑊* +為代表 InSAR 所測得的相位；𝜙𝑡 𝑝 ,𝑝𝑖 為地形誤 差造成的相位變化；𝜙 𝑒 ,𝑝𝑖 為目標物變形產生的相 位變化；𝜙𝑎𝑡 ,𝑝𝑖 為大氣影響導致的相位變化；𝜙 ,𝑝𝑖 為軌道誤差產生的相位變化；𝜙 ,𝑝𝑖 為相位噪聲。
依照公式（4），相關點 p 與相關點 q 的相位差可 寫為公式（5）：
𝜙𝑝, 𝑖 = 𝑊
{
𝜙𝑡 𝑝 ,𝑝, 𝑖 + 𝜙 𝑒 ,𝑝, 𝑖 + 𝜙𝑎𝑡 ,𝑝, 𝑖 +𝜙 ,𝑝, 𝑖 + 𝜙 ,𝑝, 𝑖
}
(5) 則在第 i 個干涉圖中的觀測可以表示為公式（6）：𝑖 =[ 𝜙𝑖 𝜙𝑖 𝜙𝑖] , 𝑖 = (6)
進一步地，所有干涉圖中的觀測為公式（7）：
（ × ）× =, - (7)
（2）軌道誤差
由於干涉圖中的軌道誤差可以用低階多項式 來擬合，考慮到干涉圖是雷達影像的線性組合，因 此干涉圖中的軌道誤差多項式亦可以看作兩個多 項式的線性組合。因此我們可以用一系列與影像相 關的多項式來組合出對應影像參與的干涉圖的軌 道誤差多項式。進一步地，考慮到多項式組合的相 對性，我們假設一幅影像的多項式各個係數為零，
即沒有軌道誤差貢獻，則其他影像的多項式可以表 示為公式（8）：
,𝑠 ,𝑝 = 𝑎 + + , = − (8) 在第 i 個干涉圖中，軌道多項式對相關點對上的相 位差的貢獻為公式（9）：
𝑠
×
= 𝑎 + + , = − (9)
其中 , 及 分別為相干點對的像素座標差。
以矩陣的形式可表示為公式（10）：
𝑠
×
= × 𝑠 , ×
, = − ( 0)
其中
× =
,
-
, 𝑠 ,× =
[
𝑎]
。 存在矩陣 A 能夠指示干涉圖 M 與雷達影像 N 的組 合關係，刪除參考影像所在的列的矩陣具有如公式（11）形式：
= [
− 0 0 0 − 0
0 0 −
]
×（𝑁−1）
(11)
則反映軌道誤差多項式與軌道誤差相位關係的設 計矩陣可表示為公式（12）：
( × )×((𝑁−1)×3)𝑜𝑟𝑏 =
×(𝑁−1) ×3 ( )
其中 為 Kronecker 張量積。與軌道誤差有關的相 位貢獻可表示為公式（13）：
( × )× = (( − )× )×
(13)
式中 包含了所有相對於參考影像的軌道誤差 多項式係數。
（3）DEM 誤差及變形速率
干涉相位與 DEM 誤差及變形速率參數之間的 關係在大多數 MT-InSAR 文獻中已被充分討論
（Zebker et al., 1997; Bombrun et al., 2009）。這裡 直接給出 M 個干涉圖中，在任意相關點 p 上相位 表達式為公式（14）：
𝑡 𝑝 + 𝑒 ,𝑝=
× [ ℎ𝑣𝑝
𝑝] (14) 式中 × 為設計矩陣反映了地形誤差與變形速率對
相位的貢獻。假設相關點對 G 與 P 個相關點直接 的關係由矩陣
×（ −1）表示，則在第 i 個干涉圖中，
由地形誤差及變形引起的相位差可以表示為公式
（15）：
𝑡𝑜𝑝𝑜 𝑑𝑒𝑓𝑜𝑖 = (
×（ −1）
⨂
𝑖1×2) 𝑝𝑎𝑟
(2×( −1)×1) , , 進一步地，所有干涉圖中所有點對的相位差可以表 示為公式（16）：
𝑡 𝑝 + 𝑒 ( × )×
= 𝑝𝑎
( × )×( ×( − ))( ×( − 𝑝𝑎 )× )
其中 𝑝𝑎 =, - 。 (4）觀測模型及參數解算
地形誤差及變形速率直接的關係可以表示為 公式（17）
( × )× = +
（ × ）× (17) 其中 = [ 𝑜𝑟𝑏 𝑝𝑎𝑟]， = [ 𝑜𝑟𝑏 𝑝𝑎𝑟]；
（ × ）×
包含了所有未參數化的相位分量，例如部份大氣延 遲及失相關噪聲。由於模型為大型稀疏線性系統，
Golub 及 Kahan（1965）提出的 LSQR（least squares residuals）方法將被用於參數解算。值得注意的是，
該模型中的觀測為未纏繞的相位資料（wrapped phase），參數解算過程需要對相位模糊進行探測 和剔除，後續透過相位差的計算，則可以推算所有 TCP 的變形量。
4. 研究 TCPInSAR 成果
TCPInSAR 所獲得為沿 LOS 方向的平均變化 率，為後續與水準測量進行比對，假設本區地表無 顯著的相對水平向變化，將 LOS 方向之變形速率 轉化為年平均垂直變形速率，並利用克立金內插法 將 TCP 點所獲得的年平均垂直變形速率，內插成 網格點，並使用 GIS 方法製作成年平均垂直變形 速率圖（圖 10），與水準測量成果（圖 4）進行比 較。
比 較 兩 種 測 量 方 式 獲 得 之 點 位 密 度 ， 由 TCPInSAR 所獲得的 TCP 密度達 196 像素/平方公 里，而水準測量樁之點位密度僅 0.2 點/平方公里，
因此圖形上也顯示 TCPInSAR 可提供較高解析度 之下陷變形量。為了有效驗證 TCPInSAR 的監測 精度，將研究區域內所有水準點資料（共 268 個）
與 TCPInSAR 成果進行比較，經計算兩者的均方 根誤差（root mean square error, RMSE）為 0.6 公分 /年，其中共有 200 個水準點（佔整體水準點數量 的 75%）的高程差異量在±1 公分/年以內，最大誤 差為 1.9 公分/年，統計 TCPInSAR 與水準測量差 異量的直方圖如圖 11，整體大略呈常態分佈。圖 12 為上述差異量的分佈圖，由圖中顯示，位於西 南角的鄉鎮（包括北港鎮、水林鄉與口湖鄉等地區）， TCPInSAR 與水準測量差異量較大。
(15)
(16)
圖 10 InSAR 成果換算為年平均垂直變形速率圖
圖 11 TCPInSAR 與水準測量差異量統計值方圖
圖 12 TCPInSAR 與水準測量差值大於±1 公分/年之水準點位置分布圖
分析可能造成 TCPInSAR 與水準測量的誤差 原因有下列三點：
（一） 衛星波長
本研究使用 ALOS 影像為 L 波段，波長為 23.6 公分，在 InSAR 衛星中屬於長波長衛 星，其優點在於高植被區域亦獲得點位，
但對於較小變形量的區域相對比較不敏 感。
（二） 施測時間
由於監測的時間有四年的時間，但本研究 僅收集到 15 幅的 ALOS 衛星影像，因此 在影像數目上，仍顯不足，同時因為測試 的影像數目不多，因此在資料的比對上會 有時間的差異，例如水準測量的期距為 2007 年 8 月~2011 年 5 月，而 InSAR 的期 距為 2007 年 3 月~民國 2011 年 3 月，由於 水準測量與 InSAR 兩者所施測的時間並 不完全一致，亦會導致兩者所獲得的下陷 量有部分差異。
（三） 水平向位移
由於 ALOS 衛星的雷達平均入射角度（約 38.7 度）比 C 波段衛星（例如 ERS 及 Envisat 衛星，約 23 度）大，因此當研究 區域內有水平向位移時，要將 LOS（Line of Sight）方向投影至垂直方向，必須考慮 水平向位移。透過公式（18），以雲林 GPS 固定站（位置參考圖 3）的資料計算，推 算本研究該項假設將導致在垂直向約有 1-3 公釐/年的誤差。
= 𝑘𝑣 𝑠 + ( 𝑠 + 𝑠𝑖 )𝑠𝑖 (18) Tk為干涉對的時間基線；v為地表線性變 形速率；θ為衛星視角；φ為衛星軌跡投 影至地面之方位角。
分析水準測量與 TCPInSAR 成果顯示，上述 成果皆具有相同的下陷趨勢，雲林嚴重下陷區整體 呈現碗狀的下陷，主要下陷中心的位置分佈在虎尾 鎮、土庫鎮、元長鄉與褒忠鄉，另外在水準測量無 法 到 達 的 區 域 ， 例 如 麥 寮 濱 海 工 業 區 ， 由
TCPInSAR 成果亦呈現南北差異沉陷之現象。
為分析上述下陷中心的垂直變形狀況，因此 對下陷中心取南北向與東西向的剖面線如圖 13。
其中圖 13（A）為剖面線的位置，圖 13（B）為南 北向剖面資料，由圖中顯示從二崙鄉往南進入虎尾 鎮之後，地層快速下陷，下陷速率由 2 公分/年增 加為 7 公分/年，其中在中科虎尾園區與高鐵虎尾 特定區附近下陷較為顯著，而土庫地區下陷速率則 在 5-7 公分/年之間，往南至元長鄉下陷逐漸減緩。
造成虎尾地區快速下陷的原因，主要在於虎尾 地區富含有薄層的易壓縮土壤（參考本文第一節說 明），依據 Terzaghi et al. （1996）的理論顯示，
土壤完成 90%的壓縮時間（）如公式（19）：
= 4 2 (19) 其中 V土壤壓縮係數，H 為易壓縮土壤厚度，TV為 常數。
因此一樣厚度的土層（相同的 CV值），若土 層中富含易壓縮土壤的厚度（H）越薄，則壓縮時 間（ ）越短，代表下陷速率越快，因此虎尾地區 的地質條件因素也是影響下陷速率的主要因素之 一。
圖 13（C）為東西向剖面資料，由圖中顯示雲 林縣從東勢鄉往東下陷逐漸增加，至土庫地區下陷 最為嚴重，年下陷速率達 6 ~ 7 公分之間，經過土 庫鎮後，下陷趨勢逐漸減緩。
5. 結論與未來工作
透過 TCPInSAR 成果與水準測量成果比對顯 示，TCPInSAR 技術在大範圍地表變形監測，可獲 得比水準測量更多點位的下陷趨勢，同時在水準檢 測網沒有涵蓋的雲林沿海工業區亦測得明顯的變 形趨勢，凸顯出 InSAR 在高空間解析度與大範圍 量測的優勢。
在精度比較方面，經由 268 個水準點成果與 TCPInSAR 進行比對，有 70%點位的下陷差異量在
±1 公分/年以內，兩者的均方根誤差為 0.6 公分/年，
顯示 TCPInSAR 技術具有高精度監測之地表變形
的能力。
由於台灣位於地殼變動與地層下陷頻繁的地 區，目前台灣地區已設置超過 400 個 GPS 固定站，
未來可利用 InSAR 技術進行大範圍的監測，並以 GPS 固定站作為控制點，修正水平向位移，同時融 合 GPS 與 InSAR 的變形監測數據，提高地表變動 區域的監測精度與可靠度，則將來可以獲得全台精 確的地表變形數據，同時降低地面監測的成本。
圖 13 雲林嚴重下陷地區東西剖面資料
參考文獻
工業技術研究院（2011），「多重感應器應用於台北、
彰化及雲林地區地層下陷監測與機制探討計 畫」，經濟部水利署。
洪偉嘉（2009），應用多重感應器監測雲林地區三 維變形。國立交通大學土木工程學系博士論 文。
陳文福、江崇榮（1999），「濁水溪扇州及鄰近地區 之沉積物分布與沉積環境」，地質第十八卷第 二期。
雲林科技大學（2009），「雲林地區多元水資源開發 需求與可行性評估」，第 2-34 頁，經濟部水利 署。
賴典章、費立沅、江崇榮（2003），台灣地下水分 區特性。水文地質調查與運用研討會論文集，
第 1-24 頁。
Abidin, H. Z., Andreas, H., Djaja, R., Darmawan, D.,
& Gamal, M. （ 2008 ） . Land subsidence characteristics of Jakarta between 1997 and 2005, as estimated using GPS surveys. GPS Solutions, 12（1）, 23-32.
Bamler, R. （ 2000 ） . Interferometric stereo radargrammetry: absolute height determination from ERSENVISAT interferograms, in International Geoscience and Remote Sensing Symposium, IGRASS'00, 742-745.
Bombrun, L., Gay, M., Trouvé, E., Vasile, G., & Mars, J. （2009）. DEM Error Retrieval by Analyzing Time Series of Differential Interferograms. IEEE Geoscience and Remote Sensing Letters, 6（4）, 830-834..
Colesanti, C., Ferretti, A., Novali, F., Prati, C., &
Rocca, F. （ 2003 ） . SAR monitoring of progressive and seasonal ground deformation using the permanent scatterers technique. IEEE Transactions on Geoscience and Remote Sensing, 41（7）, 1685-1701.
Ferretti, A., Prati, C., & Rocca, F. （ 2001 ） . Permanent scatterers in SAR interferometry.
IEEE Transactions on Geoscience and Remote Sensing, 39（1）, 8-20.
Gabrysch, R., & Ronald, J. （2005）. Measuring a century of subsidence in the Houston–Galveston region, Texas, USA. Paper presented at the Proc., 7th International Symposium on Land Subsidence, Shanghai （CN）.
Galloway, D. L., & Hoffmann, J. （2007）. The application of satellite differential SAR interferometry-derived ground displacements in hydrogeology. Hydrogeology Journal, 15（1）, 133-154.
Galloway, D. L., Hudnut, K., Ingebritsen, S., Phillips, S., Peltzer, G., Rogez, F., & Rosen, P. （1998）.
Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope Valley, Mojave Desert, California. Water Resources Research, 34（10）,
（C）
（B）
（A）
2573-2585.
Golub, G.H., Kahan, W. （1965）. Calculating the singular values and pseudoinverse of a matrix, SIAM J. Numer. Anal., 2, 205-224.
Hoffmann, J. （2003）. The application of satellite radar interferometry to the study of land subsidence over developed aquifer systems.
Hooper, A., Segall, P., & Zebker, H. （ 2007 ） . Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos.
Journal of Geophysical Research: Solid Earth
（1978–2012）, 112（B7）.
Hooper, A., Zebker, H., Segall, P., & Kampes, B.
（ 2004 ） . A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers.
Geophysical research letters, 31（23）, 611-615.
Hung, W. C., Hwang, C., Chang, C. P., Yen, J. Y., Liu, C. H., & Yang, W. H. （2010）. Monitoring severe subsidence in Taiwan by multi-sensors:
Yunlin, the southern Choushui River Alluvial Fan. Earth Science Geology, 59, 1535-1548.
Hung, W. C., Hwang, C., Chen, Y. A., Chang, C. P., Yen, J. Y., Hooper, A., & Yang, C. Y. （2011）.
Surface deformation from persistent scatterers SAR interferometry and fusion with leveling data: A case study over the Choushui River Alluvial Fan, Taiwan. Remote Sensing of Environment, 115（4）, 957-967.
Hung, W. C., Hwang, C., Liou, J. C., Lin, Y. S., &
Yang, H. L. （2012）. Modeling aquifer-system compaction and predicting land subsidence in central Taiwan. Engineering Geology, 147-148, 78-90.
Massonnet, D., Holzer, T., & Vadon, H. （1997）.
Land subsidence caused by the East Mesa geothermal field, California, observed using SAR interferometry. Geophysical research letters, 24（8）, 901-904.
Terzaghi, K., Peck, R. B., & Mesri, G. （1996）, Soil Mechanics in Engineering Practice, Wiley, New York, 592 p.
Tosi, L., Teatini, P., Carbognin, L., & Frankenfield, J.
（ 2007 ） . A new project to monitor land subsidence in the northern Venice coastland
（Italy）. Environmental Geology, 52（5）, 889-898.
Zebker, H. A., Rosen, P. A. & Hensly S. （1997）.
Atmospheric Artifacts in Interferometric SAR Surface Deformation and Topographic Maps. J.
Geophys Res., 102 （4）, 7547-7563.
Zhang, L., Ding, X., & Lu, Z. （2011）. Modeling PSInSAR time series without phase unwrapping.
Geoscience and Remote Sensing, IEEE Transactions on, 49（1）, 547-556.
Zhang, L., Lu, Z., Ding, X., Jung, H. S., Feng, G., &
Lee, C. W. （2012）. Mapping ground surface deformation using temporarily coherent point SAR interferometry: Application to Los Angeles Basin. Remote Sensing of Environment, 117, 429-439.
1 Adjunct Assistant Professor, Department of Civil Engineering, Received Date: Oct. 10, 2013
2 National Chiao Tung University Revised Date: Jul. 18, 2014
2 Master, Department of Civil Engineering, National Chiao Tung University Accepted Date: Dec. 24, 2014
3 Chair Professor, Department of Civil Engineering, National Chiao Tung University
4 Assistant Professor, Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University
*.Corresponding Author, Phone: 886- 3-5910279, E-mail: [email protected]
Monitoring Serious Land Subsidence in Yunlin Using Temporarily Coherent Point SAR Interferometry
Wei-Chia Hung 1* Da-Ren Huang 2 Cheinway Hwang 3 Lei Zhang 4
ABSTRACT
Yunlin area is suffering severe subsidence caused by groundwater withdrawal and the high compressibility of sediment. The conventional geodetic methods including GPS, leveling and multi-layer compaction monitoring well have been used to monitor land subsidence. Those sensors can only deliver point-wise displacements and the monitoring capacity of displacement is limited by point density. In this study, we applied temporarily coherent point SAR interferometry （TCPInSAR） to monitor land deformation rate over Yunlin area. The proposed method can estimate deformation parameters without the effect of unresolved phase ambiguity resolution.
The study utilizes 15 ALOS PALSAR acquisitions from March 2007 to March 2011 to derive land deformation. The TCP density over Yunlin is 196 pixels/km2, compared to 0.22 point/km2 of the leveling benchmarks. TCPInSAR yields vertical displacements matching the leveling result to 0.6 cm/year （RMSE）.
This research suggests that TCPInSAR with ALOS PALSAR data can effectively monitor land subsidence in a large area like Yunlin area.
