利用多元衛星影像監測格陵蘭Russell冰河之變動行為與消融機制分析 - 政大學術集成

全文

(1)國立政治大學地政學系私立中國地政研究所. 碩士論文. 利用多元衛星影像監測格陵蘭 Russell 冰河之變動行為與消融機制分析 A Remote Sensing Monitoring of Greenland Russell Glacier Dynamics and Analysis of Melting Mechanism. 研究生：蔡亞倫指導教授：林士淵. 中. 華. 民. 國. 一. 零. 六. 年. 五. 月.

(2)

(3) 謝誌首先感謝我的家人，媽媽、姐姐和爸爸，不論在生活或情緒上的支援與幫助，沒有你們我不可能專心而安心地忙於研究，也才有機會好好唸書並完成研究所的學業。再次感謝你們的付出與包容。也非常感謝系上的各位老師，尤其是我的指導老師－林士淵老師。若不是在大二時您應允讓我參與研究，我今天可能不會聽過合成孔徑雷達，也不會跨入研究的殿堂。也感謝老師給予我極多的機會發揮所長，我才能參與多場國內外研討會與各界先進合作討論，甚至出國學術交流並發表外文學術期刊。以上種種均讓我得以在比別人更輕的年紀就可一窺學術界的型態與規矩，令我更快對未來研究生涯有更清晰的認知。也感謝您時常給予我幫助，不僅是研究上的方向指引，更提點我如何作為一個學術者來思考各種問題。並感謝詹進發老師，您讓我自大一起就擔任系上電腦教室管理員，而有機會吸收最基礎而重要的電腦相關知識，才能在後續的研究生涯中更有效率的處理各種資料，並也讓我了解英語能力的重要性。並感謝各界諸學者的教誨與分享，包含中研院地球所林正洪老師，耐心指導我關於火山與斷層的各種知識；中興社李璟芳博士與黃韋凱，傳授我各種關於崩塌地的基本知識；並感謝國防大學環工系林玉菁老師提供本研究使用之 SARscape 軟體和學術交流的機會。最後謝謝研究室的各位同學與學長姐們，為我的研究生生涯帶來了值得回憶回味的張張風景。 In addition, I would like to thank Prof. Dr. Jung-Rack Kim for giving me lots of tasks to train my processing and analysis skills in every aspect. Without those valuable experiences, it is impossible for me to become a learning researcher. I also appreciate those stories you shared with me, telling me the tough situation but also the bright side of the academic world. 碩士畢業，是作為學者的起點，也祝我自己未來一切順利。蔡亞倫筆民國一○六年五月二十五日臺北南港 3.

(4) 摘要近年全球暖化現象日益嚴重，格陵蘭等極區融冰所造成之海平面上升將對全球人類帶來嚴重威脅。因冰層質量之改變與冰河移動速度高度相關，故可藉由監測格陵蘭冰層（Greenland Ice Sheet，GrIS）上冰河之移動推估全球暖化對其造成之影響。衛星影像因具有連續且快速獲得大範圍地表資訊之能力，且可結合各影像處理技術獲得地表變形量，故已廣泛應用於廣域冰河之監測。然不同影像與技術均有其優勢與限制，故本研究將使用合成孔徑雷達（Synthetic Aperture Radar，SAR）與光學影像，並結合合成孔徑雷達差分干涉（Differential Interferometric SAR，D-InSAR）、多重合成孔徑雷達干涉（Multi-aperture Interferometric SAR，MAI）與偏移偵測法（Pixel-offset，PO）技術獲得冰河表面於不同方向之位移向量，再整合各向量透過三維變動量解構法（3D decomposition）求解表面於三維方向之變形量。據此執行數值冰層動力模型（Numerical Ice Sheet Model，ISM），並結合模擬之冰底基岩渠道網絡、數化之冰面冰隙與冰面湖及氣象觀測資料後，參佐冰河變動理論，進一步了解格陵蘭 Russell 冰河之變動行為與機制。關鍵字：格陵蘭冰層、合成孔徑雷達、偏移偵測法、數值冰層動力模型. 4.

(5) Abstract Global warming has been a worldwide issue and significantly increasing icecap melting rate over polar area. Consequently the sea level rises continuously and poses a fundamental threat to whole human beings. Since the mass loss of Greenland ice sheet （GrIS） is highly correlated to the velocity of glacier movement, this study aims to monitor the impact of global warming by tracking glacier terminus displacement over GrIS using remote sensing techniques. As there are multiple spaceborne images of various characteristics and also multiple techniques with different functions, we proposed a monitoring strategy using Synthetic Aperture Radar （SAR） and optical images, with Differential Interferometric SAR （D-InSAR）, Multi-aperture Interferometric SAR （MAI） and Pixel-offset （PO） techniques to estimate glacier movement vectors. The vectors were then merged using 3D decomposition method to derive 3D deformation. Based on the resultant 3D deformation, the Numerical Ice Sheet Model （ISM） is conducted and then integrates with modeled subglacial drainage channel network and glaciological theories, the melting dynamics and mechanism of Russell glacier can be further understood. Keywords: Greenland Ice Sheet （GrIS）, Synthetic Aperture Radar （SAR）, Pixel Offset （PO）, Numerical Ice Sheet Model （ISM）. 5.

(6) 目錄摘要 ...................................................................................................................... 4 ABSTRACT ......................................................................................................... 5 目錄 ...................................................................................................................... 6 圖目錄 .................................................................................................................. 8 表目錄 ................................................................................................................ 14 第一章. 緒論 .................................................................................................... 16 第一節研究背景與動機 ...................................................................... 16 第二節研究目的 .................................................................................. 18. 第二章. 理論基礎與文獻回顧 ........................................................................ 19 第一節格陵蘭冰層近年位移狀態 ...................................................... 19 第二節監測冰河變形之遙測技術 ...................................................... 22 第三節合成孔徑雷達與光學影像偵測變形技術 .............................. 24 第四節三維變動量解構法 .................................................................. 34. 第三章. 冰河變動監測方法 ............................................................................ 37 第一節研究區域 .................................................................................. 37 第二節研究流程 .................................................................................. 42 第三節研究工具 .................................................................................. 45. 第四章. 冰河變動行為與成果檢核 ................................................................ 52 第一節資料選取 .................................................................................. 52 第二節使用影像 .................................................................................. 57 第三節前處理與參數測試 .................................................................. 66 第四節各方向位移量成果 .................................................................. 97 第五節三維位移量成果 .................................................................... 136 6.

(7) 第六節成果檢核 ................................................................................ 148 第五章. 冰河變動機制分析與討論 .............................................................. 155 第一節冰河結構物與變動機制理論 ................................................ 155 第二節數值冰層動力模型介紹 ........................................................ 160 第三節 RUSSELL 區域冰層構造與速度分析成果 ............................. 164 第四節數值模型模擬與分析成果 .................................................... 171 第五節冰層表面與底層之結構物分析成果 .................................... 179 第六節結構物與冰面速度與 ISM 成果之比較 ............................... 184. 第六章. 結論與建議 ...................................................................................... 189 第一節結論 ........................................................................................ 189 第二節建議 ........................................................................................ 193. 參考文獻 .......................................................................................................... 194. 7.

(8) 圖目錄圖 2 - 1 GrIS 於 2008-2009 之地表變化速度場（Rignot and Mouginot, 2012） ............................................................................................................................ 21 圖 2 - 2 冰雪覆蓋區尚無法有效偵測控制點（Gleitsmann and Kappas, 2006） ............................................................................................................................ 23 圖 2 - 3 D-InSAR 步驟示意圖（本研究繪製） ............................................. 25 圖 2 - 4 PS-InSAR 步驟示意圖（本研究繪製） ............................................ 26 圖 2 - 5 PS 點與 DS 點於真實世界之示意圖（Ferretti et al., 2009）........... 27 圖 2 - 6 以 C 波段 ERS-1 偵測西南極洲 Dotson 冰層之平均速度場（McMillan et al., 2012） ................................................................................. 28 圖 2 - 7 以 C 波段之 ERS-2 偵測美國加州 Hector Mine 地震前後地表變動（Jung et al., 2014） .......................................................................................... 28 圖 2 - 8 MAI 原理示意圖（Hu et al., 2012） ................................................. 29 圖 2 - 9 MAI 技術流程處理圖（Hu et al., 2012） ......................................... 30 圖 2 - 10 以 COSMO-SkyMed 雷達影像執行偏移偵測法以偵測阿根廷 Viedma 冰河移動速度（Riveros et al., 2013） ............................................... 32 圖 2 - 11 因水平方向移動造成的誤差示意圖（Samieie-Esfahany et al., 2009） ............................................................................................................................ 35 圖 2 - 12 以 C 波段 ERS-1/-2 偵測美國加州 Hector Mine 地震前後地表變動（Fialko et al., 2001） ....................................................................................... 36 圖 2 - 13 以 C 波段 ERS-1/-2 針對荷蘭油田開採之地層變動之解構（Samieie-Esfahany et al., 2009） .................................................................... 36. 8.

(9) 圖 3 - 1 臨海與陸地消融型冰河之構造特徵與消融機制（Chu, 2014a） .. 38 圖 3 - 2 Russell 冰河之光學影像（本研究繪製） ......................................... 38 圖 3 - 3 Russell 冰河之地表高程示意圖（本研究繪製） ............................. 39 圖 3 - 4 Russell 冰河於 2009 與 2010 冬季之平均移動速度場（Fitzpatrick et al., 2013） ................................................................................................................ 39 圖 3 - 5 Russell 冰河其上分佈的大小 SGL（Fitzpatrick et al., 2014） ........ 40 圖 3 - 6 Kangerlussuaq 位置示意圖 ................................................................. 40 圖 3 - 7 Kangerlussuaq 年均氣候示意圖（WeatherSpark，2016） .............. 41 圖 3 - 8 衛星成像幾何各角度與方向示意圖（本研究繪製） ..................... 43 圖 3 - 9 真實地表變形量於 LOS 方向之投影示意圖（TRE-ALTAMIRA， 2016） ................................................................................................................ 43 圖 3 - 10 研究流程示意圖 ............................................................................... 44 圖 3 - 11 COSI-Corr 匹配功能選單 ................................................................. 50 圖 4 - 1 Sentinel-2 影像涵蓋範圍 .................................................................... 58 圖 4 - 2 Sentinel-2 波段光譜示意圖 ................................................................ 59 圖 4 - 3 Landsat-8 影像涵蓋範圍 ..................................................................... 60 圖 4 - 4 ASTER 影像涵蓋範圍 ........................................................................ 62 圖 4 - 5 Sentinel-1 影像涵蓋範圍 .................................................................... 64 圖 4 - 6 本研究將使用影像之時間排序 ......................................................... 66 圖 4 - 7 D-InSAR 同調性圖 ............................................................................. 68 圖 4 - 8 基岩區域視衛星方向變動量 ............................................................. 69 9.

(10) 圖 4 - 9 冰雪覆蓋區域 LOS 方向變動量 ....................................................... 70 圖 4 - 10 6 天與 12 天之 D-InSAR 成果差異 .............................................. 72 圖 4 - 11 影像對基線組合 ............................................................................... 73 圖 4 - 12 時間序列分析影像對干涉圖 ........................................................... 74 圖 4 - 13 全相位回復後之干涉條紋 ............................................................... 75 圖 4 - 14 平均移動速度場（mm/yr） ............................................................ 76 圖 4 - 15 SNAP 像素偵測法成果示意圖 ........................................................ 77 圖 4 - 16 CIAS 像素偵測法成果示意圖.......................................................... 77 圖 4 - 17 不同形式與偏極之雷達影像差異 ................................................... 80 圖 4 - 18 不同形式與偏極之偏移偵測法成果 ............................................... 81 圖 4 - 19 Sentinel-2 與 Landsat-8 之波段光譜位置差異（Kääb et al., 2016） ............................................................................................................................ 82 圖 4 - 20 Sentinel-2 各波段之地貌顯示 .......................................................... 83 圖 4 - 21 Sentinel-2 最高空間解析度之各波段光譜位置 .............................. 84 圖 4 - 22 PCA 各元素成果 ............................................................................... 85 圖 4 - 23 以不同 PCA 元素進行偏移偵測法之成果 ..................................... 86 圖 4 - 24 光學影像正射化品質與偏移偵測法成果之關係 ........................... 88 圖 4 - 25 2015/3/18-2015/7/8 之偏移偵測變動量 ........................................... 89 圖 4 - 26 不同罩窗大小之 Landsat-8 偏移偵測成果差異 ............................. 90 圖 4 - 27 不同罩窗大小之 Sentinel-1 偏移偵測成果差異 ............................. 91 圖 4 - 28 MAI 同調性 ....................................................................................... 93 10.

(11) 圖 4 - 29 不同技術之同調性比較 ................................................................... 94 圖 4 - 30 MAI 變動量成果 ............................................................................... 96 圖 4 - 31 D-InSAR 同調性圖 ......................................................................... 100 圖 4 - 32 D-InSAR 干涉條紋圖 ..................................................................... 102 圖 4 - 33 D-InSAR 干涉條紋圖 ..................................................................... 104 圖 4 - 34 D-InSAR 變動量圖 ......................................................................... 107 圖 4 - 35 N-SBAS 累積變動量 ...................................................................... 110 圖 4 - 36 興趣點累積變動量 ......................................................................... 111 圖 4 - 37 剖線累積變動量 ............................................................................. 114 圖 4 - 38 LOS 方向之變動量 ......................................................................... 115 圖 4 - 39 方位方向之變動量 ......................................................................... 116 圖 4 - 40 升軌組偏移偵測法 LOS 方向成果 ............................................... 120 圖 4 - 41 升軌組偏移偵測法方位方向成果 .................................................. 122 圖 4 - 42 降軌組偏移偵測法 LOS 方向成果 ................................................ 124 圖 4 - 43 降軌組偏移偵測法方位方向成果 .................................................. 126 圖 4 - 44 各冰舌之季節性速度變化 ............................................................. 129 圖 4 - 45 Landsat-8 偏移偵測變動量 ............................................................. 131 圖 4 - 46 Landsat-8 偏移偵測日變動量 ......................................................... 132 圖 4 - 47 Sentinel-2 偏移偵測日變動量 ........................................................ 133 圖 4 - 48 Sentinel-2 錯誤示意圖 .................................................................... 134. 11.

(12) 圖 4 - 49 ASTER 偏移偵測日變動量 ............................................................ 135 圖 4 - 50 衛星成像幾何各角度與方向示意圖 ............................................. 136 圖 4 - 51 各所用技術之可使用影像時間示意圖 ......................................... 137 圖 4 - 52 Type2 變動量解構之輸入與產出 ................................................... 141 圖 4 - 53 Type4 變動量解構之輸入與產出 ................................................... 144 圖 4 - 54 Type6 變動量解構之輸入與產出 ................................................... 147 圖 4 - 55 本研究 N-SBAS 監測之速度與 ESA 成果之趨勢比較 ............... 149 圖 4 - 56 N-SBAS 各冰舌位置速度檢核成果 .............................................. 150 圖 4 - 57 各冰舌之流線位置 ......................................................................... 151 圖 4 - 58 三維變動量各冰舌位置速度檢核成果 ......................................... 153 圖 5 - 1 GrIS 特徵物結構示意圖（Zwally et al., 2002） ............................. 156 圖 5 - 2 不同冰體模擬模式幾何示意圖（Larour et al., 2012） .................. 161 圖 5 - 3 南極大陸模擬之冰層/基岩界面摩擦係數（Morlighem et al., 2013） .......................................................................................................................... 163 圖 5 - 4 分析冰層行為與變動機制之流程示意圖 ....................................... 165 圖 5 - 5 Russell 冰河移動與冰層結構 ........................................................... 167 圖 5 - 6 熱－黏性系統崩解示意圖（Colgan et al., 2015） ........................ 168 圖 5 - 7 Russell 冰河之冰層厚度與基岩高度 3D 示意圖 ............................ 168 圖 5 - 8 冰層表面移動速度與氣象資料之比較 ............................................ 170 圖 5 - 9 ISM 模擬之成果................................................................................ 177 圖 5 - 10 各波段影像與冰隙成果示意圖 ..................................................... 180 12.

(13) 圖 5 - 11 數化冰面冰隙之成果 ..................................................................... 181 圖 5 - 12 模擬冰底渠道之流程圖 ................................................................. 182 圖 5 - 13 模擬冰底基岩渠道與流域成果 ..................................................... 183 圖 5 - 14 表面冰隙分布與冰面流線及冰底渠道之關係 ............................. 185 圖 5 - 15 冰面流線與基岩渠道之差異 ......................................................... 187 圖 5 - 16 數值冰層動力模型模擬之速度差異與冰底渠道之比較 ............. 188. 13.

(14) 表目錄表 3 - 1 D-InSAR 軟體比較（本研究整理） ................................................. 48 表 3 - 2 時間序列分析軟體比較（本研究整理） ......................................... 49 表 4 - 1 合成孔徑雷達與光學影像比較（本研究整理） ............................. 54 表 4 - 2 使用之 Sentinel-2 影像資訊 .............................................................. 57 表 4 - 3 Sentinel-2 影像時間基線 .................................................................... 58 表 4 - 4 Sentinel-2 光譜資訊 ............................................................................ 59 表 4 - 5 使用之 Landsat-8 影像資訊 ............................................................... 60 表 4 - 6 Landsat-8 影像時間基線 ..................................................................... 61 表 4 - 7 Landsat-8 光譜資訊 ............................................................................. 61 表 4 - 8 使用之 ASTER 影像資訊 .................................................................. 62 表 4 - 9 ASTER 影像時間基線 ........................................................................ 62 表 4 - 10 ASTER 光譜資訊 .............................................................................. 63 表 4 - 11 使用之 Sentinel-1 影像資訊 ............................................................ 64 表 4 - 12 Sentinel-1 影像時間基線 .................................................................. 64 表 4 - 13 光學影像與雷達影像之特性差異 ................................................... 78 表 4 - 14 三維變動量求解可能組合 ............................................................. 137 表 4 - 15 Type2 變動量解構之影像組合 ....................................................... 139 表 4 - 16 Type4 變動量解構之影像組合 ....................................................... 142 表 4 - 17 Type6 變動量解構之影像組合 ....................................................... 145 14.

(15) 表 4 - 18 ESA CCI 冰河表面速度產品資訊（Nagler et al., 2015） ............ 148 表 5 - 1 不同冰體流動方程式之差異 ........................................................... 161 表 5 - 2 ISM 模擬輸入之資料........................................................................ 171 表 5 - 3 不同模式輸出成果之比較 ............................................................... 172 表 5 - 4 使用影像之詮釋資料 ....................................................................... 179. 15.

(16) 第一章. 緒論. 本章共有三節，分別為「研究背景與動機」、「研究目的」與「文章架構」，第一節回顧研究背景說明動機並提出過往研究缺口；第二節針對該缺口提出本研究預計完成之目標；第三節則概略彰示本篇論文之章節架構。. 第一節研究背景與動機近年來在全球暖化（Global warming）現象不斷持續的情形下，全球海平面上升為伴隨而來的問題之一。全球約有六億人口居住在海平面高 10 公尺以下的海岸地區（McGranahan et al., 2007），故全球海平面上升將會對人類生活造成嚴重的影響。針對此議題，聯合國跨政府氣候變遷小組（Intergovernmental Panel on Climate Change, IPCC）以及許多科學研究（Anthoff et al., 2010; Lichter et al., 2010; Hallegatte et al., 2013; Nicholls et al., 2014）均已明確指出全球暖化與海平面上升之間的關聯性，並透過電腦模擬氣溫升高造成之海平面上升將會對於社會經濟等面向造成巨大衝擊。而導致海平面上升之主因，源自於北極極區與格陵蘭島（Greenland）冰床以及高山冰川的融解與移動（Jevrejeva et al., 2014）。格陵蘭島位於北極與北極海之間，面積約一百七十萬平方公里，其上蘊含之冰量體積約兩百九十萬立方公里，若透過約百支的冰河全數流入並融化於海中，將造成海平面升高約七公尺（Weidick et al., 1995）。然近年來，格陵蘭冰層（Greenland ice sheet，GrIS）上的冰量已快速減少（Rignot et al., 2008），其原因來自許多連接海洋的冰川於近十幾年來快速向海移動崩解與冰體表面消融所致（Joughin et al., 2004; Rignot and Kanagaratnam, 2006; Csatho et al., 2008; van den Broeke et al., 2009）。因冰層質量與冰河移動速度極為相關（Dunse et al., 2015），冰層移動速度的增快或減緩均會影響冰層質量的變化（Strozzi et al., 2008）。因此為了解格陵蘭冰層消融情況，對各主要冰川的移動速度進行監測確屬必要，. 16.

(17) 然傳統測量方式常受限於可及性低、精度不足、觀測範圍小、施測成本高、時間解析度低等限制，並無法有效率獲得連續高精度成果。為克服以上缺點，過往許多研究使用衛星影像，如合成孔徑雷達（Synthetic Aperture Radar， SAR）與光學影像進行監測，以達到長時期、大區域、高精度之目標。但利用衛星影像估算地表幾何變形量之各項技術，亦有其誤差因素與自然限制，更因冰河區域之快速移動、地貌均質度高且位於高緯度區域等狀態，使得各技術應用於冰河監測更具挑戰。如合成孔徑雷達差分干涉（Differential Interferometric SAR, D-InSAR）雖為普遍的地表變形偵測技術，然其易受到大氣延遲誤差（Hanssen, 2001; Ding et al., 2008）、高緯度之電離層干擾、（Massonnet et al., 1993; Fujiwara et al., 2000）等限制；光學影像之雲霧遮蔽區無法獲得地表資訊；均質性高之地貌令影像對位與匹配準確率大幅降低等。故整合雷達與光學影像，並以各自技術獲得之不同方向冰河變形量確實有其重要性。而整合後之冰河表面三維方向變動量，可進一步用於了解冰河於空間與時間之變動特性，基於此，本研究將使用多元雷達與光學影像，使用各項干涉與時間序列分析技術，求解各方向變形後整合為三維方向變形，進而求解冰河表面位移之狀態。同時亦評估應用各項遙測技術於監測冰河動態行為之可行性，最後提出一理想之遙測監測策略。針對造成冰河變動之機制，近年研究認為除針對表面位移進行監測外，亦需考慮冰河表面特徵物與冰層底部結構之關聯，才能完全掌握形成變動之機制。據此本研究將回顧相關冰河變動理論後，由光學影像繪製冰河表面特徵物並以基岩高程結合水文模擬估計冰層底部渠道網絡，再與動力模型反演之摩擦力資訊比較，進而了解該區域冰河之動力行為與消融機制。. 17.

(18) 第二節研究目的根據上述之研究缺口，本研究之主要目的即在於如何整合使用衛載合成孔徑雷達與光學影像，並透過各自之影像處理技術求得格陵蘭冰河表面之各方向變動量，進而求解三維方向變動量，並進行數值冰層動力模型模擬後，與觀察之冰河表面與底層結構比較以了解冰河移動之機制。欲達成之階段性目標列點如下：一、. 以多元衛載雷達與光學影像結合各類處理技術，進而計算冰河表面之位移量. 二、. 歸納各類影像與處理技術，針對快速移動之冰河區域偵測變動量之條件. 三、. 整合不同方向之位移量以求解三維方向變形量，並比較不同組合與求解條件對最終成果品質影響. 四、. 利用觀察所得之冰河表面移動速度，分析格陵蘭 Russell 冰河之變動動力行為與機制. 18.

(19) 第二章. 理論基礎與文獻回顧. 本研究為整合衛載合成孔徑雷達與光學影像偵測地表位移量之相關技術，包含傳統干涉技術與時間序列分析等，並結合近年之偏移偵測法、多重合成孔徑干涉技術等方式，以監測格陵蘭冰河於不同方向之變動量。故先回顧格陵蘭冰層近年位移狀態與其他遙測技術之限制後，介紹合成孔徑雷達相關技術並說明近年新技術之理論。. 第一節格陵蘭冰層近年位移狀態格陵蘭島位於北極與北極海之間，面積約一百七十萬平方公里，其上蘊含之冰量體積約兩百九十萬立方公里，若透過約百支之冰河全數流入並融化於海中，將造成全球平均海水面（Global mean sea level，GMSL）升高約七公尺（Weidick et al., 1995）並導致海洋環流異常（Fichefet et al., 2003; Marsh et al., 2010; Straneo et al., 2010）。然近年來，格陵蘭冰層（GrIS， Greenland ice sheet）上之冰量已快速融化而減少（Rignot et al., 2008; Bhattacharya et al., 2009）而向海之融水增加（Shepherd and Wingham, 2007; Bamber and Riva, 2010），其原因主要來自許多連接海洋的冰川於近十幾年快速向海移動崩解造成（Luckman et al., 2006; Rignot and Kanagaratnam, 2006; Howat et al., 2007; Csatho et al., 2008），然部分冰川亦出現河道縮展之情形（Pritchard et al., 2009），故知此並非唯一原因，另約五成的冰量減少係源自冰體表面消融所致（van den Broeke et al., 2009）。而各冰川移動之速度自公分至公尺等級均有，因此為了瞭解格陵蘭地區冰量減少情況，對各主要冰川之監測確屬必要。藉由對於冰川變遷監測，可了解冰量在空間與時間上變化趨勢與變動情形，並推估全球暖化對冰川造成之影響。冰層變化之時間尺度原應為百年或千年（Alley and Whillans, 1984），然由於全球暖化作用，格陵蘭島之冰河融化與移動速度均有顯著提升（Shepherd et al., 2012; Khan et al., 2014）。其作用機制為全球暖化導致空氣溫度與海水溫度增加，進而使與其接觸之冰層表面逐漸消融，且降雪增加 19.

(20) 速度不及消融速度（Box et al., 2006; Ettema et al., 2009），故於兩者作用下使得冰層質量失衡且冰層厚度緩慢改變（Holland et al., 2008）。經觀察 GrIS 發現，通常在七月至九月之夏季有較明顯移動速度（Mohr et al., 1998; Joughin et al., 1999; Luckman et al., 2006; Joughin et al., 2008; Moon et al., 2014）。監測冰層移動速度對於監測冰層狀態而言十分重要，因冰層質量與冰河移動速度極為相關，冰層移動速度之增快或減緩均會影響冰層質量變化（Strozzi et al., 2008; Khan et al., 2014）。近十年來 GrIS 消融對於海平線上升之貢獻已增加兩倍以上，自 1996 年時 0.23 mm/yr 至 2005 年時增加為 0.57 mm/yr （Rignot and Kanagaratnam, 2006），而許多較大型的冰河移動速度亦顯著提高（Howat et al., 2005; Rignot and Kanagaratnam, 2006）。且許多研究亦顯示格陵蘭島區域之冰雪累積量速度不足以補充冰河消融造成之質量損失（Velicogna and Wahr, 2006; Howat et al., 2011; Moon et al., 2012）。關於 GrIS 冰量流失情況造成格陵蘭島重力改變之情況，許多研究均透過衛星重力系統（Luthcke et al., 2006; Velicogna and Wahr, 2006; Khan et al., 2010; Schrama et al., 2011）進行觀測；亦有人以衛載及空載測高雷達（Krabill et al., 2004）、全球定位系統（Zwally et al., 2002）、數值地形模型變遷結合光學影像（Stearns and Hamilton, 2007）與雷達干涉技術進行觀測（Rignot et al., 2008），眾研究均指出近年來冰量消融之速度高過以往。而甚至透過 GRACE 衛星系統亦指出 GrIS 質量減少之速度仍在加速中（Fichefet et al., 2003; Velicogna, 2009），且經 GPS 觀測可知，高速消融區域由過往之南部區逐漸延伸至更北處（Rignot et al., 2008; Khan et al., 2010）。而 Rignot 等人將全格陵蘭島之冰河以地形、速度結構等因子（Weidick et al., 1995）分為六大區域，並對全島以三種不同雷達於 2008 至 2009 年期間進行變遷觀測，其結果顯示冰河變動方向與地形並不一致，係因內陸地區主要由冰川移動決定，而近海地區受到地形控制所致。而其成果亦顯示不同區域之冰河有不同移動速度與方向，其中以 Jakobshavn Isbræ 冰河有最 20.

(21) 快之移動速度（13.2km/yr）（Rignot and Mouginot, 2012），如圖 2 – 1 所示。. 圖 2 - 1 GrIS 於 2008-2009 之地表變化速度場（Rignot and Mouginot, 2012） GrIS 因溫室效應造成冰層消融速度加劇（Greenland, 1997），且因冰層不穩定導致大量冰棚崩落與分解之事件（Heinrich event），不僅造成海水面上升，亦因格陵蘭鄰近北極海區域，而可能影響全球環流與氣候變化（Bond et al., 1993），故為研究全球暖化之重要區域（Bretherton et al., 1990）。. 21.

(22) 第二節監測冰河變形之遙測技術因冰層質量消融速率與冰河移動速度極為相關（Dunse et al., 2015），冰層移動速度之增快或減緩均會影響冰層質量變化（Strozzi et al., 2008），故為了解格陵蘭冰層消融情況，對各主要冰川移動速度進行監測確屬必要，然因冰河區域之快速移動、地貌均質度高且位於高緯度區域等狀態，使得各技術應用於冰河監測更具挑戰。針對各項監測技術於本節中分析之。一、全球導航衛星系統（GNSS）與地面觀測傳統之地面觀測技術，包含傳統以全測站（Total Station）及全球導航衛星系統進行測量，均需使用大量人力與時間成本（Strozzi et al., 2008; Ke et al., 2013），然其成果之空間涵蓋僅為點狀，無法進行有效率之觀測（Barrand et al., 2009）。且因格陵蘭等極區之冰河涵蓋範圍極大且多處地形不適宜進行設站或人員進入，將嚴重限制現地觀測之可達性（Karpilo Jr et al., 2009）。二、攝影測量（Photogrammetry）攝影測量為一歷史悠久之常用並可獲得相當高精度之技術，亦被用於解算監測冰河移動速度（Kääb et al., 1997; Favey et al., 1999; Kääb and Funk, 1999）。然由各文獻之結論與成果均可看出其應用於冰河區域時有極大限制與挑戰。如陡峭之地勢將嚴重影響正射化品質、過度均質而低對比度之地貌特徵物於冰雪覆蓋區難以尋找（Kääb et al., 2002）（如圖 2 – 2）、匹配錯誤難透過演算法剔除而需手動處理（Baltsavias et al., 2001）、地面控制點難以設定或不足數量（Quincey et al., 2005）等，進而導致最終成果精度不高，如利用 ASTER 影像產製之 DEM 其垂直向精度僅達 18 公尺（Kääb et al., 2002）。. 22.

(23) 圖 2 - 2 冰雪覆蓋區尚無法有效偵測控制點（Gleitsmann and Kappas, 2006）三、雷射掃描（Laser Scanning）雖有部分研究成功以雷射掃描技術取得冰河表面資訊（Echelmeyer et al., 1996; Adalgeirsdottir et al., 1998），並受惠於該技術可取得較高解析度之成果且可偵測均質地貌資訊而不需控制點（Favey et al., 1999），然其亦受到多項限制。如地面雷射掃描與實地調查相同，需要大量人力與時間投入方可取得大面積資訊（Bauer et al., 2003）；空載雷射掃描則為確保解析度而無法以過高航高飛行而增加成本，且雷射掃描雖於高反射率地表如均質性高之冰雪區有良好成果，然於反射率較差如碎屑覆蓋等地貌則易有不精確之成果。且部分成果中之錯誤可能源自於掃描定向之錯誤（Baltsavias et al., 2001），以上限制均導致雷射掃描用於冰河偵測之精度較差（Kennett and Eiken, 1997）。故回顧以上技術，可知地面觀測、傳統攝影測量與雷射掃描技術均有其技術限制而不完全適用於冰河區域之地貌資訊監測，因此許多研究均使用衛星影像，如合成孔徑雷達（Synthetic Aperture Radar，SAR）與光學影像進行監測，以達到長時期、大區域、高精度之目標，以下針對各項技術進行介紹。. 23.

(24) 第三節合成孔徑雷達與光學影像偵測變形技術一、合成孔徑雷達干涉技術合成孔徑雷達（Synthetic Aperture Radar，SAR）為一種主動式微波遙感探測系統，可由儀器自行發射雷達波並接收雷達回波，並以振幅（Amplitude）與相位（Phase）資料紀錄地表粗糙度、地物物徵性質等資訊。而所謂合成孔徑，係指虛擬天線藉由在沿軌方向移動之較短真實天線合成而得，此法可顯著提升方位方向（Azimuth）解析度（Rogers and Ingalls, 1969; Curlander and McDonough, 1991; Chan and Peng, 2003）。因其相較於被動式偵測而具有可穿透雲霧、日夜皆可施測、較短時間解析度與快速偵測大範圍等特點，目前已廣泛運用於地表變形之偵測（Anderssohn et al., 2009; Zhou et al., 2009; Cheng et al., 2012; Bonforte et al., 2013; Parrella et al., 2016）。針對 SAR 影像之各種特性，過往已建立完整理論，包括以衛星軌道升降分為升軌（Ascending）與降軌（Descending）兩形式；不同波長與頻率對雲霧氣候與地表偵測能力之差異（ Tsang et al., 1985; Iisaka, 1998; Campbell, 2002; Jensen, 2009）；垂直與水平發收之偏極（Polarization）（Lillesand and Kiefer, 1994; Novak et al., 1995; Werner et al., 2007）等，均為現今各項 SAR 影像應用之基礎核心。而利用 SAR 影像紀錄之相位資訊，結合光學干涉原理（Young, 1805），許多技術應運而生，列舉如下： 1.. 合成孔徑雷達干涉技術（Interferometric Synthetic Aperture Radar，InSAR）以同一地點不同時間拍攝之兩幅以上 SAR 影像，藉由相位差、振幅差. 產生干涉圖，以產製數值地形模型（Digital Elevation Model，DEM）（Rogers and Ingalls, 1969）。 2.. 合成孔徑雷達差分干涉（Differential Interferometric SAR, D-InSAR）利用兩張以上不同時期之 InSAR 成果干涉圖，其中一張干涉圖代表地. 24.

(25) 表高程，即地形對；另一張干涉圖包含地表變形及地形效應，即變形對。將兩干涉圖進行差分處理後即可移除地形效應，進而獲得地表變形干涉圖。（Gabriel and Goldstein, 1988），詳細步驟如圖 2 - 3。. 圖 2 - 3 D-InSAR 步驟示意圖（本研究繪製） 3.. 永久散射體合成孔徑雷達差分干涉技術（Persistent Scatter Interferometric SAR, PS-InSAR）為得到時間序列（Time-series, TS）之地表形變量，藉由偵測 SAR 影. 像中具有強而穩定回波之永久散射體（Persistent Scatter），解算其相位值，進一步移除地形誤差、大氣誤差及軌道誤差以解算真實地表形變速度（Ferretti et al., 2000），步驟如圖 2 - 4。 25.

(26) 圖 2 - 4 PS-InSAR 步驟示意圖（本研究繪製） 4.. 短基線差分干涉技術（Small Baseline Subset, SBAS）與傳統 PS 方式不同，選擇採用較短之空間與時間基線於 SAR 影像中. 配對，可有較高同調性（Coherence），以獲得良好之干涉條紋圖（Berardino et al., 2002）。 5.. 新式短基線差分干涉技術（New Small Baseline Subset, N-SBAS）與傳統 SBAS 相比，可允許不連續之像對集（Disconnected Subsets）. 與部分影像中同調（Coherent）之像元，並使用非統一之同調性門檻以篩選分布散射體（Distributed Scatter，DS），以增加最終成果點之數量與覆蓋率（López-Quiroz et al., 2009; Doin et al., 2011）。由以上技術演進可見，雷達干涉技術於時間方面，由單一時間段之處理進步為長時間序列分析；於空間方面則由僅包含於全數 SAR 影像上有強而穩定回波之永久散射體 PS，增加涵蓋分布散射體 DS，如圖 2 – 5。綜合此兩面向之改進，可知 SAR 干涉技術可高效率地取得大範圍之長時間連續變化速度，故近年已廣泛應用於各式地表變形之監測。. 26.

(27) 圖 2 - 5 PS 點與 DS 點於真實世界之示意圖（Ferretti et al., 2009）二、多重合成孔徑雷達干涉技術而於 2006 年，由 Bechor 與 Zebker 學者提出多重合成孔徑雷達干涉技術（Multi-apeture InSAR，MAI），自傳統 InSAR 干涉圖中分離出前視與後視之相位差值干涉圖，以觀測方位方向變形量，而此方法精度亦較偏移量追蹤法（Offset-tracking）為佳（Bechor and Zebker, 2006）。如 ERS 高同調性干涉圖利用偏移追蹤法可得變動精度為 12 到 15 公分，而利用 MAI 可得之精度為 6.3 公分（Bechor and Zebker, 2006; Jung et al., 2013）。而後續亦有學者提出改正其干涉圖相位雜訊、地面高程誤差、基線誤差與電離層誤差方法（Jung et al., 2009; Jung et al., 2011），至此 MAI 技術始廣泛應用於火山、冰河與地滑等偵測，如圖 2 – 6、圖 2 – 7 所示（Gourmelen et al., 2011; Hu et al., 2012; Hu et al., 2012; Jebur et al., 2013; Jung et al., 2014）。. 27.

(28) 圖 2 - 6 以 C 波段 ERS-1 偵測西南極洲 Dotson 冰層之平均速度場（a） MAI 技術取得之沿軌方向變動量（b）傳統 InSAR 技術堆疊之橫軌方向變動量（c）結合上兩圖之平均地表變形速度（McMillan et al., 2012）. 圖 2 - 7 以 C 波段之 ERS-2 偵測美國加州 Hector Mine 地震前後地表變動（a）傳統 InSAR 干涉圖（b） MAI 技術之干涉圖（Jung et al., 2014）傳統合成孔徑雷達虛擬天線長，係由在沿軌方向移動之較短真實天線合成而得，此法可顯著提升方位方向解析度（Liao and Lin, 2003）。而 MAI 技術即在零都卜勒位移之前後視 LOS 方向收集影像，原理如圖 2 – 8，或亦可藉由方位共同波段濾波（Azimuth common band filtering，ACBF）處理單觀點複數影像（Single look complex image，SLC）而得（Wegmuller et al., 2006; Wegmüller et al., 2009）。. 28.

(29) 圖 2 - 8 MAI 原理示意圖（Hu et al., 2012）共同波段濾波原本用於補償距離方向（Range）因地形高差造成之位移（Gatelli et al., 1994），而於 MAI 技術中用於減少兩影像在方位方向上都卜勒中心（Doppler centroid frequency）差異量，即進行相位波段對位方法。而藉由 ACBF 方法，MAI 技術會先將原始 SLC 影像利用都卜勒相位中心分為兩均分之前視與後視影像，然兩者之波段寬均為原始影像的一半，此舉亦將方位方向解析度降低量減少。而若對於干涉圖影像對，利用 ACBF 法可得到四張 SLC，即前視主影像、後視主影像、前視從影像與後視從影像。分別利用前視影像對與後視影像對各自進行傳統 D-InSAR 處理，得到兩張干涉圖後，方位方向變動量 Φ𝑀𝐴𝐼 即為此兩干涉圖之相位差值，即為 MAI 干涉圖。 Φ𝑀𝐴𝐼 = Φ𝑓 − Φ𝑏 = −. 4𝜋 𝑙. ∙ 𝑛 ∙ 𝑑𝑎𝑧. （2-1）. 而方位方向變動量𝑑𝑎𝑧 可藉由Φ𝑀𝐴𝐼 計算而得（Bechor and Zebker, 2006）： 𝑑𝑎𝑧 = Φ𝑀𝐴𝐼 ∙. 𝑙. ∙. 1. 4𝜋 𝑛. （2-2）. 式中之 l 為有效天線長，n 為分解之合成孔徑雷達波束天線寬（Beam width），通常為 0.5。以 ALOS PALSAR 與 ENVISAT ASAR 為例，其天線 29.

(30) 長度分別為 8.9 與 10 公尺，故 MAI 干涉圖上一組波紋2𝜋即表示方位方向上有 8.9 與 10 公尺的變動。而其中Φ𝑀𝐴𝐼 是由前後視共軛相乘產生，故不需要額外對位步驟，詳細之處理流程如圖 2 - 9。. 圖 2 - 9 MAI 技術流程處理圖（Hu et al., 2012）三、偏移偵測法雖然 D-InSAR 技術為現今普遍之地表變形偵測技術，然其具有空間與時間不相關（Zebker and Villasenor, 1992）、大氣延遲誤差（Hanssen, 2001; Ding et al., 2008），且只可觀測視衛星方向（Line of sight，LOS）即斜距方向（Slant range）一維變動量（Massonnet et al., 1993; Fujiwara et al., 2000; Wright et al., 2004）等限制。隨技術演進，前兩者限制已可大大解決，然最後一點則為該技術之原始限制（Ferretti et al., 2001; Berardino et al., 2002）。且其中 LOS 方向與衛星飛行方向相關，若衛星為北－南方向飛行，則 D-InSAR 可利用升軌與降軌兩方向取得兩 LOS 方向，然此兩方向測得之變動量僅對上－下與東－西方向較為敏感（Rocca, 2003）。因此，以現今 SAR 設計，D-InSAR 技術並無法取得北－南方向，即沿軌方向的變動量。為克服此限制，Michel 等人（1999）提出次像元相關法（Subpixel correlation methodology，SPC），又稱偏移量追蹤法（Offset-Tracking method）， 30.

(31) 或偏移偵測法（Pixel-offset，PO）。其以對位方式（Coregistration）偵測兩 SAR 影像振幅波段偏移量，進而估計方位方向（Azimuth）（沿軌方向）地表變形量。雖然 SPC 法精度較傳統 D-InSAR 技術為低，然其具有多項優勢。除因其相關圖（Correlogram）對於地表相位差較為敏感，故可克服對於高地表變形量限制；且亦可透過方法估計並移除電離層誤差（Raucoules and De Michele, 2010）。且因方位方向與 LOS 方向垂直，故此方法可補足 D-InSAR 技術不足，而被廣泛應用於地震偵測與冰河移動等（Fialko et al., 2001; Luckman et al., 2007; Yun et al., 2007; Hu et al., 2008; Giles et al., 2009; Grandin et al., 2009; Schellenberger et al., 2015; Ou et al., 2016）。像素偵測法（Pixel-tracking algorithm）即偵測兩影像中對應物體位移差（Offset），可針對冰河整體移動速度進行偵測。然此方法精度與 SAR 影像之空間解析度及對位精度極相關，軌道誤差、長空間基線（超過 300 公尺）導致之地形誤差、拍攝角誤差或對位失準都會形成對同一地區影像上之位移與扭曲（Werner et al., 2005）。為了得到真實冰河移動速度，如何於對位過程中，除去其他相關誤差即為重要課題。而若可正確移除誤差，針對冰河移動速度之偵測精度通常可達 1/10 甚至 1/20 像素大小（Fialko et al., 2001; Strozzi et al., 2002）。而對於偵測冰河表面移動速度，亦可使用回波強度（Intensity）、相關性（Coherence）、振幅（Amplitude tracking）追蹤法進行，此些方法均可得到 LOS 方向與方位方向之地表變動量，進而更完整地監測冰河移動，如圖 2 - 10（Gray et al., 1998; Rott et al., 1998; Gray et al., 2000; Gray et al., 2001; Strozzi et al., 2002; Riveros et al., 2013）。. 31.

(32) 圖 2 - 10 以 COSMO-SkyMed 雷達影像執行偏移偵測法以偵測阿根廷 Viedma 冰河移動速度（Riveros et al., 2013）回波強度偵測法，又稱互相關優化處理（Cross-correlation optimization procedure）（Rott et al., 1998; Gray et al., 2000; Gray et al., 2001; Strozzi et al., 2002），其利用標準化交叉相關（Normalized Cross Correlation, NCC）法偵測兩影像回波強度中幾近相似的物徵。而為了提高估計之精度，將利用過採樣技術（Oversampling）與二維迴歸模式求出相關性最大值，並以信噪比（Signal-to-noise ration，SNR）作為偏移量精度指標。其最大優點在於其可在同調性較差條件下，如時間基線較大之山區等條件下進行偵測。相關性偵測法，又稱可視條紋演算法（Fringe visibility algorithm）或相關性優化處理（Coherence optimization procedure）（Derauw, 1999），利用一系列小型干涉圖計算相關性並取得偏移量。而亦利用過取樣技術與二維迴歸模視求出相關性最大值，並以最大相關度與平均相關度差異量，作為精. 32.

(33) 度評估標準。然此些方法估計出之像素偏移量，包含因衛星不同位置與兩影像時間中地面變動所造成之變形量，其中軌道偏移與基線相關，而應以地面參考點計算並移除之（Strozzi et al., 2002）。且此些方法所得之 LOS 方向與方位方向變動量僅為空間中三維方向之兩向量，故可引入數值地形模型（Digital elevation model，DEM）並於冰河流向平行於表面之假設下求出移動方式（Joughin et al., 1998; Mohr et al., 1998）。若無 DEM 則亦可以升軌與降軌 SAR 影像進行，即先利用次像元相關比較兩 SLC 影像，產生一相關圖（Correlogram）後，即可取得 LOS 與方位方向之變形量利用矩陣運算可獲得空間三維方向移動向量（De Michele et al., 2010）。為了解並評估上述各方法於冰河監測之可行性，本節所整理之差分干涉、多重差分干涉、新式最小基線時間序列分析與偏移偵測法技術均會使用。並透過下述介紹之三維變動量解構法求解三維方向變形，以下介紹之。. 33.

(34) 第四節三維變動量解構法由以上合成孔徑雷達與光學影像各技術獲得之冰河變動量方向略有差異，如 InSAR 技術僅可求得 LOS 方向變形；MAI 與 PO 可獲得方位方向變動，然其中雷達影像之 LOS 與沿軌方向均會受到衛星飛行軌道航向方位角和地面入射角之不同而因影像而異，使得不同影像對之成果無法直接比較，且 LOS 方向僅為真實變形方向於衛星－地表連線之投影，故不便判釋地表真實變形方向。故為取得三維方向之變動量，Wright 提出可由至少三個不同 LOS 幾何方向之 InSAR 觀測量，分解出北、東、上三方向之變動量元素以克服（Wright et al., 2004），然此法因受限 SAR 影像拍攝方法而不可行。而其他學者亦提出斑點追蹤法（Speckle-tracking method）取得方位方向之變形（Derauw, 1999; Fialko et al., 2005），然此法僅適用於大規模變形或高同調性條件下。而 Hassen 等人於 2009 年提出三維變動量解構法（3D decomposition），建立在水平變動量與傾斜度（即垂直方向變形量之一階導數）成比例之假設下（Kratzsch, 1986），可僅以一幾何方向之 InSAR 成果中取得地表於三方向之變動量，進而克服雷達變動量偵測技術之限制（Samieie-Esfahany et al., 2009）。其中文章中提到，InSAR 技術所偵測之 LOS 方向地表變形量實際上是北、東與上三方向向量投影至 LOS 方向的向量總合，故雖然 InSAR 主要對於上下方向之變動較為明顯，亦不可忽略其他兩方向之向量隱含於其內，否則會因水平變動之方向不均，造成上下方向變動量局部高估或低估，如圖 2 – 11 所示。. 34.

(35) 圖 2 - 11 因水平方向移動造成的誤差示意圖（A）垂直方向變形量（B）平行方向變形量，箭頭表示變形方向（C）忽略水平方向變形對垂直方向變形之誤差（D）衛星飛行方向與方位方向（Samieie-Esfahany et al., 2009）而實際之計算流程，則是先將 InSAR 取得 LOS 方向變動量轉換為垂直變動量場，並利用克利金內差法（Kriging interpolation）將稀疏點資料建立連續變動場（Continuous subsidence field）後，以坡度核（Gradient kernel）與該網格進行二維摺積計算各位置一階空間導數值，再以此求出各位置之傾斜度與水平變動量場，並與初始垂直方向變動量場合併計算。此步驟會重複迭代執行，直至垂直變動量場有收斂解（Samieie-Esfahany et al., 2009）。除此方法外，亦有許多學者提出以不同數目或來源資料之結合計算出三維變動量值，如結合兩 D-InSAR 干涉圖與方位方向偏移量（Fialko et al., 2001）、不同攝影角度與升降軌（Wright et al., 2004; Jian-bao et al., 2006）、InSAR 與 GPS （Gudmundsson et al., 2002）、振幅偵測法（Michel Rémi et al., 1999; Hu et al., 2010）等，如圖 2 – 12、圖 2 - 13。. 35.

(36) 而本研究為檢核各方向變動量之精度與分析冰河變動方向，亦將利用多元影像配合多元技術所估算之不同方向偏移量，進而計算該區域之三維方向變動量。. 圖 2 - 12 以 C 波段 ERS-1/-2 偵測美國加州 Hector Mine 地震前後地表變動（a） 99/11/21-95/11/12 升軌干涉圖（b） 99/10/20-99/09/15 降軌干涉圖（c）降軌之方位方向變動量（箭頭表示衛星橫軌方向）（d）-（f）向上、東、北之三維地表變形分量（Fialko et al., 2001）. 圖 2 - 13 以 C 波段 ERS-1/-2 針對荷蘭油田開採之地層變動之解構（A）向上變動量（B）向北變動量（C）向東變動量（藍色箭頭代表水平之變動方向）（Samieie-Esfahany et al., 2009） 36.

(37) 第三章. 冰河變動監測方法. 本章研究方法可分為三節，首先介紹研究區域之特性；第二節介紹本研究之流程設計，並提出流程圖；第三節則說明本研究預計使用之軟體。. 第一節研究區域 Russell 冰河位於格陵蘭西南部（中心經緯度 67°N，49°W），其與一般消融於海中之冰川（Marine-terminating）不同，屬於陸地消融型之冰川（Land-terminating），故其冰量減損係因季節性表面融水之侵蝕作用所致，而非裂解為冰山（Zwally et al., 2002; Joughin et al., 2008; Van de Wal et al., 2008; Bartholomew et al., 2010）。該區域冰河大致包含 Isunnguata Sermia、 Russell、 Leverett、Ø rkendalen 與 Isorlersuup 冰河，如圖 3 – 2。而其表面上之冰面湖（Supraglacial lakes，SGL）與複雜河道景觀亦為此冰川之特點，因其廣泛分佈之融冰湖，使此處成為研究 SGL 水文動力極佳之實驗區。（Shepherd et al., 2009; Selmes et al., 2011; Andrews et al., 2014）。NASA OIB （Operation IceBridge）於 2010 與 2011 年間，利用高密度雷達進行偵測，足以彰顯此處研究價值。且又因其二十公里以西即為格陵蘭島最主要之機場樞紐，坎克爾路斯斯雅克（Kangerlussuaq）民航機用機場，故更說明監測本冰河之重要性。而此區域之冰河構造與速度等相關資訊，可見圖 3 – 1 至圖 3 – 5。. 37.

(38) 圖 3 - 1 臨海與陸地消融型冰河之構造特徵與消融機制（Chu, 2014a）. 圖 3 - 2 Russell 冰河之光學影像（本研究繪製）. 38.

(39) 圖 3 - 3 Russell 冰河之地表高程示意圖（本研究繪製）. 圖 3 - 4 Russell 冰河於 2009 與 2010 冬季之平均移動速度場（底圖為 Landsat TSX 於 2008 年七月之影像）（Fitzpatrick et al., 2013）. 39.

(40) 圖 3 - 5 Russell 冰河其上分佈的大小 SGL （Fitzpatrick et al., 2014）由文獻可知，因氣候條件改變造成之峽灣融水（Fjord water）將造成冰河移動速度與消融速率增快（Mair et al., 2001; Motyka et al., 2003; Scambos et al., 2004 ），且氣候條件之變動造成之地貌狀態改變亦會影響執行 D-InSAR 與 PO 時之成果好壞（Strozzi et al., 1999; Strozzi, Luckman, Murray, Wegmuller, et al., 2002），故知了解冰河區域之氣候狀態亦為分析時重要項目。針對 Russell 冰河之氣候條件，因該處並無測站，故採用西邊 25 公里之 Kangerlussuaq（丹麥語 Søndre Strømfjord）區域，即圖 3 – 6 紅點處，其為長約 175 公里之峽灣（Fjord），亦為格陵蘭最主要之航空樞紐與最大之商用機場，周遭亦為格陵蘭最重要之生態棲地，包含麝牛、海東青等（Eisner et al., 1995）。該區域具有高密度之自動氣象測站（Automatic Weather Station， AWS），提供豐富之氣象歷史資料（As et al., 2012）。. 圖 3 - 6 Kangerlussuaq 位置示意圖 40.

(41) 針對該區域之長期氣候型態，根據 1974 至 2012 之觀測資料顯示，如圖 3 – 7 所示，最熱之季節約為六至八月（最熱 16°C），最冷則為一至二月（最冷-25°C），同時各自亦為降雨與降雪季節，終年有約六成以上之雲霧覆蓋機率（WeatherSpark，2016）。（a）. （b）. （c）. 圖 3 - 7 Kangerlussuaq 年均氣候示意圖（a）最高低溫度（b）降雨降雪機率（c）雲霧覆蓋機率（WeatherSpark，2016） 41.

(42) 第二節研究流程為完整監測格陵蘭西南側之 Russell 冰河表面三維變動情形，並整合各自獲得位移向量之成果，本研究先以不同影像處理技術獲得不同方向之冰河表面位移量，簡要分述如下：一、合成孔徑雷達影像 A. 最小基線法技術取得一維 LOS 方向之平均速度場 B. 差分干涉技術取得一維 LOS 方向變形量 C. 多重差分干涉技術取得一維沿軌（Azimuth）方向變動量 D. 偏移偵測法取得沿軌與 LOS 方向變動量二、光學影像 A. 偏移偵測法取得沿軌與橫軌方向變動量由以上可知，各技術獲得之冰河變動量方向略有差異，包含 LOS、衛星沿軌與橫軌方向，各方向於衛星幾何之意義可見圖 3 - 8。然其中雷達影像之 LOS 與沿軌方向均會受到衛星飛行軌道航向方位角和地面入射角之不同而因影像而異，使得不同影像對之成果無法直接比較，且 LOS 方向僅為真實變形方向於衛星－地表連線之投影，故不便判釋地表真實變形方向，如圖 3 - 9。綜上所述，為比較不同感測器、影像、時間與技術獲得之地表變形量，需再透過三維變動量解構法以投影矩陣之方式，將各方向之變形量轉為三維方向（ENU）之向量。. 42.

(43) 圖 3 - 8 衛星成像幾何各角度與方向示意圖（本研究繪製）. 圖 3 - 9 真實地表變形量於 LOS 方向之投影示意圖（TRE-ALTAMIRA，2016）然為使得最後之三維方向變動量成果合理可靠，故需先以地形構造與位移量值等條件判斷各技術取得之變形量是否合理，經篩選後之合理變動量方可輸入三維變動量解構法，否則將造成最終成果之品質劣化。綜合以上所述，本研究之研究流程圖如圖 3 - 10。 43.

(44) 圖 3 - 10 研究流程示意圖. 44.

(45) 第三節研究工具因本研究將執行多項不同技術與影像處理，故使用軟體之選擇為一重要課題，除最終之三維變動量解構法將自行撰寫程式處理外，針對雷達與光學影像之各項技術處理可用之軟體分述如下。關於傳統雷達干涉與時間序列分析技術，目前已有眾多軟體供研究者或專業公司使用，但因其開發單位與目的之不同，造成各軟體間有極大差異存在，包含免費開源與否、支援之 SAR 影像感測器種類、干涉功能與演算法、作業系統環境等均大相逕庭，因此在研究選擇使用工具時，必須依照研究目的並考量所用影像與目標區域之特性後，綜合參考各軟體特性差異而定，回顧目前常用之軟體，分為 D-InSAR 與時間序列兩大項說明如下：一、D-InSAR 軟體 1. Doris（Delft object-oriented radar interferometric software） Doris 為荷蘭 TU Delft 大學地球觀測與空間系統研究所（Delft Institute of Earth Observation and Space Systems, DEOS）開發之免費且開源軟體，並支援使用者自行下載附屬套件以增減處理步驟（Kampes and Usai, 1999; Kampes et al., 2003）。其為最早之 SAR 影像處理軟體之一，然因近年並無針對新型感測器進行更新而無法處理之，故使用者數目大幅減少。 2. ROI_PAC（Repeat Orbit Interferometry PACkage） ROI_PAC 為美國加州理工學院（Caltech）與噴射推進實驗室（Jet Propulsion Laboratory, JPL）合作開發之免費軟體，亦提供使用者自行擴充套件（Rosen et al., 2004）。近年雖仍有使用者針對新型感測器進行相關套件開發，然因功能不齊全，故近年使用者已無過往多。. 45.

(46) 3. ISCE（InSAR Scientific Computing Environment） ISCE 與 ROI_PAC 之開發群相同，然為針對新型感測器影像而生，其內提供多項功能供使用者調整處理步驟，包含自動下載 DEM、裁切執行範圍等，為近年使用者較多之軟體（Rosen et al., 2011）。其另一優點為具有豐富資源之使用者討論區，其上並有開發團隊人員定期回覆各項問題，便於各項錯誤與疑問之解決，本研究將使用之。 4. GMTSAR GMTSAR 為美國 San Diego State University 海洋工程系所開發之軟體，亦支援最新之感測器影像，其特色為透過著名之出圖軟體 GMT （Generic Mapping Tools）語法以控制干涉流程（Sandwell et al., 2011b, 2011a）。然針對軟體之操作說明文件過於簡略，且並無任何討論區供發問，導致使用較不友善。 5. SNAP（Sentinel Application Platform） SNAP 為歐洲太空總署 ESA 針對其新發射之 Sentinel 系列衛星開發之軟體，然亦支援其他最新感測器影像，為舊版 STEP（Science Toolbox Exploitation Platform）軟體之延伸。其另二優點為除與 ISCE 相同具有使用者與開發團隊交流之線上討論區外，亦因其為近期開發之軟體而至今仍持續增加功能、移除軟體錯誤之更新版本，本研究將使用之。 6. SARscape SARscape 為由瑞士 sarmap SA 公司開發之商業軟體，其特色為整合於專業 ENVI 遙測分析軟體上，使得影像裁切、轉檔或濾波等功能均可快速執行，其高頻率更新更確保針對新感測器影像之處理能力，為現代最常用之商用軟體之一，本研究將使用之。. 46.

(47) 7. GAMMA GAMMA 為由瑞士 GAMMA Remote Sensing and Consulting AG 公司開發之商業軟體，亦針對各項感測器影像提供干涉處理支援。二、時間序列分析 D-InSAR 軟體 1. StaMPS（Stanford Method for Persistent Scatterers） StaMPS 原為美國 Stanford 大學地球物理系開發之免費軟體，其整合多項軟體以進行時間序列分析，包含 ROI_PAC、Doris 與 triangle 等，變動點則以挑選永久散射體為主（Hooper et al., 2012）。 2. GIAnT（Generic InSAR Analysis Toolbox） GIAnT 為美國加州理工學院開發之新時間序列分析開源免費軟體，其支援輸入其他軟體之 D-InSAR 成果，包含 ROI_PAC、DORIS、ISCE 與 GMTSAR，並提供相關組件進行大氣改正（Agram et al., 2013）。因其提供多種時間序列分析演算法，本研究將使用之。 3. SARscape 該軟體亦提供時間序列分析之模組，然其可調整之參數較 GIAnT 為少。針對上述軟體於開發團隊、免費開源與否、使用者使用之程式語言、作業系統支援、演算法差異等比較綜合如表 3 – 1 與表 3 - 2：. 47.

(48) 軟體名稱開發者免費開源程式語言*2. DORIS TU Delft O O C, C++. 表 3 - 1 D-InSAR 軟體比較（本研究整理） ROI_PAC ISCE GMTSAR SNAP JPL, CalTech O O C, Fortran, Perl. JPL, CalTech O O Python. SARscape. GAMMA. ESA O O JAVA. Sarmap X X N/A. Gamma RS X X N/A. O O O O O O. O O O*3 O O O. O O O O O O. O O X. O O O. O O O. NCC. NCC. NCC. N/A. SNAPHU. Region Growing MCF. Region Growing MCF. UC SanDiego O*1 O C, GMT. 影像支援類型 ERS ASAR Sentinel-1 PALSAR-1 PALSAR-2 COSMO-SkyMed. O O X O X O. O O X O X O. O O O O O O. Linux Mac OS Windows. O O X. O O X. O O X. O O O O O O. 作業系統支援類型 O O O*4. 對位演算法 N/A. Ampcor*5. Ampcor Denseap Nstage. SNAPHU. GRASS ICU SNAPHU. GRASS ICU SNAPHU. 全相位回復演算法. *1:GMTSAR 需安裝 MATLAB *2:使用者所用. SNAPHU. *3: 直至 5.2.1 版本無法讀取 48. *4:Win32-based. *5:包括 NCC.

(49) 表 3 - 2 時間序列分析軟體比較（本研究整理）名稱. 開發者. 免費開源作業系統程式語言. 其他功能. StaMPS. GIAnT. SARscape. Stanford University, University of CalTech Sarmap Iceland, TU Delft, University of Leeds O*6 O X O O X Linux Linux Windows/Linux MATLAB Python N/A 時間序列分析演算法 SBAS PS-InSAR N-SBAS PS-InSAR SBAS TimeFun SBAS MTI MInTS 結合 PyAPS 編輯影像組合以移除大氣誤差裁切執行範圍移除大氣誤差 deramp 移除編輯影像組合軌道誤差. *6: StaMPS 需安裝 MATLAB 三、MAI 軟體針對執行多重合成孔徑雷達干涉技術之軟體，雖然部分文獻指出可透過傳統 InSAR 軟體進行處理（He et al., 2014），然並無法得知可利用何軟體，而其他公開之免費軟體下載連結亦已失效（Barbot et al., 2008），故本研究將使用 SARscape。四、PO 軟體偏移偵測法之原理歷史悠久，最早為透過人為選擇共軛點比較物體變化量（Lucchitta and Ferguson, 1986），而該技術於早期即應用於冰河移動之監測（Stream, 1991; Scambos et al., 1992），故發展至今已有極多功能相似但演算法不同之軟體，整體而言可分為七大類（Heid and Kääb, 2012）：（1）於空間域執行 NCC （2）於頻率域執行 CC （3）於頻率域執行相位相關 49.

(50) （Phase Correlation）（4）於頻率域執行方向性 CC （5）於頻率域執行方向性相位相關（6）COSI-Corr 之相位相關（7）最小二乘法匹配。其中最後方式需設定初始值，故不常用於影像變形偵測；方法二至五均透過相同之傅立葉匹配技術但有不同之正規化方式，COSI-Corr 軟體雖亦為傅立葉處理方式，然其於頻率域匹配完後並不轉換回空間域（Leprince et al., 2007; Ayoub et al., 2009）。比較上述各演算法於冰河移動偵測之效果，COSI-Corr 除可於低地貌對比度、雲霧及雪覆蓋影像均有較其他演算法為好之成果外，亦於全球五個冰河測試區中獲得最高偵測變動量準確率（Heid and Kääb, 2012）。此可歸因於 COSI-Corr 匹配罩窗之權重為常態分配，故可避免邊緣效應；並以穩健性函數及迭代權重最小誤差法排除錯誤離群值，進而取得精度較高之位移量，故本研究中針對光學影像之偏移偵測法將使用 COSI-Corr 軟體。其匹配功能之選單如圖 3 – 11。. 圖 3 - 11 COSI-Corr 匹配功能選單針對雷達影像之偏移偵測法，因該種影像成像與紀錄之資訊與光學影像大相逕庭，使得影像具有不同之扭曲如前坡倒轉等成像限制（Lillesand and Kiefer, 1994），且受到區域性回波隨機建設與破壞性干涉造成之雜訊. 50.

(51) （Speckle）大幅降低影像品質（Chan and Peng, 2003），故需採取與光學影像不同之策略。回顧文獻，前人多利用自行開發（然並未公開）之軟體執行偏移偵測法（Nagler et al., 2015），或使用商用之 GAMMA 軟體進行處理（Strozzi, Luckman, Murray, Wegmuller, et al., 2002; Wang et al., 2015; Schellenberger et al., 2016）。免費軟體方面，ESA 之 SNAP 雖亦可解算並與 Sentinel 衛星有良好相容性，然直至 2016/12 發布之 5.0.0 版本其偏移偵測成果仍僅可顯示整體平面位移速度，並無法拆分為縱與橫方向，故分析上較有限制。然受限於 GAMMA 軟體所費不貲，本研究將使用 ISCE 軟體中承自 ROI_PAC 軟體差分干涉前處理之主從影像對位之 ampcor 演算法進像素偏移偵測。. 51.

(52) 第四章. 冰河變動行為與成果檢核. 依照本研究處理流程，本章第一節先介紹所欲使用之影像與數值地形模型；第二節詳列研究所用影像之詮釋資料；第三節為各方向偵測變形技術之前處理與參數測試成果；第四節展示以各項技術求解之各方向位移量；第五節則由各向位移量組合解算最終之三維變動量成果；最後一節則進行成果檢核。. 第一節資料選取針對本研究預計使用之合成孔徑雷達影像、光學影像與數值地形模型，簡要說明比較如下：一、雷達影像 . Sentinel-1 Sentinel-1 系列衛星係由同軌道但角度差 180 度之衛星群（Satellite. constellation）組合之系統，其意在於歐洲太空總署（ESA）接續 ERS 與 ENVISAT 系統，其上搭載 C 波長 SAR 設備，為哥白尼計畫（Copernicus）衛星，S1A 於 2014 年 4 月發射，S1B 於 2016 年 4 月發射，至今均服役中。二、光學影像 . Landsat-8 Landsat-8 為美國太空總署（ National Aeronautics and Space. Administration，NASA）Landsat 計畫於 2013 年 11 月發射之第八顆衛星，其上酬載之 Operational Land Imager （OLI）感測器確保各波段之品質精度與穩定性。 . Sentinel-2 Sentinel-2 為歐洲太空總署（ESA）之哥白尼計畫（Copernicus）發射. 之第二組衛星，與 Sentinel-1 相同為兩衛星組成，S2A 已於 2015 年 6 月發 52.

(53) 射，其上酬載之 Multispectral Imager （MSI）感測器提供各波段之紀錄。 . ASTER （Advanced Spaceborne Thermal Emission and Reflection Radiometer） ASTER 感測器裝載於 NASA 於 1999/12 發射之 Terra 衛星上，較同年. 代衛星（如 Landsat-7）相比有更高之解析度，然其重複拍攝範圍則較少而稀疏。其包含三種儀器針對不同波段進行拍攝，可見光與近紅外（Visible and Near Infrared，VNIR）、短波紅外（Shortwave Infrared，SWIR）與熱紅外（Thermal Infrared，TIR），其中 SWIR 感測器自 2008/4 起因儀器溫度異常而關閉至今。而本研究所用之影像與現今其他常用影像之比較可見表 4 – 1。. 53.

(54) 表 4 - 1 合成孔徑雷達與光學影像比較（本研究整理）名稱. ASAR. 類別管理機構免費下載. ESA X. 發射日期. 2002/3/1. 結束日期. 2012/4/8. 衛星群軌道類型太陽同步時間軌道涵蓋再訪週期（天）. 1 太陽同步 10:00 a.m. 全球 35 58~110 (SM) 405 (Scan) N/A 30~150 (SM) 150 (Scan). 像幅寬（公里）波段數空間解析度（公尺）輻射解析度（位元）. Sentinel-1. PALSAR-1. PALSAR-2. ASTER. Landsat-8. JAXA X. JAXA X. NASA O. Optical NASA O. 2006/1/24. 2014/5/24. 1999/12/18. 2013/2/11. 2011/5/12. N/A. N/A. N/A. 1 太陽同步 10:30 a.m. 全球 10（5）. 2 太陽同步 10:00 a.m. 全球 16. ESA O 2015/6/23 TBA N/A N/A 1 太陽同步 10:30 a.m. 56°S - 84°N 10（5）. 60. 185. 290. 14. 11. 13. 15/30/90. 15/30/100. 10/20/60. 8/12. 12. 12. SAR. N/A. ESA O 2014/4/3 2016/4/25 N/A N/A 2 太陽同步 10:00 a.m. 全球 12（6） 250 (IW) N/A 5~20 (IW) N/A. Sentinel-2. 1 1 太陽同步太陽同步 10:30 a.m. 11:59 a.m. 全球全球 46 14 40/70 (SM) 50/70 (SM) 250/350 (Scan) 350/490 (Scan) N/A N/A 7~88 (SM) 3~10 (SM) 100 (Scan) 100 (Scan) N/A. N/A. 54.

(55) 三、數值地形模型 . SRTM （Shuttle Radar Topography Mission） SRTM 全球數值高程模型即是使用 InSAR 技術產製之成果，此任務由. NASA 噴射推進實驗室（Jet Propulsion Laboratory，JPL）、德國太空中心及義大利太空中心合作，此計畫使用架置在太空梭上之雙雷達天線，產製北緯 60 度至南緯 57 度間，涵蓋全球面積約 80%（不包含格陵蘭），解析度為 30 公尺與 90 公尺之數值地形模型。 . ASTER （Advanced Spaceborne Thermal Emission and Reflection Radiometer）是利用 NASA 於 1999 年發射之 EOS 衛星（Earth Observing System）. 取得全球高程之資料，其涵蓋區域自北緯 83 度至南緯 83 度，包含約 99% 之地表，其解析度為 30 公尺與 1 度。 . GIMP （Greenland Ice Mapping Project） GIMP 為利用 ASTER 與 SPOT-5 之高程，與 ICEsat 於 2003 至 2009 間. 測得之高程對位而得，其解析度為 30 公尺，坐標系統為建立於 WGS84 上之 ESPG 3413（Howat et al., 2014）。. 55.

(56) 第二節使用影像為針對 2016 當年之 Russell 冰河移動情形進行監測，本研究將使用 2016 年可用之近代免費高解析度衛載感測器影像，以下針對研究中使用之雷達與光學影像群之詳細資訊與將使用之技術，分組說明：一、光學 Sentinel-2 影像選擇適用之 2016 之 Sentinel-2 影像時，雖因其軌道設計與目標區域較高緯度之條件故可較赤道區域之 10 天再訪週期有更密之影像拍攝，然仍需篩選目標冰河區域內幾乎無雲霧覆蓋之影像，故僅可挑選出時間涵蓋 2016/7/15 至 2016/10/29 之 14 張影像，其詮釋資料資訊如表 4 - 2。由表中可看出該些影像均已輻射並幾何校正之地理對位正射影像。其影像涵蓋與時間基線可見圖 4 - 1 與表 4 - 3。表 4 - 2 使用之 Sentinel-2 影像資訊 Orbit Number 68 Tile Number T22WEV Orbit Direction Descending Processing Level LEVEL-1C （ radiometric and geometric corrections along with ortho-rectification and spatial registration） Product Format JPEG2000 Datum WGS84 UTM Zone 22N Resolution （m） 10, 20, 60. 57.

(57) 30 km. 圖 4 - 1 Sentinel-2 影像涵蓋範圍表 4 - 3 Sentinel-2 影像時間基線影像時間時間基線（天） 2016/7/15 0 2016/7/21 6 2016/7/25 4 2016/7/28 3 2016/7/31 3 2016/8/7 7 2016/8/10 3 2016/8/14 4 2016/8/17 3 2016/8/20 3 2016/8/30 10 2016/9/3 4 2016/9/26 23 2016/10/29 33 而 Sentinel-2 衛星上之 Multi-Spectral Instrument （MSI）感測器可記錄自可見光至紅外線之 13 個波段，而各波段之空間解析度均不同，如圖 4 – 2 與表 4 - 4。. 58.

(58) 圖 4 - 2 Sentinel-2 波段光譜示意圖表 4 - 4 Sentinel-2 光譜資訊 Spatial Resolution Band Number （m） 10. 20. 60. Central Wavelength （nm）. Bandwidth （nm）. 2. 490. 65. 3. 560. 35. 4. 665. 30. 8a. 842. 115. 5. 705. 15. 6. 740. 15. 7. 783. 20. 8b. 865. 20. 11. 1610. 90. 12. 2190. 180. 1. 443. 20. 9. 945. 20. 10. 1380. 30. 59.

(59) 二、光學 Landsat-8 影像選擇適用之 2016 之 Landsat-8 影像時，因其軌道安排故目標區域雖於較高緯度然仍為赤道區域之 16 天再訪週期，且與 Sentinel-2 相同仍需篩選目標冰河區域內幾乎無雲霧覆蓋之影像，故僅可挑選出時間涵蓋 2016/2/1 至 2016/9/12 之 7 張影像，其詮釋資料資訊如表 4 - 5。由表中可看出該些影像均為最高校正品質等級之 L1T，與僅透過軌道資訊校正之 L1G 與加上 DEM 校正之 L1Gt 等級不同，為精確幾何校正之地理對位正射影像。其影像涵蓋與時間基線可見圖 4 - 3 與表 4 - 6。. WRS Path WRS Row Nadir/Off Nadir Orbit Direction Data Type Level 1. Image Quality. 表 4 - 5 使用之 Landsat-8 影像資訊 7 13 Nadir Descending L1T （Precision Ortho-Corrected Product） (most accurate level of processing, using GCPs and a DEM to provide systematic geometric and topographic accuracy) 9. 40 km 圖 4 - 3 Landsat-8 影像涵蓋範圍 60.