Sentinel-1 SAR數(shù)據(jù)在四川茂縣山體滑坡災(zāi)害監(jiān)測中的應(yīng)用研究

高思遠 ，黃詩峰 ，孫亞勇 ，王 慧

（1.中國水利水電科學研究院，北京 100038；2.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038）

1 研究背景

我國是世界上滑坡災(zāi)害最為嚴重的國家之一，2016年全國共發(fā)生地質(zhì)災(zāi)害9710起，其中滑坡7403起，占地質(zhì)災(zāi)害總數(shù)的76.2%，與上年同期相比，地質(zhì)災(zāi)害發(fā)生數(shù)量、造成死亡失蹤人數(shù)和直接經(jīng)濟損失均有所增加。及時快速獲取滑坡災(zāi)害的相關(guān)信息，實現(xiàn)滑坡災(zāi)害預(yù)防、監(jiān)測、治理及災(zāi)后評估，對災(zāi)中有效地組織救援工作、最大限度地挽救人民生命、降低經(jīng)濟損失具有重要的意義。在滑坡災(zāi)害中，快速高效地獲取滑坡信息是一切工作的基礎(chǔ)，對防災(zāi)、減災(zāi)和救災(zāi)發(fā)揮著至關(guān)重要作用。衛(wèi)星遙感技術(shù)具有快捷高效、大范圍對地觀測的優(yōu)勢，已成為一種有效的自然災(zāi)害監(jiān)測手段。根據(jù)傳感器光譜范圍不同，可分為光學遙感和微波遙感兩類。由于滑坡災(zāi)害發(fā)生時往往伴隨陰雨等天氣條件，光學遙感應(yīng)用受到很大限制。合成孔徑雷達（Synthetic Aperture Radar，SAR）作為微波遙感的一種類型，不受云雨霧等自然條件影響，能夠全天時、全天候的獲取數(shù)據(jù)，在滑坡監(jiān)測應(yīng)用中得到了學者的高度重視。Weirich［1］利用多幅航拍的SAR數(shù)據(jù)以及Spot5數(shù)據(jù)進行滑坡敏感制圖比較和研究。Noferini［2］利用SAR數(shù)據(jù)進行滑坡災(zāi)害的緩慢變形的監(jiān)測。Chihira等［3］人對汶川地震發(fā)生嚴重地形形變的區(qū)域進行滑坡判斷，并分析了滑坡的分布情況。VanWestern等［4］人以在雷達圖像等多種衛(wèi)星數(shù)據(jù)的基礎(chǔ)上，實現(xiàn)了汶川地震區(qū)域內(nèi)共60000處地震滑坡的的分布圖。相比國外，我國對于利用SAR圖像提取滑坡信息的研究起步較晚，成功開展應(yīng)用研究的案例不多，數(shù)據(jù)處理過程復雜。因此，針對滑坡災(zāi)害的應(yīng)急監(jiān)測，如何利用SAR數(shù)據(jù)快速及時提取滑坡體信息成為目前廣泛關(guān)注的焦點。本文根據(jù)滑坡體快速提取的需求，以滑坡體等地物的雷達后向散射原理為基礎(chǔ)，提出雷達比值指數(shù)法，并結(jié)合災(zāi)前和災(zāi)中SAR圖像統(tǒng)計的滑坡區(qū)及非滑坡區(qū)的實際雷達特征，分析雷達比值指數(shù)的合理性，最后，針對單景災(zāi)中SAR圖像，運用雷達比值指數(shù)法開展滑坡體的快速提取。

2 研究區(qū)與數(shù)據(jù)

2.1 四川省茂縣區(qū)概況茂縣位于青藏高原向川西平原過度地帶，屬四川省阿壩藏族羌族自治州，處于岷江上游中段。地形以高山峽谷為主，地勢西北高、東南低；水系為岷江水系，其中岷江由北向南貫穿全境。境內(nèi)自然地貌特征呈現(xiàn)高山林立、河流深切。茂縣地層以馬爾康分區(qū)地層為主，主要有千枚巖、砂巖、板巖、白云巖、碳酸鹽、泥頁巖分布。境內(nèi)茂汶斷裂、龍門山斷裂地質(zhì)活動較頻繁，巖石節(jié)理深度發(fā)育，地質(zhì)條件脆弱。氣候干燥多風，降雨年內(nèi)分布不均，多集中在5-10月，最大日降水量可達75.2 mm，瞬時降雨大，導致茂縣地區(qū)滑坡、崩塌和泥石流等地質(zhì)災(zāi)變現(xiàn)象普遍發(fā)生。2017年6月24日5點45分，四川省阿壩州茂縣疊溪鎮(zhèn)新磨村新村組富貴山山體突發(fā)高位垮塌，該滑坡具體位于岷山南部的富貴山西坡，夾于松坪溝河谷和岷江河谷之間［5］，持續(xù)時間100s。茂縣滑坡是2008年“5·12”汶川地震后發(fā)生的最大規(guī)模的滑坡，垮塌方量約1000萬立方米，有62戶100余人被掩埋，堵塞河道約2公里。

2.2 Sentinel-1 SAR影像及預(yù)處理Sentinel-1星座以歐洲及世界其他地區(qū)的自然災(zāi)害和環(huán)境事故快速監(jiān)測為主要目標，由歐洲委員會（Council of Europe，UC）出資、歐洲航天局（European Space Agency，ESA）設(shè)計與開發(fā)，位于同一軌道平面內(nèi)軌道相位差為180°的A和B兩顆衛(wèi)星組成，屬于歐洲哥白尼計劃（Copernicus Programme；又稱Global Monitoring for Environment and Security，GMES，全球環(huán)境與安全監(jiān)測計劃）的中重要組成部分。Sentinel-1A于2014年4月3號發(fā)射，單顆星軌道周期為12 d。后續(xù)A、B兩顆星常規(guī)重復周期為6 d，高緯度地區(qū)最快重訪時間為1-3 d。Sentinel-1A/B搭載一臺C波段（5.404GHz）合成孔徑雷達（Synthetic Aperture Radar，SAR），包含4種成像模式：SM（Stripmap）、IW（Interferometric Wide swath）、EW（Extra-Wide swath）和 WV（Wave）模式，并且 SM，IW 和EW模式分為單極化（HH/VV）和雙極化（HH+HV/VV+VH）成像，最大覆蓋幅寬可達400km［6］。其中IW模式影像為Sentinel-1的主要模式拍攝數(shù)據(jù)，采用了TOPSAR（Terrain Observation with Prigressive Scanning SAR）成像技術(shù)，由3個條帶拼接而成，地距分辨率為20m×22m，覆蓋幅寬為375 km。

通過利用2017年6月24日滑坡災(zāi)害后單景IW模式、GRD類型的VH和VV雙極化方式SAR影像，對茂縣地區(qū)的滑坡體進行快速提取。針對Sentinel-1A SAR數(shù)據(jù)，利用歐空局發(fā)布的SNAP2.0軟件進行數(shù)據(jù)預(yù)處理，主要包括定標、熱噪聲去除、軌道糾正、濾波、相互配準和正射校正等，最后生成空間分辨率為20m×20m的雷達后向散射系數(shù)數(shù)據(jù)。有關(guān)SAR數(shù)據(jù)介紹如表1所示。

表1 Sentinel-1A SAR數(shù)據(jù)

2.3 其他數(shù)據(jù)除了SAR數(shù)據(jù)外，為了對Sentinel-1數(shù)據(jù)進行正射校正，本研究收集了研究區(qū)30 m分辨率ASTER DEM數(shù)據(jù)。且為了對比分析滑坡體提取結(jié)果，收集了災(zāi)后無人機航拍影像。

3 滑坡體雷達遙感監(jiān)測方法

3.1 雷達后向散射機理森林與裸土作為自然界的兩種典型地物類型，與雷達入射波的相互作用具有不同的散射機理。由于入射雷達波性質(zhì)、植被區(qū)地形、植被的覆蓋度等多種因素影響，雷達波與植被相互作用的機理非常復雜。森林區(qū)雷達波散射既包含冠層頂部的面散射，也包含經(jīng)冠層衰減的地面反射的面散射和穿透植被冠層后樹干與地面之間的二次散射。此外，還有更加復雜的植被冠層體散射回波［7］。與裸土表面相比，森林冠層的體散射具有較強的多次散射和去極化過程。其中，去極化過程是指散射的雷達波的極化方向與入射波極化方向相比發(fā)生偏轉(zhuǎn)，極化特性改變，呈現(xiàn)交叉極化回波強度的增強的現(xiàn)象。裸土區(qū)土壤的散射特性與雷達波參數(shù)（入射角、頻率和極化方式）、土壤的含水量及土壤的粗糙度等因素相關(guān)，主要呈現(xiàn)為粗糙面散射。裸土區(qū)粗糙面散射以單次散射為主，多次散射為輔，雷達回波的去極化現(xiàn)象相對較弱。并且隨著粗糙度的增加，后向的單次散射和多次散射增大，土壤粗糙度與后向散射系數(shù)間存在近似指數(shù)的依賴關(guān)系［8-9］，這種依賴關(guān)系隨著入射角的增加更為顯著［10-11］；同時裸土區(qū)的雷達后向回波強度與土壤含水量存在正比關(guān)系，隨著土壤含水量的增加，裸土區(qū)后向散射系數(shù)增大［12］。

3.2 基于雷達比值指數(shù)的滑坡體提取方法根據(jù)植被與裸土區(qū)的散射原理分析可知，與植被區(qū)相比，裸露地表因無植被冠層覆蓋，體散射較少，雷達回波的去極化現(xiàn)象較弱。一般裸土區(qū)的VH極化后向散射強度較小，VV極化散射強度一般較大，且不存在因植被冠層遮擋引起的地面散射衰減影響。但是，在植被覆蓋的山坡發(fā)生滑坡時，滑坡區(qū)土層裸露，地表粗糙度和土壤含水量增加，VH極化和VV極化后向散射系數(shù)可能同時增大，降低了植被區(qū)與非植被區(qū)的單極化后向散射強度差異。因此，在山體滑坡區(qū)，相比單極化方式的VH或VV后向散射系數(shù)，VV極化與VH極化的后向散射系數(shù)比值應(yīng)該更能突出裸露地表與植被覆蓋區(qū)的雷達圖像差異特征。為此，針對Sentinel-1雷達C波段的VV和VH雙極化數(shù)據(jù)，提出基于雷達比值指數(shù)（Ratio Index，RI）的滑坡體快速監(jiān)測法。雷達比值指數(shù)RI，計算公式如下：

式中：σVV為VV極化后向散射系數(shù)，linear形式；σVH為VH極化后向散射系數(shù)，linear形式。

4 茂縣滑坡體散射特征分析

在茂縣滑坡災(zāi)害發(fā)生前后，滑坡區(qū)的地表特征發(fā)生變化，由森林覆蓋轉(zhuǎn)為粗糙裸露地表。在滑坡前滑坡區(qū)雷達散射為森林冠層表層面散射、體散射及經(jīng)冠層散射的地表面散射，并且森林冠層的面散射和體散射對后向散射強度貢獻最大。在滑坡發(fā)生后，滑坡區(qū)的植被覆蓋消失，呈現(xiàn)裸露地表，且地表粗糙度增大，對后向散射強度的貢獻主要為裸土地表的粗糙面散射。除此之外，滑坡災(zāi)害發(fā)生時常伴隨持續(xù)降水，土壤含水量的增加進一步增大滑坡區(qū)的雷達回波強度的變化。因此，為了說明基于雷達比值指數(shù)滑坡體監(jiān)測法的合理性，針對Sentinel-1 SAR數(shù)據(jù)開展滑坡體的后向散射特征分析。

由于2017年6月12日與2017年6月24日前后都伴隨降雨，因此選取與滑坡區(qū)植被和地形坡度類似的非滑坡區(qū)作為對比，分別統(tǒng)計6月12日與6月24日非滑坡區(qū)與滑坡區(qū)兩類統(tǒng)計區(qū)的VH極化、VV極化后向散射系數(shù)的差異，如表2所示。

表2 災(zāi)前與災(zāi)中Sentenel-1 SAR圖像后向散射特征統(tǒng)計

由表2可知，對于同類型統(tǒng)計區(qū)、同時相SAR數(shù)據(jù)，VV極化后向散射系數(shù)均大于VH極化后向散射系數(shù)，且以6月24日災(zāi)中滑坡區(qū)SAR數(shù)據(jù)不同極化的差異性最為顯著。另外對于非滑坡區(qū)同極化不同時相SAR數(shù)據(jù)，6月24日VH極化、VV極化后向散射系數(shù)平均值僅分別略大于6月12日VH極化、VV極化后向散射系數(shù)，且差異較??；而對于滑坡區(qū)，6月24日VH極化后向散射系數(shù)平均值略大于6月12日VH極化后向散射系數(shù)，6月24日VV極化后向散射系數(shù)平均值明顯大于6月12日VV極化后向散射系數(shù)。根據(jù)公式（2）可以計算不同時相SAR數(shù)據(jù)同極化后向散射系數(shù)變化幅度，計算結(jié)果見表3。

式中，pq為極化方式。

表3 針對災(zāi)前與災(zāi)中SAR數(shù)據(jù)的不同類型統(tǒng)計區(qū)的同極化后向散射系數(shù)變化幅度

由表3可知，對于災(zāi)前和災(zāi)中兩時相的SAR同極化數(shù)據(jù)，滑坡區(qū)后向散射系數(shù)變化幅度明顯大于非滑坡區(qū)。進一步分析可知，相對災(zāi)前，災(zāi)中非滑坡區(qū)兩種不同極化方式后向散射系數(shù)增加幅度都不足10%，表明在森林冠層覆蓋的非滑坡區(qū)，降雨引發(fā)的土壤水含量增大對非滑坡區(qū)總后向散射值影響很小。因此，相對災(zāi)前，災(zāi)中非滑坡森林區(qū)后向散射特征未發(fā)生顯著變化；而相對災(zāi)前，災(zāi)中滑坡區(qū)VV極化后向散射系數(shù)增加幅度明顯大于VH極化后向散射系數(shù)增加幅度。VH極化后向散射系數(shù)增加而幅度較小，表明盡管滑坡發(fā)生后森林冠層的減少甚至消失，會使VH極化總后向散射系數(shù)減小。但降雨引發(fā)的土壤水含量增大以及滑坡后土壤粗糙度增大，對后向散射系數(shù)復合貢獻略大于森林冠層減少的貢獻。綜上分析，森林冠層的減少、土壤粗糙度和含水量增大均致使VV極化后向散射系數(shù)增大，而VH極化后向散射系數(shù)影響復雜部分情況下增大、部分情況下減少。因此，從地物的散射特征的角度，VV極化后向散射系數(shù)與VH極化后向散射系數(shù)比值更能突出滑坡區(qū)的圖像特征。

5 茂縣滑坡體監(jiān)測

2017年6月24日5點45分，四川省阿壩藏族羌族自治州茂縣疊溪鎮(zhèn)新磨村富貴山山體頂部震裂山體突然滑動，高速撞擊下方坡體，并沿程鏟刮坡積塊石土層，碎裂解體并轉(zhuǎn)化為碎屑流，協(xié)同下部老滑坡堆積體，形成擴散型碎屑流引發(fā)高位順層巖質(zhì)滑坡。利用Sentinel-1數(shù)據(jù)，采用雷達比值指數(shù)監(jiān)測方法，對四川省阿壩藏族羌族自治州茂縣新磨村地區(qū)滑坡災(zāi)害進行了應(yīng)急監(jiān)測。在滑坡災(zāi)害發(fā)生前，滑坡體及周圍地區(qū)植被完全覆蓋；在災(zāi)害發(fā)生后，滑坡體地區(qū)高位垮塌，植被覆蓋摧毀，伴隨著降雨影響的裸露土壤、巖石碎屑物的后向散射特征體現(xiàn)出來。與災(zāi)前相比，災(zāi)后的滑坡體區(qū)雷達散射特征改變，滑坡體的雷達比值指數(shù)增大，雷達比值指數(shù)結(jié)果對比如圖1（c）（d）所示，（c）圖表示災(zāi)后雷達比值指數(shù)圖，（d）圖表示災(zāi)前雷達比值指數(shù)圖。圖1（a）（b）分別表示災(zāi)后VV、VH極化后向散射強度圖，滑坡區(qū)均未有顯著的體現(xiàn)。圖1（a）（b）（d）分別與圖1（c）對比可以看出，無論相對單極化的雷達圖像或災(zāi)前雷達比值指數(shù)圖，災(zāi)后雷達比值指數(shù)圖表征的滑坡體特征更突出，高亮度顯示。

為了從圖像特征角度評價雷達比值指數(shù)對滑坡體與非滑坡體的識別效果，采用M.Shimada等［13］提出的分離參數(shù)R來分析雷達比值指數(shù)的優(yōu)勢和合理性，其中分離參數(shù)R計算公式如下：

式中：R為分離參數(shù)，R值越大說明分離效果越明顯；μ1為6月24日（滑坡后）目標參數(shù)在森林區(qū)的平均值；μ2為6月24日（滑坡后）目標參數(shù)在滑坡區(qū)的平均值；σ1為6月24日（滑坡后）目標參數(shù)在森林區(qū)的標準差；σ2為6月24日（滑坡后）目標參數(shù)在滑坡區(qū)的標準差。

表5表明，與雷達單極化的數(shù)據(jù)相比，雷達比值指數(shù)RI更能突出滑坡體與周圍森林區(qū)的差異，以此對滑坡體進行快速提取更為合理。最后，根據(jù)多次試驗確定將RI≥8的部分作為滑坡體范圍。另外，提取的茂縣滑坡體范圍與四川測繪信息局的茂縣滑坡無人航拍圖對比結(jié)果，如圖2所示，兩者幾何形狀相似且范圍基本吻合，滑坡體長度為2378 m，最大滑坡寬度為1135 m，滑坡體面積為1.3 km2（無人機航拍數(shù)據(jù)，空間分辨率0.1 m，監(jiān)測的滑坡體面積為1.62 km2，滑坡體長度約為2658 m，最大滑坡寬度為 1200 m）［14］。

圖1 利用雷達比值指數(shù)的滑坡體區(qū)Sentinel-1 SAR圖像對比（GCP1、GCP2、GCP3表示滑坡體所在位置）

表5 災(zāi)前與災(zāi)中滑坡區(qū)的Sentenel-1 SAR圖像后向散射特征統(tǒng)計

6 結(jié)論

本文以森林和裸土等地物的雷達后向散射機理分析為基礎(chǔ)，提出基于VV極化和VH極化后向散射系數(shù)的比值指數(shù)法，根據(jù)單景2017年6月24日災(zāi)中Sentinel-1 SAR圖像，運用雷達比值指數(shù)法開展滑坡體的快速提取和災(zāi)害監(jiān)測。監(jiān)測和研究結(jié)果顯示：（1）基于Sentinel-1 SAR數(shù)據(jù)快速提取的茂縣滑坡體結(jié)果與四川測繪信息局的茂縣滑坡無人航拍圖監(jiān)測結(jié)果基本一致，新磨村嚴重受災(zāi)。（2）與航拍圖相比，通過精度分析可知，滑坡體面積相差19.7%，滑坡體長度相差10.53%，最大滑坡寬度相差5.42%，總體精度均在80%以上，滿足快速提取精度要求。（3）研究結(jié)果表明：雷達比值指數(shù)能夠顯著的突出滑坡體的圖像特征，該方法在針對Sentinel-1 SAR的滑坡提取與監(jiān)測上具有較大的應(yīng)用潛力和優(yōu)勢。（4）采用本文提出的滑坡體提取方法可快速獲取滑坡體信息，為搶險救災(zāi)和災(zāi)后評估提供信息支持。

圖2 機載航拍的光學圖像的滑坡體與Sentinel-1 SAR圖像提取的滑坡體對比