國立交通大學
土木工程學系
碩 士 論 文
以時域相關點雷達干涉量測研究彰化、雲林與嘉
義地區之地層下陷
Detection of Land Subsidence by Temporarily Coherent Point
SAR Interferometry in Changhua, Yunlin, and Chiayi
研究生：黃大任
以時域相關點雷達干涉量測研究彰化、雲林與嘉義地區之地
層下陷
Detection of Land Subsidence by Temporarily Coherent Point
SAR Interferometry in Changhua, Yunlin, and Chiayi
研究生：黃大任 Student：Da Ren Huang
指導教授：黃金維 Advisor：Dr. Cheinway Hwang
國 立 交 通 大 學 土 木 工 程 學 系
碩 士 論 文
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
July 2013
Hsinchu, Taiwan, Republic of China 中華民國一 O 二年七月
以時域相關點雷達干涉量測研究彰化、雲林與嘉義地區之地層下陷
研究生：黃大任 指導教授：黃金維
國立交通大學土木工程學系
摘要
台灣西部過去因不當抽用地下水造成大範圍的地層下陷情形，傳統的監測方 法採用點狀式的大地監測方法（包括水準測量、GPS），上述方法會因監測密度 的 高 低影 響監 測的 精度 。本 研 究以 時 域相關點雷達 干涉技 術（Temporarily Coherent Point SAR Interferometry, TCPInSAR）對地表進行監測，該方法不僅可獲取大範圍地表變形資訊，同時因為不需經由相位解纏（Phase Unwrapping），避 免可能存在相位模糊（Phase Ambiguities）的錯誤解算。 本研究使用了 ALOS 軌道 446 從 2006 至 2011 年 19 幅衛星影像獲取彰化地 區 TCP 密度約為 361 pixel/km2，高於水準點密度 0.2 point/km2，以及軌道 447 從 2007 至 2011 年 15 幅衛星影像獲取雲林與嘉義地區 TCP 密度分別為 196、250 pixel/km2，而該區域水準點密度分別為 0.3 與 0.2 point/km2。本研究將所獲得衛
星視距方向（Line of Sight, LOS）變化量，結合 GPS 資料簡單內插模擬區域內水 平速度場影響後，予以扣除水平向變化，轉換成垂直向變化，比對 TCP 成果與 水準量測成果差異之均方根誤差（Root Mean Square Error）分別為 0.9、0.6、0.8 cm/year，顯示該方法可有效地應用在地表的變形監測。
II
Detection of Land Subsidence by Temporarily Coherent Point SAR
Interferometry in Changhua, Yunlin, and Chiayi
Student：Da Ren Huang Advisor：Dr. Cheinway Hwang
Department of Civil Engineering
National Chiao Tung University
Abstract
Western Taiwan is suffering severe subsidence caused by groundwater withdrawal. The conventional geodetic sensors, including leveling and GPS, 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 derive long-term deformation rates over Changhua, Yunlin, and Chiayi. The proposed method can estimate deformation parameters without the effect of unresolved phase ambiguity resolut.
The study uses 19 ALOS PALSAR images along track 446 from 2006 to 2011 to derive line of sight (LOS) deformation in Changhua. The pixel density over Changhua is 361 pixel/km2, compared to 0.2 point/km2 of the leveling benchmarks. The study also utilized 15 ALOS PALSAR images at 447 track from 2007 to 2011 derive LOS deformation in Yunlin and Chiayi. The pixel density over Yunlin and Chiayi are 196 and 250 pixel/km2, compared to 0.3 and 0.2 point/km2 of the leveling benchmarks respectively. We use GPS data to reduce the effect of the horizontal displacement on the LOS deformation to generate vertical displacement. The vertical displacements from TCPInSAR matches the leveling result to 0.9, 0.6, and 0.8 cm/year (RMSE) in Changhua, Yunlin, and Chiayi. This research suggests that TCPInSAR with ALOS PALSAR data can effectively monitor land deformation in a large area like western Taiwan.
誌 謝
論文最後一頁、最後一夜，感到無比放鬆和愉悅，碩班兩年的充實，或喜或 悲，肯定是此生難以抹去的一段，同時，也將進入人生下個階段。 研究的過程，幸運地得到許多人的幫助，不管是研究相關、生活調劑，還是 精神支持，只有萬分感謝，單憑一己之力，不可能完成這本論文，要感謝的人太 多，若有遺漏還請見諒。感謝指導老師黃金維教授，放手讓我選擇感興趣的方向， 並在觀念、應用上的諸多提點與討論可能性，使我在研究時得到各種啟發，從老 師身上看到研究態度與精神，策勵著我向前；感謝洪偉嘉學長，領著我進入 InSAR 領域，從理論的教導與技術討論中，我才能一點一滴收穫臻至完整，更在待人處 世分享自身經驗，於公於私我獲益良多；感謝中央大學張中白教授，研究上給我 各種意見，從地球科學的角度開拓我的視野，讓我清楚僅在數字上精進無法完整 研究，背後的現象與解釋也是重要環節；感謝東華大學顏君毅教授，不厭其煩的 幫助我解決技術性問題，也在研究討論中給予諸多寶貴意見與經驗；感謝香港理 工大學張磊教授，於研究技術上傾全力無私教導，使我順利完成研究，後續更樂 於討論與提點成果中的任何問題；感謝中央大學陳怡安學長，讓我從毫無悉知的 門外漢進入狀況，也耐心地解決我許多的問題。感謝文化大學陳柔妃教授，豐富 的地質故事讓我了解研究不應缺乏解釋的部分，避免最後僅得到數字成果的空洞； 感謝交大史天元教授與張智安教授認真的教授課程與解惑；工研院地陷團隊的上 智、小光、秋香、怡姿，各種資源提供與打氣支持著我，謝謝各位。 這兩年間，也謝謝研究室的各位，鄭博、綉雯、訓訓哥、強哥、亨利、彥杕、 詠升、邦和、宜珊、雅琦、紫漪、元旎、健輝、俊銘、鎔壑、琨原、書涵，各種 歡樂活動調劑苦悶的研究生活；同學們，朔哥、葉董、YINO、胖、鬼鬼，不管IV
目錄
摘要 ... I Abstract ... II 致謝 ... III 目錄 ... IV 圖目錄... VII 表目錄... X 第一章 緒論 ... 1 1.1 研究動機與目的 ... 1 1.2 文獻回顧 ... 1 1.3 論文大綱 ... 6 第二章 研究區域概論... 8 2.1 彰化、雲林地區 ... 8 2.2 嘉義地區 ... 14 第三章 原理與方法介紹 ... 18 3.1 雷達 ... 18 3.1.1 雷達發展 ... 19 3.1.1.1 航空側視雷達系統 ... 19 3.1.1.2 合成孔徑雷達 ... 20 3.1.2 雷達訊號特性... 21 3.1.2.1 偏極模式 ... 21 3.1.2.2 入射角與觀測幾何 ... 223.1.2.3 雜訊與斑駁現象 ... 24 3.1.2.4 目標物特徵 ... 24 3.2 InSAR 技術 ... 25 3.3 DInSAR 技術 ... 26 3.3.1 二、三、四軌跡差分干涉法... 29 3.3.2 DInDAR 誤差來源 ... 31 3.4 PSInSAR 技術 ... 35 3.4.1 PS 選點 ... 37 3.4.2 相位解算 ... 40 3.4.3 SBAS ... 41 3.5 TCPInSAR ... 42 3.5.1 TCP 選取 ... 42 3.5.2 TCP 網型 ... 45 3.5.3 TCP 參數估計 ... 46 3.5.3.1 先驗精度 ... 48 3.5.3.2 初步參數估計 ... 49 3.5.3.3 剔除相位模糊 ... 49 3.5.3.4 方差-協方差矩陣估計 ... 50 3.5.3.5 參數估計 ... 52 第四章 雷達干涉技術研究流程與成果分析... 53
VI 4.4 TCP 成果分析 ... 56 4.4.1 彰化地區 ... 56 4.4.2 雲林地區 ... 61 4.4.3 嘉義地區 ... 67 4.5 TCPInSAR 成果分析 ... 71 4.5.1 衛星波長 ... 71 4.5.2 衛星影像數量... 72 4.5.3 水準測量施測時間 ... 72 4.5.4 水平向位移 ... 72 4.5.5 視角誤差 ... 72 第五章 結論與建議 ... 74 參考文獻 ... 77 附錄 ... 84
圖目錄
圖 1.2.1 雲林地區 1996 年至 1999 年地表形變速率(LOS 向)(Tung et al.,2012) ... 5 圖 1.2.2 濁水溪沖積扇 2006 年至 2008 年 PSInSAR 地表垂直下陷速率(Hung et al., 2011) ... 5 圖 1.2.3 濁水溪沖積扇 2006 年至 2008 年地表垂直下陷速率 PSInSAR 與水準測
量 Fusion 成果(Hung et al., 2011) ... 6 圖 2.1.1 彰雲地區之區域地質圖(中央地質調查所，2000) ... 8 圖 2.1.2 彰化地區民國 81 年至 100 年累積下陷量圖(工業技術研究院，2011) .... 9 圖 2.1.3 彰化地區民國 81 年至 90 年累積下陷量圖(工業技術研究院，2011) ... 9 圖 2.1.4 彰化地區民國 90 年至 100 年累積下陷量圖(工業技術研究院，2011) .. 10 圖 2.1.5 彰化地區民國 96 年至 100 年平均下陷速率圖(工業技術研究院，2011) ... 10 圖 2.1.6 彰化地區民國 98-99 年與 99-100 年平均下陷速率等值線比對圖(工業技 術研究院，2011) ... 11 圖 2.1.7 雲林地區民國 81 年至 100 年累積下陷量圖(工業技術研究院，2011) .. 12 圖 2.1.8 雲林地區民國 81 年至 88 年累積下陷量圖(工業技術研究院，2011) .... 12 圖 2.1.9 雲林地區民國 88 年至 100 年累積下陷量圖(工業技術研究院，2011) .. 13 圖 2.1.10 雲林地區民國 96 年至 100 年平均下陷速率圖... 13 圖 2.1.11 雲林地區民國 98-99 年與 99-100 年平均下陷速率等值線圖(工業技術研 究院，2011) ... 14 圖 2.2.1 嘉義地區區域地質圖(中央地質調查所，2000) ... 14 圖 2.2.2 嘉義地區民國 80 至 100 年累積下陷量圖(工業技術研究院，2011) ... 15
VIII 術研究院，2011) ... 17 圖 3.1.1 航空側視雷達示意圖... 20 圖 3.1.2 合成孔徑陣列示意圖... 21 圖 3.1.3 合成陣列天線示意圖(陳卉瑄，2001) ... 21 圖 3.1.4 入射角與地形效應對雷達回波強度影響示意圖... 22 圖 3.1.5 側視雷達觀測幾何... 23 圖 3.1.6 雷達波與相異目標物特徵交互作用... 24 圖 3.2.1InSAR 成像幾何示意圖(修改自謝嘉聲，2006) ... 25 圖 3.3.1DInSAR 幾何示意圖(修改自洪偉嘉，2009) ... 27 圖 3.3.2 大氣壓力與濕度變化對不同波段之相位計算影響(Zebker et al., 1997) .. 32
圖 3.3.3L 波段相對溼度與大氣延遲造成 DInSAR 誤差示意圖(Zebker et al., 1997) ... 32 圖 3.3.4DInSAR 成果軌道誤差示意圖... 33 圖 3.4.1 加州長谷火山口地區平均地表位移量(Hooper et al., 2004) ... 36 圖 3.4.2 加州長谷火山口地區於時間序列之地表變形量(Hooper et al., 2004) ... 36 圖 3.4.3 永久散射體雷達回波示意圖(Hooper et al., 2007) ... 37 圖 3.4.4 相位解纏示意圖... 41
圖 3.5.1 香港機場區域計算影像偏移量之示意圖(Zhang et al., 2011a) ... 43
圖 3.5.2 罩窗尺寸與像素相關性改變在估計偏移量標準差的影響 (Zhang , 2012b) ... 44
圖 3.5.3 罩窗尺寸對像素相關性與偏移量計算影響(Zhang, 2012b) ... 45
圖 3.5.4Delaunay 三角網建立方法(Zhang et al., 2011b) ... 46
圖 3.5.5 殘差離群檢測(Zhang et al., 2012a) ... 50
圖 4.1.1TCPInSAR 解算流程 ... 53
圖 4.4.1 彰化地區影像對空間、時間基線資訊... 57
圖 4.4.3 彰化地區 TCPInSAR 原始成果(LOS 向) ... 58 圖 4.4.4 彰化地區 TCPInSAR 成果 ... 59 圖 4.4.5 彰化地區 2007 至 2011 年水準測量平均下陷速率成果 ... 59 圖 4.4.6 彰化地區 TCPInSAR 與水準測量成果差異量分布直方圖 ... 60 圖 4.4.7 彰化地區 TCPInSAR 成果與水準測量成果差異點位 ... 60 圖 4.4.8 雲林地區影像對空間、時間基線資訊... 62 圖 4.4.9 雲林地區 TCP 選取示意圖 ... 62 圖 4.4.10 雲林地區 TCPInSAR 原始成果(LOS 向) ... 63 圖 4.4.11 雲林地區西側 TCPInSAR 成果 ... 64 圖 4.4.12 雲林地區 2007 至 2011 年水準測量平均下陷速率成果 ... 64 圖 4.4.13 雲林地區 TCPInSAR 與水準測量成果差異量分布直方圖 ... 65 圖 4.4.14 雲林地區西側 TCPInSAR 成果與水準測量成果差異點位 ... 65 圖 4.4.15 濁水溪沖積扇示意圖(Hung et al., 2010) ... 66 圖 4.4.16 濁水溪沖積扇剖面地質圖(剖線位置如圖 4.4.15)(Hung et al, 2010) ... 66 圖 4.4.17 嘉義地區影像對空間、時間基線資訊... 67 圖 4.4.18 嘉義地區 TCP 選取示意圖 ... 68 圖 4.4.19 嘉義地區 TCPInSAR 原始成果(LOS 向) ... 69 圖 4.4.20 嘉義地區 TCPInSAR 成果 ... 69 圖 4.4.21 嘉義地區 2007 至 2011 年水準測量平均下陷速率成果 ... 70 圖 4.4.22 嘉義地區 TCPInSAR 與水準測量成果差異量分布直方圖 ... 70 圖 4.4.23 嘉義地區西側 TCPInSAR 成果與水準測量成果差異點位 ... 71
X
表目錄
表 3.1.1 雷達常用波段資料... 18 表 4.2.1ALOS 影像列表 ... 54 表 4.4.1 彰化地區解算參數說明... 56 表 4.4.2 雲林地區解算參數說明... 61 表 4.4.3 嘉義地區使用參數說明... 67 表 5.1 取樣密度－與水準測量比較... 74 表 5.2 研究區域下陷分析... 75 表 5.3 水平位移改正... 75第一章 緒論
1.1 研究動機與目的 台灣西部地區因各產業發展，用水需求大幅提升，故抽取地下水因應地表水 資源不足與降低用水成本，導致許多地區地下水過度抽取而引起地層下陷災害。 而隨著台灣高速鐵路的興建與通車，因其沿線穿越台灣西部平原地層下陷嚴重區 域，包含雲林虎尾鎮與元長鄉地區，成為高鐵行車安全的一大隱憂，使地層下陷 的嚴重性受到重視。 傳統使用的地層下陷監測方法仍以點狀式的大地監測為主，包含 GPS 與水 準測量，上述點狀測量會受到點位密度高低影響檢測的精度，且增加監測點位時 又需增加大量的人力、時間與金錢成本。近年來衛載 SAR 應用於地表變形監測 的技術逐漸成熟，其優勢為可在短時間獲取大範圍之地表變形資訊，且技術演進 使精度可提升至公分甚至釐米等級，故使用 SAR 影像來監測大範圍的地層下陷 成為趨勢。 目前台灣西部平原的嚴重下陷區分布於彰化、雲林、嘉義、台南與屏東地區， 多屬於農業或養殖漁業重要發展地區，因超抽地下水而有大範圍的地層下陷，故 本研究使用 ALOS 衛星影像，以彰化、雲林與嘉義地區進行時域相關點雷達干 涉技術(Temporarily Coherent Point SAR Interferometry, TCPInSAR)對地表變形進 行監測，再嘗試於資料後處理時加入 GPS 資料以消除水平場影響量，以水準測 量成果進行比對，評估 ALOS 影像以 TCPInSAR 技術解算用於彰雲嘉地層下陷 監測的可行性。1.2 文獻回顧
有對月球的觀測研究(Zisk,1972)。Graham(1974)在傳統的 SAR 系統上加裝一橫軌 方向的天線，利用兩個天線同時獲得的相位值形成干涉效果，以類比的方法計算 得到地面高程值。 隨著雷達影像資料量增大，美國噴射推進實驗室利用數位方式記錄雷達影像 的震幅(Amplitude)與相位(Phase)資訊，其包含了每一點的複數資料，可直接計算 地面的高程資訊，此研究使用側視空載 SAR 系統成功的得到了解析度 10 公尺的 數值地形模型(Zebker and Goldstein,1986)。
Gabriel et al. (1989)證實了利用 DInSAR 技術能夠獲得公分級精度的變形量， 其使用三張 Seasat 影像以 DInSAR 技術獲得兩幅干涉條紋圖(interferogram)，並 得到地表變形量。Massonnet et al. (1995)使用 ERS 於地震前後的兩幅影像產生干 涉圖，再使用現有 DEM 資料扣除地形效應後，得到僅受地震影響的干涉圖，成 功的偵測出 1992 年於美國加州 Lander 地區的 7.3 級地震產生的地表變形量。而 Zebker et al. (1994)則使用三幅 ERS 影像產生兩幅干涉圖後進行差分計算其變形 量，該法也成功的偵測出 Lander 地震造成的地表變形量，並使用地真資料進行 比對，也得到同樣為公分級的精度。後續有更多關於地震變形的研究被發表，也 證實了此技術應用於震後形變的實用性。
DInSAR 技術除了用於地震後形變監測外，也用於火山變形監測(Massonnet et al.,1995；Lu et al.,1998)；而 Joughin et al. (1995) and Rignot et al. (1995)則利用 DInSAR 技術作冰河漂移監測；以及 Massonnet et al. (1997)監測加州 East Mesa 地熱區的地層下陷(Subsidence)，且因地層下陷問題為全球許多地區遭遇的問題， 陸續有相關研究應用 DInSAR 於地層下陷監測 (Galloway et al.,1998；Amelung et al.,1999；Hoffmann et al., 2001)。
然而 DInSAR 較適用於平原都會地區，在山區與植被較厚的地區，因地表特 徵物會隨時間變化及受到大氣效應影響，使得干涉成果的精確性較差。而雷達波 雖然可以穿透大氣中的雲、霧及水氣，但可能使雷達波傳遞遲滯，進而造成計算 誤差(Pathier, 2003)。
為改善 DInSAR 技術的限制，Ferretti et al. (2000)提出永久散射體(Permanent Scatterer, PS)的概念與 PSInSAR 處理程序，PSC(PS Candidate)選取是基於多張 SAR 影像的強度值(Amplitude)與其變異程度，強度值高的像素在影像上較亮且 容易辨識，當多張影像中像素的標準差與像素強度平均值之比值小於一門檻時， 則將該像素視為 PSC。基於像素強度值進行 PSC 挑選需要多於 30 張影像(Ferretti et al.,2001)才具有統計意義，最後由 PSC 中萃取出 PS 點並得到 PS 點的形變速率、 DEM 誤差與大氣效應誤差。而後多位學者改善並發展出永久散射體干涉技術 (Persistent Scatters SAR Interferometry, PSI)，Hooper et al. (2004) and Kampes (2006) 根據相位的穩定性來辨識 PSI 點，能夠找出具有較低振幅但變化穩定之 PS 點。 PSInSAR 將原先 DInSAR 處理中，因 DEM、大氣遲滯效應與地表物非同調性引 起的誤差降低，進而提高干涉測量成果之精確度與可信度並獲得時間序列上的地 表變形量。
PSInSAR 在進行處理時，相位解纏(phase unwrapping)是一項關鍵的工作， 過去主要採用兩種方法，一種是以經典的離散最小費用流算法(sparse MCF)對每 一幅干涉圖解纏。因相關點無法被精確選定，在某些相關點上可能含有較大的相 位噪聲，使得這些點的整數模糊度解算失敗，並導致整體的解纏誤差，對於大區 域一般採用分塊(patches)處理，但塊與塊之間的相位跳變容易估算錯誤，造成解 纏後的相位不連續，即為解纏誤差；另一常用的方法是構建一時域相關指標 (temporal coherence)，透過解空間搜索(solution space searching)的方式從纏繞相位 中獲取變形參數與 DEM 誤差，但在參數搜索過程中會出現多組參數(即形變與 DEM 誤差)使得參數解不唯一，屬於對整數模糊度(phase ambiguity)的錯估。為了 解決上述問題，Zhang et al. (2011a；2011b；2012a；2012b)提出了 TCPInSAR (Temporarily Coherence Point SAR Interferometry)技術，使用最小二乘法進行參數
同調性時會因罩窗尺寸改變，此情形於低相關區域更為明顯，改以像素強度值估 算像素偏移量標準差，可在低振幅區域找到更多偏移量變化小且一致的 TCP，亦 與永久散射體的概念相似。
台灣亦有許多 InSAR 的相關研究，近年來對於台灣西部彰化、雲林與嘉義 地區的研究則包含 DInSAR 和 PSInSAR 相關技術。盧玉芳(2007)以 PSInSAR 對 雲林縣土庫與元長兩鄉鎮進行監測，發現於 1996 至 1999 年，此地區最大下陷可 達 7 cm/year，與水準測量相比精度應可達公分等級。童忻(2008)同樣研究 1996 至 1999 年雲林地區於內陸地區最大下陷量達 6 cm/year，沿海地區則有 1.58 cm/yaer 下陷。Hung et al. (2010)則應用了 DInSAR 監測雲林地區的地層下陷情形， 成果與水準測量之 RMSE 約 1 至 2 cm/year，並使用 Fusion 法結合水準測量與 DInSAR 資料，使 RMSE 提升至 0.7 至 0.9 cm/year。而 Hung et al. (2011)亦使用 PSInSAR 技術對濁水溪沖積扇進行監測，於彰化和雲林地區都觀察到最大下陷 7 cm/year，甚至於六輕工業區部分可達下陷速率 8cm/year，與水準測量的 RMSE 為 0.6 cm/year(圖 1.2.2)，使用 Fusion 結合水準測量與 PSInSAR 成果後，RMSE 提升至 0.4 cm/year(圖 1.2.3)。楊佳祥(2011)提出了改良式 PSInSAR 技術，使用 短基線進行影像配對，並測試不同同調性門檻，得到較多的 PSC 點位，而解算 成果精密度 RMS 為 1.52 cm/year。蘇柏宗(2012)則使用多組不同主影像的
PSInSAR 成果，以多餘觀測的概念解算出一組平均成果，成果精度 RMS 由 1.53 cm/year 提升至 0.71 cm/year。
圖 1.2.3 濁水溪沖積扇 2006 年至 2008 年地表垂直下陷速率 PSInSAR 與水準測 量 Fusion 成果(Hung et al., 2011)
1.3 論文大綱 第一章：緒論 包含研究動機與目的、文獻回顧與本文章節架構說明。 第二章：研究區域概論 本章分為兩個小節，分別簡述濁水機沖積扇與嘉義地區過去與近年來的地層 下陷情形。 第三章：原理與方法 本章分為五個小節，說明了雷達系統演進、InSAR 技術、DInSAR 技術、 PSInSAR 技術與 TCPInSAR 技術。 第四章：雷達干涉技術研究流程與成果分析
本章分為五個小節，第一小節說明本研究流程、第二小節說明研究使用的資 料來源、第三小節說明初步成果後處理的方式、第四小節說明各研究區域的成果 與地真資料比較、第五小節用於分析討論成果。
第二章 研究區域概論
2.1 彰化、雲林地區 彰化與雲林地區位於台灣西部平原中段，平原區域為濁水溪沖積而成，又稱 為濁水溪平原，其地面地質大部分為現代沖積層(如圖 2.1.1)，砂礫質沉積物為 河道狀分布，由扇頂向西呈輻射展開；扇央與扇尾部分顆粒較細，含泥質較多， 其地層多半為細砂與黏土交錯層組成(工業技術研究院，2011)。 彰化地區過去因不當抽取地下水造成大範圍的地層下陷，由水準測量成果可 發現，早期主要下陷區域集中於沿海地區的大城鄉，而內陸於芳苑鄉、二林鎮、 竹塘鄉、埔鹽鄉、西湖鎮、埤頭鄉與溪州鄉亦有下陷(圖 2.1.2)，然民國 90 年後， 彰化地區的下陷中心向內陸移動，故將彰化地區民國 81 年至 100 年之累積下陷 量分為民國 81 年至 90 年(圖 2.1.3)與民國 90 年至 100 年(圖 2.1.4)兩階段觀察， 發現於民國 81 年至 90 年間，下陷區域集中於大城鄉，民國 90 年至 100 年時大 城鄉下陷情況趨緩，而內陸的二林鎮、溪湖鎮、溪州鄉與埤頭鄉因產業發展增加 抽水量，造成下陷情況。圖 2.1.5 為彰化地區民國 96 年至 100 年的年平均下陷 速率圖，可看出過去四年間彰化地區地下陷區域集中於二林鎮、溪湖鎮與溪州鄉 三個區域。彰化地區各年度的最大下陷速率與嚴重下陷面積數據參照附錄 一。 近年彰化地區下陷情形，由圖 2.1.6 民國 98 年至 99 年與 99 年至 100 年的 平均下陷速率等值圖，觀察出下陷速率大於 3 公分的面積大幅減少，顯見彰化地 區的地層下陷情形逐漸趨緩。 圖 2.1.1 彰雲地區之區域地質圖(中央地質調查所，2000)圖 2.1.4 彰化地區民國 90 年至 100 年累積下陷量圖(工業技術研究院，2011)
圖 2.1.6 彰化地區民國 98-99 年與 99-100 年平均下陷速率等值線比對圖(工業技 術研究院，2011) 雲林地區為台灣重要的農業區之一，因民國 70 年代農產品價格欠佳收益不 良，加上養殖業利潤優厚之利誘下，導致沿海地區農地大量變更為漁業養殖。惟 養殖事業須仰賴大量淡水以保持魚池水質潔淨，在沿海地區地面水源缺乏的情況 之下，養殖業者轉而抽汲地下水(工業技術研究院，2011)。 根據水準測量成果，雲林地區早期於民國 80 年至 88 年，下陷中心為沿海的 台西鄉、麥寮鄉與內陸的元長鄉(圖 2.1.7)，而後雲林下陷中心逐漸向內陸移動， 從民國 88 年至 100 年的累積下陷圖(圖 2.1.9)可觀察出下陷中心集中於褒忠鄉、 土庫鎮、虎尾鎮與元長鄉，沿海地鄉鎮的下陷情形趨緩。圖 2.1.10 則為民國 96 年至 100 年年平均下陷速率圖，其嚴重下陷區與民國 88 年至 100 年的累積下陷
林地區的地層下陷情形逐漸趨緩。
圖 2.1.7 雲林地區民國 81 年至 100 年累積下陷量圖(工業技術研究院，2011)
圖 2.1.9 雲林地區民國 88 年至 100 年累積下陷量圖(工業技術研究院，2011)
圖 2.1.11 雲林地區民國 98-99 年與 99-100 年平均下陷速率等值線圖(工業技術 研究院，2011) 2.2 嘉義地區 嘉義地區位於台灣本島西岸之南部，地面地質大部分為現代沖積扇(如圖 2.2.1)，區域內的沖積層為河道堆積物、海岸風機砂與潟湖淤泥等物質組成。而 區域內多條河流，於上游多泥質岩層，故沖積層中缺少粗鬆的沉積物，其厚度在 東部邊緣最薄，向西逐步增厚(工業技術研究院，2011)。 圖 2.2.1 嘉義地區區域地質圖(中央地質調查所，2000)
嘉義地區過去的土地利用主要為農田，而後於沿海地區逐漸興起養殖漁業， 然各產業用水隨著產業發展而提升，地下水抽取則為普遍現象，便造成大範圍的 地層下陷情形。圖 2.2.2 為民國 80 年至 100 年的累積下陷量，可看出嘉義地區 20 年來的嚴重下陷區域涵蓋了東石鄉、布袋鎮、朴子市、義竹鄉等鄉鎮，而主 要下陷中心為東石鄉與布袋鎮。而圖 2.2.3 為民國 93 年至 100 年的累積下陷量， 其下陷區域大幅縮小，但仍分布於布袋鎮、東石鄉與義竹鄉之部分區域。嘉義地 區逐年最大下陷速率與嚴重下陷面積參考附錄 三。 而近年來嘉義地區的下陷情形，由圖 2.2.4 與圖 2.2.5 比較民國 96 年至 100 年平均下陷速率圖、民國 98 年至 99 年與 99 年至 100 年的平均下陷速率等值圖， 可觀察出下陷速率於 3 公分以上的面積大幅減少，顯見嘉義地區的地層下陷情形 逐漸趨緩。 圖 2.2.2 嘉義地區民國 80 至 100 年累積下陷量圖(工業技術研究院，2011)
圖 2.2.3 嘉義地區民國 93 至 100 年累積下陷量圖(工業技術研究院，2011)
圖 2.2.5 嘉義地區民國 98 至 99 年與 99 至 100 年平均下陷速率等值線圖(工業技 術研究院，2011)
第三章 原理與方法介紹
目前的大地測量法常用的為水準測量和 GPS 測量，兩者皆屬於點的量測， 對於大範圍的量測需耗費大量的人力或資源提升觀測點位密度，又或是使用光達 技術獲取區域內精細的地表資訊，光達技術雖可獲得高密度的成果，但其單次量 測即耗費大量資金，若進行長期的監測是一沉重的負擔。因此目前於大範圍的量 測多使用雷達干涉技術進行，其具有量測範圍廣大、相較上述量測法低成本，且 僅需少數人力配置即可完成。 而雷達干涉技術從合成孔徑雷達的發展後，後續提出了 DInSAR、PSInSAR、 SBAS 與 TCPInSAR 等技術，本章節將分別介紹雷達干涉技術的原理及其後續發 展的技術。 3.1 雷達雷達(Radar)為 radio detection and ranging (無線偵測與距離量測)的縮寫，意 指透過一定向天線對目標物發射一系列雷達波，經無線電波接觸物體產生反射波 再由接收器接收後，量測其反射波強度與時間間隔而獲取方向與距離量測，屬於 主動式系統(active system)。 雷達系統的常用波段如表 3.1.1，各波段的電磁波對物體的反射程度不同， 應依據目標物性質選擇適合的波段，目前雷達衛星使用的波段以 C、L、X 三種 波段最為常見，此類波段屬於長波長訊號，在大氣中具有較佳的穿透性，因此主 動式雷達可以全天候進行量測，並且不受天候影響，部分波長甚至能穿透低矮植 被。 表 3.1.1 雷達常用波段資料 波段 頻率 波長 載台 L 1-2GHz 15-30cm JERS、SEASAT、ALOS C 4-8GHz 3.75-7.5cm RadarSAT-1/2、ERS-1/2、 Envisat X 8-12GHz 2.5-3.75cm SRTM、TerraSAR、 COSMO-SkyMed
3.1.1 雷達發展
3.1.1.1 航空側視雷達系統
航空側視雷達(Side-Looking Airborne Radar system, SLAR)(如圖 3.1.1)運作 上以側視(side looking)進行觀測主要原因有二 : 第一，飛機於行進中接受回波訊 號時，航行方向上的訊號會受到都卜勒效應(doppler effect)影響，造成回波頻率 改變，而在側向接收回波則可以令都卜勒效應造成影響最小化。第二，雷達影像 是藉由接收回波的先後時間組成，使用側視方向接收回波能夠分辨出距離的遠近， 因而提高雷達影像解析力。而 SLAR 的影像解析度又可分為兩種：(1)距離解析 度(range resolution)、(2)軌道方向解析度(azimuth resolution)
(1)距離解析度 距離解析力ΔRg是定義此雷達系統能夠辨識測距方向上兩個目標物的 最小距離，受到發射波段的脈衝持續或脈衝寬度而決定，理論上為脈衝寬度 一半，可用下式表示： 式 3-1 中的 τ 為脈衝發射的持續時間(pulse duration)，單位為微秒(10-6 sec)， c 為光速，θ 為雷達對目標物的俯角(平視時 θ 為 90∘)。正常情況中，雷達 發射器能夠發射的脈衝持續時間τ 很短，造成能量不足以達到合理的訊雜比
(Signal to Noise Ratio，SNR)，因此需藉由脈衝壓縮技術(pulse compression) 獲取較高解析度與 SNR。
(2)軌道方向解析度
軌道方向解析力ΔRa是雷達掃描地表時於軌道方向的寬度，若要於軌道
具其解析度尚可接受，但若使用衛星作為載具，因 r 和 λ 因衛星軌道高度與 波段固定屬於常數無法變動，若想提升軌道方向解析度，則需加大天線長度， 所以在衛載上 SLAR 系統的解析度受到限制，此後則發展出合成孔徑雷達 (Synthetic Aperture Radar, SAR)解決此限制。
圖 3.1.1 航空側視雷達示意圖
3.1.1.2 合成孔徑雷達
合成孔徑雷達(Synthetic Aperture Radar, SAR)是用以改進 SLAR 無法提高解 析度的缺點，當載具航高與使用雷達波段固定時，僅能增加天線長度來提高軌道 方向解析度，但載具的負重與尺寸均有嚴格的限制，因此才發展出在天線長度受 限時仍能提升影像解析度的合成孔徑雷達。 合成孔徑雷達是利用合成天線陣列的想法來提高軌道方向解析度，如圖 3.1.2，透過載具飛行時雷達不斷發射脈衝並接收回波訊號，在多個位置能看見同 一目標物時，因雷達天線至目標物的距離不相等，使得回波產生不同的相位移， 為了獲得正確數據，雷達將這些回波儲存於不同的記憶體位置，再經由相位平移 器給予相反方向的相位平移補償相位差，使該目標物的所有回波相位趨於一致， 最後將不同目標物的樹出成果總合，獲得現有的 SAR 影像，處理流程如圖 3.1.3。
圖 3.1.2 合成孔徑陣列示意圖
圖 3.1.3 合成陣列天線示意圖(陳卉瑄，2001)
3.1.2 雷達訊號特性 3.1.2.1 偏極模式
極，而 ALOS PALSAR 則根據拍攝方式有多偏極與全偏極選項。在觀測物語偏極 方向類似的平面上，雷達回波訊號較弱，如：拍攝垂直的稻作時，HH 偏極模式 所得影像會比 VV 偏極模式還要亮(較亮表示訊號較強)。 3.1.2.2 入射角與觀測幾何 入射角(incidence angle)指目標物位置的垂直方向與雷達波束入射方向的夾 角。回波強弱與雷達波束的入射角有極大相關，若雷達波和目標物接觸面呈直角 或高入射角，回波訊號強度會比低入射角強；以低入射角感測時，雷達波能量反 射回感測器較低，因此在影像上的色調偏暗。另外，目標物的局部入射角(local incidence angle)也是影響雷達回波強度的主因，其定義為目標物所在的當地地形 之垂直方向和雷達波束入射方向的夾角，性質相近的目標物會因所在位置地形不 同、相對於感測器的角度方位不同而使回波的成像強度有所差異(如圖 3.1.4)。 圖 3.1.4 入射角與地形效應對雷達回波強度影響示意圖 而地形坡度的差異，使雷達訊號的相對入射角不同，除了回波強度有差異之 外，此觀測幾何(viewing geometry)也可能對訊號產生幾何變形。因 SAR 屬於側 視感側成像系統，記錄目標物至感測器的斜距(slant range)，但成像依據接收回波 的先後順序而定，因此在地形起伏較大時，容易產生幾何變形，造成成像與實際 情形不符，如圖 3.1.5a，而幾何變形有以下三種：
(1)前波縮短(foreshortening) 當面向坡(foreslope)的坡度小於雷達波束與地面的夾角時(α<θ)，則產生 前波縮短現象，如圖 3.1.5b 所示，因雷達波接觸到面向坡山腳 a 和山頂 b 的時間相近，因此面向坡的回波訊號會在a’b’成像，造成影像中的特徵比實 際情況縮短且亮(a’b’<ab)。 (2)疊置(layover) 當面向坡的坡度大於雷射波束和地面的夾角時(α>θ)，產生疊置現象， 如圖 3.1.5c 所示，因山頂 b 的雷達回波最先被接收，造成背向坡(backslope) 和面向坡的回波訊號在影像上重疊，主動疊置(Active Layover)區 b’’a 以及遠 距(far range)的被動疊置區 ab’’的回波訊號在影像 a’b’上重疊，使疊置區域特 別亮，如：陡峭的山脊、稜線在影像上特別亮。
(3)陰影(shadow)
當背向坡過於陡峭( )，雷達波被遮蔽無法到達，該區將 無回波資料，且除了背向坡(bc)，後方部分區域(cd)也會無資料，在影像上 b’d’ 僅有雜訊，如圖 3.1.5d。
3.1.2.3 雜訊與斑駁現象 雷達成像在同一像元內的不同目標物所反射的雷達回波訊號可能為同相位 (in-phase)，強化了該像元的回波強度，但也可能是因為反相位(out-of-phase)而使 回波強度相互抵消而減弱，因此雷達影像常呈現粒狀(grainy)或木紋狀，難以分 辨特徵，需要使用適當的濾波器凸顯影像特徵，增加其應用價值。 3.1.2.4 目標物特徵 地表物的幾何形狀及其各種物理性質皆會影響雷達波吸收及反射，其地表反 射物可分為三種(如圖 3.1.6)： (1)鏡面反射體(specular reflector) 以水體為例，當雷達波入射至無風雨之平滑如鏡面的水體表面時，雷達 波會做全反射向外散出，導致感測器無法接收到回波訊號。 (2)直角反射體(corner reflector) 直角反射體可以產生最強的雷達回波，因直角反射體是由三個互呈直角 的平面構成，入射的雷達波會經兩次全反射後以相同角度反射回波至感測器， 因此此類的回波訊號較強。 (3)散射體(diffuse reflector) 若目標物表面凹凸不平或地勢崎嶇，入射的雷達波接觸到目標物後會向 周圍散射，其中部分能量會反射回感測器，但回波訊號較微弱且不均勻。以 植被為例，地表植被分布區域容易使雷達波產生散射作用，造成影像的同調 性降低無法產生干涉條紋。 圖 3.1.6 雷達波與相異目標物特徵交互作用
3.2 InSAR 技術
合成孔徑雷達干涉(Synthetic Aperture Radar Interferometry, InSAR)技術為利 用不同時間或不同位置的雷達天線所獲得的兩幅或多幅複數影像進行干涉計算， 藉由相位值差異得到地表三維資訊。為了解其原理，將 InSAR 技術之簡單成像 幾何關係表示如圖 3.2.1，其中 A1、A2為獲取地表資訊時天線的空間位置、B 為 兩天線之基線長、H 為載具航高、θ 為視角、r 為第一個天線至目標物的距離、δ 為第二個天線至目標物距離減去 r 後的距離差值、α 為基線與水平線之夾角。 圖 3.2.1InSAR 成像幾何示意圖(修改自謝嘉聲，2006) 兩天線接收同一目標物訊號相位差值(ψ)，與目標物距離變化(δ)關係為： 因雷達成像是由天線發射雷達波後經目標物再折返至天線接收，訊號傳播路 徑差為兩倍的目標物距離變化(2δ)，代入上式可得到接收訊號時的相位差值： 因此可求出距離差值(δ)為：
經式 3-4 可得到： 移項後得： 經式 3-8 可導出 r 為： 由圖 3.2.1 中，目標物高程為： 並將式 3-5 與式 3-9 代入式 3-10 中得目標物高程值為： 根據上述方法推導至式 3-11，可了解 InSAR 技術只要能確定軌道參數值： 航高、視角、基線長與基線和水平線之夾角等資料，已知雷達波長時，便可計算 目標物的高程值。 3.3 DInSAR 技術
合成孔徑雷達差分干涉(Differential Synthetic Aperture Radar Interferometry, DInSAR)技術係用於量測地表微小變動，利用不同時間、相同區域的兩張 SAR 影像進行干涉，並移除原始地形效應後，即可獲得在干涉影像時間對內純粹由地 表變形所產生的相位干涉圖。但在差分干涉技術中，僅能由相位資訊獲得影像對 時間內的地表變形產生相位變化，須經由相位回復才可得到沿雷達視距方向
(Line of Sight, LOS)之地形變化量，圖 3.3.1 為 DInSAR 技術之成像幾何示意圖：
圖 3.3.1DInSAR 幾何示意圖(修改自洪偉嘉，2009)
圖 3.3.1 中，A1、A2與 A3為三個不同時間天線之位置，同 InSAR 原理，A1
與 A2兩天線所獲得的影像可形成第一組干涉對，假設於此兩張影像的時間間隔 內，目標點由 P 移動至 P’(P、P’與 A1三點共線)，其中，A1至 P 的距離為 r1=r2+δ1、 A2至P’點的距離為 R2+Δr，Δr 為目標物沿 S2 視距方向地表位移所產生之位移向 量。 套合兩張 SAR 影像後，根據 InSAR 原理，S1與 S2兩影像間的 P 與 P’相位 差ψ 為： 式 3-12 中，ψ 為相位回復後的相位差，由圖 3.3.1 可知(Zebker et al.,1994)： 由式 3-12 可觀察出，相位差 ψ 是由地形效應產生的相位差 和地表變形
為了移除因地形效應產生的相位差 ，假設天線 A1與 A3影像對拍攝 期間地表沒有變形，即無地表變形產生的位移向量，則此第二組干涉影像對可獲 得純粹由地形效應產生的相位差： 且由圖 3.3.1 可知： 因此由式 3-13、式 3-15 與式 3-16 可得： 再將式 3-17 代入式 3-14 得： 若 可求得，則地表的每一點地表位移(ΔR)即可求出。由式 3-13 與式 3-16 可得： 式 3-19 中 θ 為視角，θ 是根據每一點的地形與成像幾何而得，若要直接使用 式 3-18 與式 3-19，必須先給予此區域的高程資訊，可以選用現有之處值高程模 型資料(二軌跡差分干涉法，詳見 3.3.1 節)，或其他干涉對所生成的高程資料(三、 四軌跡差分干涉法，詳見 3.3.1 節)。 InSAR 之應用除了地表變形監測之外，亦可建立術值高程模型，以下將比較 干涉成果中，地表變形與地表起伏所產生之相位差靈敏度。 由圖 3.3.1 已知 h=H-rcosθ，其中 h 為目標物高程、H 為載具航高、r 為天線 至目標物距離、θ為視角。 由式 3-10 可得 dθ對 h 之影響 dh 如下： 由式 3-12 與式 3-13 可知：
而由式 3-21 可知 dθ對ψ之影響 dψ如下： 由式 3-20 與式 3-22 可得 dh 對ψ之影響 dψ如下： 由式 3-21 可得到 d(Δr)對ψ之影響 dψ如下： 對衛載的 SAR 而言，航高遠大於基線長(r>>B)，由式 3-23 與式 3-24 可知， 對等量變化的地表起伏(dδ)與地形變化(Δr)，地形變化(Δr)產生之相位變化大 於地表起伏(dδ)產生的相位變化，表示 DInSAR 技術量測地形變化(Δr)造成的 相位變化較地表起伏(dδ)相位變化敏感，所以利用衛載的 DInSAR 技術製作術 值高程模型時，精度為公尺等級；而應用於監測地形變化時，精度可達公分等級。 3.3.1 二、三、四軌跡差分干涉法 DInSAR 需要一組干涉成果作為地形對，代表地表高程，另一組干涉成果作 為變形對，而變形對資料包括了地形效應與地形變化，將兩組干涉成果進行差芬 處理後，可移除變形隊中的地形效應，僅保留因地形變化造成的相位差值。而地 形對資料依據使用資料來源不同有不同的差分法，以下個別介紹其方法、原理與 優缺點。 (1)二軌跡差分干涉法
此法使用現有的數值高程模型(Digital Elevation Model, DEM)模擬為干 涉圖做為地形對；另外使用兩張 SAR 影像拍攝期間發生地形變化產出的干
及 DEM 資料是否能與雷達影像精確套合(Zebker et al., 1994)。而此法最大的 限制則為研究區是否有數值高程模型資料。 (2)三軌跡差分干涉法 此法選用三張 SAR 影像，令第一張 SAR 影像作為主影像，且假設與第 二張 SAR 影像拍攝期間地表無變形發生，故使用第一、二張 SAR 影像產生 的干涉圖作為地形對；而若地表變形發生於第二張影與第三張影像拍攝期間， 使用第一、三張 SAR 影像產生的干涉圖作為變形對，此干涉圖包括了地形 效應與地形變化造成的相位值，再將兩幅干涉圖進行差分處理後，即可得到 純粹由地形變化造成的相位值。 三軌跡差分干涉法的優點在於不需使用現有 DEM 資料，因此也不會因 DEM 品質不佳導致成果不好，且兩幅影像都計算至同一參考坐標系統，不 需再進行套合的處理(Zebker et al., 1994)。但在處理過程中，此法必須先將 地形對進行全相位回復，而處理過程中容易引起錯誤，且以干涉圖作為地形 對會有許多地方因相關性不好而沒有干涉條紋(Massonnet and Feigl, 1998)， 使純粹地形效應的干涉圖求得困難。 (3)四軌跡差分干涉法 此法需使用四張 SAR 影像，令第一張 SAR 為主影像，且假設與第二張 SAR 影像(副影像)於拍攝間隔內地表無變形，故使用第一、二張 SAR 影像 產生的干涉圖作為地形對，當地表變形發生於第三張與第四張影像拍攝間隔 內，因此用地三、四張 SAR 影像產生的干涉圖作為變形對，此干涉圖包括 了地形效應與地表變形所造成的相位值，將兩幅干涉圖進行差分處理後，即 可獲得純粹由地形變化產生的相位值。 而四軌跡差分干涉法中，影響干涉處理品質好壞的因素有：利用干涉法 產生之地形對必須進行相位回復，過程中容易產生錯誤、以干涉圖作為地形 對會因部分區域相關性低而無法產生干涉條紋、兩幅干涉影像使用不同的主 影像，可能產生套合誤差與純粹的地形效應干涉圖求得困難。
3.3.2 DInDAR 誤差來源 使用 DInSAR 計算地形變化量，其相位差結果(Δψ)代表解析度內包含各種 誤差之綜合雷達視距方向地表變形量(Pathier et al, 2003)： ：綜合雷達視距方向地表變形量 ：雷達視距方向地表變形量(每一周期 2π代表位移量λ/2) ：大氣效應 ：地形誤差 ：軌道誤差 ：幾何及時間產生之不相關 ：雜訊及資料處理過程產生之不相關 DInSAR 技術誤差來源大致分為五類，以下就其特性與影響個別說明： (1)大氣效應 雷達訊號於發射和接收都需經過大氣層，而大氣對訊號的影響包含：訊 號衰減與訊號延遲，造成影像的雜訊增加、降低像元同調性並且導致相位計 算誤差。一般來說，當影像獲得時間之大氣分布均勻時，其影響量可以忽略； 而當大氣分布不均或影像取樣前後時間的大氣條件差異過大時，則須將大氣 效應影響量去除。
Zebker et al. (1997)以 SIR-C/X-SAR 的三種波段測試大氣效應中，大氣 壓力變化與濕度變化對相位計算的影響。測試成果顯示，濕度變化對相位延 遲的影響量遠大於大氣壓力變化造成的影響量，且波長越短的訊號受到傳遞 時間延遲的現象越顯著，如圖 3.3.2。此外，大氣延遲現象對二軌跡法與三 軌跡法造成的斜距變形誤差估算式如下：
三軌跡法： 以 L 波段在相對濕度 20%的條件為例，大氣延遲分別對二軌跡法與三 軌跡法造成約 10 公分與 14 公分之斜距變形誤差(圖 3.3.3)。 而 Fujiwara et al. (1998)的實驗成果顯示，部分地區干涉圖之相位差與地 表變形無關，而是因為大氣中的水蒸氣造成雷達穿透時間延遲，導致相位計 算誤差。 圖 3.3.2 大氣壓力與濕度變化對不同波段之相位計算影響(Zebker et al., 1997) 由圖 3.3.2 可知，濕度變化對相位延遲的影響量遠大於大氣壓力變化造 成的影響，且波長越短訊號受到的延遲影響越顯著(X 波長為 3.1 公分、C 波 長為 5.67 公分、L 波長為 23.6 公分)(Zebker et al., 1997)。
L 波段於相對濕度 20%時，大氣延遲對二軌跡法與三軌跡法造成約 10 公分與 14 公分的雷達斜距變形誤差(Zebker et al., 1997) (2)地形誤差 地形誤差之影響量隨著空間基線的長短而有所差異，當基線越短時期影 響量越小；反之，基線越長時地形誤差對雷達差分干涉成果影響越大。而在 二軌跡差分干涉法中，當使用的數值高程模型(DEM)資料有空間性誤差時， 也會造成相位計算誤差。 (3)軌道誤差 干涉途中若有軌道誤差，則干涉條紋會產生趨勢性與軌道水平方向之平 行條紋，如圖 3.3.4。 圖 3.3.4DInSAR 成果軌道誤差示意圖 (4)成像幾何與時間產生之不相關 影像間的同調性表示兩張影像中地形與地物的相似程度，其值介於 0
進行干涉處理，稱為不相關(decorrelation)或是無相關性(incoherence)。 可能造成影像同調性不佳的原因可分為成像幾何與時間間隔，其中成像 幾何造成同調性不佳的原因為雷達系統本身拍攝與接收訊號造成，詳細情形 請參考 3.1.2.2 節；而時間間隔所造成的不相關是因兩幅影像拍攝時間內地 形、地物有所改變，可能是時間間隔過長，或是季節週期性的植被與土壤含 水量改變，此類變化於熱帶雨林區較為明顯，而在沙漠地區則無明顯變化 (Zebker and Villasenor, 1992)。
(5)雜訊與資料處理過程產生不相關 (a)雜訊產生不相關 位於同一像元內的不同目標物對於雷達回波訊號可能為同相位 (in-phase)，強化 ccccc 了該像元的回波強度，但也可能因反相位 (out-of-phase)使回波強度相互而減弱，故雷達影像常呈現粒狀(grainy) 或木紋狀，難以分辨特徵，需要使用適當的濾波器突顯影像特徵。此外， 載具運行的穩定度也會影響系統相關性與訊號的精確度。 (b)資料處理過程產生不相關 (i)影像套合 同一地面點在兩幅雷達影像因幾何條件不同會產生成像差異， 其中因合成孔徑雷達係利用都普勒原理提升軌道方向解析度，因此 在軌道方向之成像幾何中，為使兩幅影像正確對應，須採用同一都 普勒中心值來產生像元之複數影像值。兩幅 SAR 影像進行差分干 涉處理時，必須正確計算影像對內每一像元之正確對應位置。像元 位置正確對應，共軛像元間才有明顯相關性；若套合精度不好，則 像元間之相關性會明顯降低，甚至不相關。根據實際分析顯示，當 套合精度達到 1/10 像元時， 可以產生高達 0.98 之相關値(Hanssen, 2001)。 (ii)多觀點平均
雷達影像中，藉由增加觀點數(the number of looks)以進行多觀 點平均可有效減少雷達之斑駁(speckle)現象，但進行多觀點平均時 會失去部分相位資訊。因此，藉由多觀點平均可有效改善訊號品質 (Zebker et al., 1994)， 但亦會減低空間解析度。 (iii)不確定性 干涉圖中雷達與地表間距離，僅其分數值保存於相位；而其餘 整數部份為不確定值。此不確定值須藉相位回復其整數值。最簡單 的回復方法為：沿著一條路徑計算條紋數。雖然目前已有許多自動 方法可以進行相位回復，但是無法確保最終結果都是正確的。 (iv)相位變化值符合要求 干涉的條件中，鄰近像元間之相位變化不能超過一條干涉條紋 (2π)， 若超過此臨界點，像元間就會產生不相關。以 Envisat 衛星 為例，進行偵測變形時，一個條紋代表變形 2.8 公分，若鄰近像元 的變形量於雷達視距方向變化超過 2.8 公分，則像元間產生不相關 且變形量無法形成干涉條紋。 3.4 PSInSAR 技術 DInSAR 技術應用於大範圍地表形變的測量相較於其他監測技術有其優勢， 但 DInSAR 處理過程常因影像間的拍攝間隔過長使地表特徵改變或空間性誤差 如：軌道效應，導致低同調性(decorrelation)，造成 DInSAR 技術僅適用於平原都 會區，而山區與植被較厚的地區無法形成干涉條紋。Ferretti et al. (2000)提出永久 散射體(Permanent ScatterersTM )的概念，是於多張 SAR 影像中找出具有較高訊雜 比(SNR)的共同固定點，提取出此點位於時間序列影像上的相位值進行相位解算 (unwrapping)，便可得到點位的長期地表位移量。
(Persistent Scatterers InSAR, PSInSAR)，此法不需要給予研究區域的線性平均速 度場，亦可得到非線性變形量，如：火山地區的不規則變形(圖 3.4.1、圖 3.4.2) 永久散射體的概念是找出雷達波對地表反射回波中具有穩定的散射特性(圖 3.4.3)，而根據實際情形發現，具有高相關性的點位多為建築物、橋墩、基座或 裸露的岩石，此類物體皆具有長時間不易改變之特性，所獲得的雷達回波訊雜比 較高，因此稱此類點位為永久散射體。得到這些永久散射點於時間序列中的相位 值後，根據衛星基線與地型效應誤差的相關性、與假設一特徵尺度範圍內大氣效 應相似的情況下，可消除對 DInSAR 成果影響最多的兩個因素，達到更精密的觀 測，且於一般情形，都市地區每平方公里有數十個 PS 點、而郊區也有數個點， 此資料密度已遠高於 GPS 站、水準點密度。 圖 3.4.1 加州長谷火山口地區平均地表位移量(Hooper et al., 2004) 圖 3.4.2 加州長谷火山口地區於時間序列之地表變形量(Hooper et al., 2004)
圖 3.4.3 永久散射體雷達回波示意圖(Hooper et al., 2007) (a)一般散射體於一個像素內雷達回波分佈；(b)永久散射體在一像素內雷達回波 分佈，提供穩定訊號 3.4.1 PS 選點 於 DInSAR 成果中包含了影像的強度資訊(amplitude)、影像同調性資訊 (coherence)與相位資訊(phase)，PS 點為影像中相位較穩定的點位，其中相位穩定 與雷達回波的強度與影像同調性為正相關，因此須分別從強度影像及同調性影像 中挑出 PS 候選點(PS Candidate, PSC)後，再取交集分析其相位，得到於時間序列 內向為穩定之點位(Hooper et al., 2004)。 Ferretti et al. (2001)提出從強度影像中篩選合適的 PSC，其公式為： 其中 為影像中像素的強度分布指標； 為像素於時間序列內的強度標準差； 則為像素於所有影像中的平均強度。當一像素的 值高於設定門檻，表示此 像素於時間序列內屬於高強度且穩定的訊號，故選為 PSC；而同樣的方法使用於 同調性影像分析後亦可選出一組 PSC，並將強度影像與同調性影像取出的兩組 PSC 取交集，再進行相位穩定性分析。
為地表於 LOS 方向的變形產生之相位變化； 為大氣效應造成的相位
變化； 為軌道誤差造成的相位變化； 為 DEM 誤差造成的相位變化；
為雜訊，包含了地表散射體的熱效應(thermal noise)與影像套疊誤差
(coregistration errors)等；而 W 則代表將上述分量轉換至 間的運算子(phase wrapped operator)。 Hooper et al. (2004)以 x 像素為中心，設一半徑 R，假設半徑內的 PSC 雜訊 屬隨機分布，計算平均值求 中非雜訊向的估計值，此法實際為利用一低通濾 波(Low-pass Filter)降低雜訊影響，但因雜訊隨不同區域與不同資料有所變化，因 此仍須先確定於半徑 R 的大小中雜訊成隨機分布。 Hooper et al. (2007)改善此法，對干涉圖重新網格化取樣，利用二維傅利葉 轉換(2D-FFT)至頻率域後，利用代通濾波(band-pass filter)與低通濾波決定各網格 內的主頻波段與通過區段，屬一自適性濾波(adaptive filter)，因此得具有空間相 關性的 ，再將式 3-29 減去 得： 、 、 三項屬於高空間相關，因此一 PSC 相位差值減其周圍濾 波後之相位差值時，此三項的值甚小，合併以 表示，故將式 3-30 改寫為： 亦被稱為地形殘差(topography residual)，是因兩幅影像拍攝時衛星的軌道位 置不同，扣除地形效應時未能完全移除的量，由 Hooper et al. (2007)的推導可表 示為： 為雷達波長； 為垂直基線長； 為衛星視角； 則為衛星視角差異造成之誤 差，此誤差可表示為： 為 DEM 與實際拍攝情形間的高程誤差； 為相位中心與像素於測距(range)的
水平距離， 為第 i 幅影像的衛星視線與地表之夾角；r 為主影像與地表的視線 距離。此項誤差中因 DEM 與實際高程間差異難以估計，且假設 DEM 無誤差時， 仍存在視角誤差與像素側向偏移的貢獻量，僅有部分研究利用不同的載波頻率得 以分離 和 的貢獻量(Kampes, 2005)。上述 與 的線性關係因 和 頻率特徵相似，故可沿用並代入式 3-31： 因 、 與 可分別由影像中的相位資料、衛星軌道資訊中獲取，且 、 和 無相關性，因此可利用最小二乘法(least squares inversion)估算出 ，
且因 DInSAR 屬單一主影像，故來自主影像的 與 皆為定值，根據最小二 乘法計算得殘差回推後可求出，再與估算得到的 代入式 3-31 得： 式 3-35 中僅剩下屬於雜訊之殘差，假設 且 ，代入 Hooper 於 2004 年提出各像素的像為 SNR 公式得： 即為第 x 像素的相位餘時間序列中的穩定程度；N 為干涉圖之數量。經多次迭 代計算後 收斂。與前人研究比較， 不考慮振幅強度，而考慮相位於時間中的 變異，若一區域的相位經 2D-FFT 後仍呈現失相關(decorrelation)，則建立於空間 相關性以估計誤差的 將不適用。 代表相位穩定之參數 可作為選取 PS 點的標準，而其門檻值選定則是加入 振幅參數 計算 PSC 為 PS 點之機率，同時也重新評估 的像素，此法利
識別像素是否為 PS 的 q 值則可定義為： 其中 與 呈線性相關，利用最小二乘法計算 ，求得 的 最佳解後，選擇 之像素作為 PS 點，並可迭代計算取得精度更高之 PS 點(Hooper et al., 2007)。 3.4.2 相位解算 篩選出 PS 點後，須進行相位解算(phase unwrapping)(圖 3.4.4)，但因 PS 點 位分布零散，必須建立點與點之間的關聯性，因此使用 Delaunay 三角網，其根 據任三點外接圓不可包含第四點、當有兩種以上的分割三角形時最大化最小角， 兩種條件生成唯一的 Delaunay 三角網，再進行後續相位解算。
PSInSAR 在相位解算與 DInSAR 相異處在於 PSInSAR 屬於三維資訊(二維空 間與時間)，因此 DInSAR 的解算為二維相位解算，Hooper and Zebker(2007)改良 了原有的二維相位解算法(2-D unwrapping)－最小費用流算法(Minimal Cost Flow, MCF)，加入了時間域的三維相位解算(3-D unwrapping)是使用 Chen and Zebker 所發表的 SNAPHU 函式(statistical-cost, network-flow algorithm for phase
unwrapping)，先計算 PS 點於時間序列的相位變化 ，扣除 DEM 誤差 與 主影像上的非空間相關分量 (Chen et al., 2000；2002)，使式 3-30 改為： 而式 3-38 經 3D-unwrapping 處理後，各 PS 點可表示為： 為 相位解算所得的位移量， 為未知的整數模糊度 (integer ambiguity)， 是在ㄧ解空間搜索(solution space searching)得一唯一值，
其利用鄰近像素相減可扣除。
三維相位解算後仍包含了許多的空間誤差，但因為 對時間為高相關，因
估算出空間誤差，由於先前計算中對相位的穩定性已進行篩選，所以 DEM 誤差 與雜訊誤差造成的影響甚小可忽略不計，而軌道誤差則利用精密軌道(precise orbits)扣除，最後剩下的空間誤差則應屬於大氣效應造成的影響。將空間誤差剔 除後，則可獲得 PS 點位在 LOS 方向上的真實地表變化量。 圖 3.4.4 相位解纏示意圖 經相位解纏步驟將原始相位的整數週波值計算後還原相位 3.4.3 SBAS PSInSAR 技術雖利用篩選雷達回波穩定的 PS 點大幅改進 DInSAR 的精度， 得到較佳的地表變形監測成果，但影像的軌道垂直基線長度過長導致幾何形狀不 佳可能會降低雷達回波的相關性，因此在散射體密度較低的區域不易選出 PS 點。 短基線差分干涉技術(SBAS)是由 Berardino et al. (2002)提出，限制影像間垂直基 線長度以避免影像幾何的去相關問題，並使用多主影像的概念，增加影像間的相 關性，且確保影像對無孤立之群組，進而降低兩影像間的物理效應造成的誤差(包 含垂直基線、時間間隔與都卜勒效應)，經 PSInSAR 解算後的 SLC 影相進行干涉， 軌道(azimuth)方向上過濾並剔除無部分重疊之都卜勒波譜(doppler spectrum)，側 距(range)方向減去幾何影像造成的低相關雜訊。而影像對的選取模型為：
同樣使用 PSInSAR 之選點方法(式 3-28)先選出候選點，再根據候選點的相 位資訊篩選出較佳的點位，而點位篩選模型同式 3-36，此類於短時間間隔內通過 空間的高相關濾波之像素稱為 SDFP 點位(Slowly-Decorrelation Filtered Phase Pixels)。雖然 SDFP 點與 PS 點的篩選條件有所不同，但都是針對相位穩定的點 位進行選取，因此多為重複選取，亦表示重複選取之點位具有更高的可信度，所 以後續計算中，因 SDFP 的點位篩選門檻較高，因此可將其視為精度較高的點位， 並與 PSInSAR 成果進行合併，提高成果精度。 3.5 TCPInSAR 於 PSInSAR 與 SBAS 技術中，相位解算是個關鍵的步驟，用以估計相位整 數波數，過去主要採用兩種方法，一種是原有的二維相位解算法－離散的最小費 用流算法(sparse minimal cost flow)對每一幅干涉圖進行相位解算，但因影像間相 關點無法精確選定，於部分相關點位可能有較大的相位噪聲，導致此類點位的整 數模糊值無法正確估計，進而造成相位解算誤差；另一目前常用的方法則為先行 建置一時域相關指標(temporal coherence)，經由解空間搜索後估計得形變量與 DEM 誤差，但在參數搜索的過程中可能因多組參數(即形變量與 DEM 誤差)造成 時域相關指標出現最大值，導致此估計方法得到不唯一解，即對相位模糊度 (phase ambiguity)的估計錯誤。為了改正相位解算可能失敗的因素，Zhang et al. (2011a；2011b；2012a；2012b)提出了時域相關點雷達干涉(Temporility Coherence Point InSAR, TCPInSAR)技術，利用影像對間相關點的相位差向量進行最小二乘 法解算，藉由殘差的大小區分此相關點的相位序列有無相位模糊度存在，取代了 可能產生錯誤或失敗的相位解算步驟，再估算出地表形變速率。
3.5.1 TCP 選取
TCPInSAR 的處理需延續 DInSAR 的成果，利用 DInSAR 步驟形成的強度影 像(intensity)，開一罩窗對各組主副影像對的每個像素進行估算，此像素於兩幅 影像偏移量標準差：
表示第 x 組影像對中第 i 個像素的偏移量標準差；N 為罩窗尺寸； 為過取樣 參數，於計算中可視為常數； 為影像對於此罩窗的交互相關性，公式為： 根據此方式 Zhang 估算影像對中香港機場的偏移量，如圖 3.5.1，於海面或 山區中，因無穩定的反射體，估算出的偏移量較大且方向隨機，而在地表有人工 建物或有較強反射體處，幾乎無偏移現象。對所有影像對進行影像間偏移量計算 後，對每組影像對所得成果點取交集，得到在時間序列中偏移量較小的點位，即 為 TCP 候選點。 使用此法估計偏移量標準差時，其中式 3-41 中的變數罩窗尺寸與像素相關 性的關係如圖 3.5.2，可看出像素的相關性低時，罩窗尺寸的改變會使其偏移估 計量變化大，反之在像素的相關性高時，罩窗尺寸對偏移量估計無顯著影響。 (Zhang, 2012b)
圖 3.5.2 罩窗尺寸與像素相關性改變在估計偏移量標準差的影響 (Zhang , 2012b) 過去僅使用像素相關性作為判斷其是否為強散射體的依據，但像素相關性的 計算受到了罩窗尺寸的影響，只能以操作者的經驗來設定罩窗尺寸。但在低相關 區域中，罩窗尺寸的改變造成像素相關性計算產生偏差量(bias)，Touzi et al. (1999) 提出了估計相關性偏差量的公式：
B 為偏差量； 為 Gamma 函式； 為 Generalized hypergeometric 函式；N 為
罩窗尺寸； 為像素相關性。當假設 為 0.1、0.5 與 0.8 時，可得圖 3.5.3 成果， 從圖 3.5.3 (a)可看出固定像素相關性時，罩窗尺寸的改變會產生相關性偏差量， 且相關性低時，罩窗尺寸改變所造成的相關性偏差量越顯著；反之對高相關性的 像素，改變罩窗尺寸後計算得相關性仍一致。因此於過去的強散射體篩選時，於 非高相關區域(水體)仍可能有相關性高的點位，即為罩窗尺寸導致選點錯誤。但 用於像素偏移量的估計(圖 3.5.3 (b))時，可觀察出其偏移量趨勢與相關性偏差量 十分相似，即改變罩窗尺寸時，相關性低的像素計算出的偏移量較大，而相關性 高的像素計算出的偏移量不顯著，並隨著罩窗尺寸增加時其偏移量趨於穩定，故 此特性顯示偏移量可正確的辨識出強散射體。
圖 3.5.3 罩窗尺寸對像素相關性與偏移量計算影響(Zhang, 2012b) (a)罩窗尺寸對像素相關性計算影響；(b)罩窗尺寸對像素偏移量影響 根據上述演算法對每組主副影像對計算偏移量後，可用一矩陣形式存取副影 像中每個像素 的偏移量標準差為： 表影像中任一像素之偏移量，其中 i=1,2,3,…, 、j=1,2,3,…, 。利用直方圖統 計(histogram)可得到影像的偏移量峰值 ： A 為初步門檻值，剔除偏移量過大的像素，以確保大部分的強散射體被選取。因 先前提及增加罩窗尺寸時估算的偏移量會趨於穩定，故後續進行迭代計算，增加 罩窗尺寸(由 4 4 至 64 64)，則每個候選點可得到一組偏移量序列，再計算此序 列的標準差，當標準差小於套合精度(如 0.1 pixel)時，保留為 TCP 候選點，其餘 剔除。最後使用一多項式擬合此 TCP 候選點於時序列的偏移量，若時序偏移量 中有無法與其他時間段偏移量擬合至一平滑的多項式時，亦將此像素剔除，得到 最終的 TCP，各幅影像的偏移量計算得，亦用於每幅影像上的 TCP 套合。
長度，於點位的關聯上產生部分空洞，因此 TCPInSAR 改採用區域型的 Delaunay 三角網，於圖幅中設置多個網格點，對網格點一定範圍內的 TCP 進行 Delaunay 三角網組成，可避免三角網形邊長過長、部分區域空洞的情形，且可於加強鄰近 點位的關聯性(如圖 3.5.4)。
圖 3.5.4Delaunay 三角網建立方法(Zhang et al., 2011b) (a)為全域型三角網 (b)為剔除含相位模糊的全域型三角網 (c)為區域型三角網 (d)為剔除含相位模糊的區域型三角網 3.5.3 TCP 參數估計 進行參數估計前，須根據 DInSAR 成果取出 TCP 的各種資訊，如：空間垂 直基線、時間基線、影像對組成等。假設時間序列中有 J+1 幅影像，進行干涉後， 在干涉圖 i 中可得到 TCP(l,m)在主影像(M)與副影像(S)間的 LOS 方向形變量： (l,m)為相片座標； 與 為拍攝主副影像時感測器至目標之斜距； 為在 至 時間序列內，所有 SLC 影像的數量。當變形量為線性時，則可假設 ，則其形變相位可寫為：
為雷達波長； 為對應之係數， 。 為兩影像組成的影像對矩陣， 。各 TCP 像素的觀測相位值為： W 為代表將上述分量轉換至 間的運算子； 為地形誤差造成的相位變 化； 為目標物形變產生的相位變化； 為大氣影響導致的相位變化； 為軌道誤差產生的相位變化； 為都卜勒中心頻率誤差產生的相位 變化； 為剩餘雜訊。其中 與 DEM 高程誤差的關聯為： 為影像對之間的垂直基線長度； 為主影像至目標物的斜距； 為衛 星視角；簡化後的 。兩個 TCP 間的相位差 可寫為： 因大氣影響量屬空間高相關，鄰近的 TCP 受到大氣影響量很小、經差分計算後， 都卜勒中心頻率誤差產生影響亦可降低，其影響量甚小、軌道誤差經差分計算後 對鄰近 TCP 相位差有相似的特徵，因此期望值 應為零。故對三角網 組成的每條弧(arc)之觀測方程式可列為：
為影像對中兩相鄰像素的相位差之矩陣；A 為一設計矩陣，由高程轉相位轉 換參數、時間參數組成；K 為一整數向量，代表其相位模糊(phase ambiguity)； W 為隨機參數向量。 使用區域型的 Delaunay 三角網選取弧(arc)時，因弧長受到限制，理論上短 弧中無相位模糊情形，故設 K 為零，式 3-51 中相位差矩陣可改為： 雖然採用短弧，但式 3-52 中的觀測相位( )仍可能存在整數模糊度，將於 後續利用殘差檢測器找出含有整數模糊度的弧並剔除。 3.5.3.1 先驗精度 過去的 InSAR 計算中，會將觀測點位視為等權，但此假設與實際情形不符， 於不同時間所獲的的 SAR 影像可能受到相異的大氣延遲與雜訊影響，以 SAR 影 像的隨機雜訊方差-協方差矩陣(Qnoise)可寫為： 其表示每一幅影像的平均雜訊精度，而根據誤差傳播，進行干涉計算後，每幅干 涉圖的隨機雜訊方差-協方差矩陣可寫為： D 表示為干涉圖是由某兩張 SAR 影像組成。而在計算每條弧的相位差時，同樣 根據誤差傳播，每條弧的隨機雜訊方差-協方差矩陣可寫為：
但大氣延遲與軌道誤差於弧上造成的相位誤差因短弧可假設其弧的兩端大氣延 遲與軌道誤差情況相當，因此 可視為零；而干涉圖中，所有的點應有相 同的固有雜訊(Kampes, 2006)，故可假設 與 相等，因此可將每條弧的 隨機雜訊方差-協方差寫為： 雖以影像平均雜訊作為每條弧的雜訊方差計算不夠精確，但仍考慮各幅影像品質， 故可作為先驗精度之依據，根據方差-協方差矩陣，得先驗權矩陣： 3.5.3.2 初步參數估計 根據最小二乘法計算，相位差的解為： 其中 表示為估計量； 為初步計算之殘差。而各參數的方差矩陣可寫為： 3.5.3.3 剔除相位模糊 使用最小二乘法估計參數時的基本假設為觀測中無粗差(outlier)與系統誤差 存在，於先前的計算中，假設所有弧皆無相位模糊影響，此為明顯的錯誤假設， 因此應於後續將含有相位模糊者篩選並剔除。
