基于SBAS-InSAR技術的深切割高山峽谷區(qū)滑坡災害早期識別

周定義，左小清，喜文飛,2，肖 波,3，畢 瑞，范 馨

（1.昆明理工大學國土資源工程學院，云南 昆明 650093；2.云南師范大學地理學部，云南 昆明650500；3.云南交通職業(yè)技術學院公路學院，云南 昆明 650500；4.成都理工大學地球科學學院，四川 成都 610059）

0 引言

中國是地質災害高發(fā)的國家，各類地質災害給人們的生命財產(chǎn)造成了巨大的損失[1]。滑坡作為一種主要的地質災害，具有隱蔽性強、危害性大、突發(fā)性高等特點，廣泛分布于中國山區(qū)和峽谷地帶[2?4]。近年來，高山峽谷滑坡頻頻發(fā)生。因此，對高山峽谷地區(qū)的滑坡災害開展早期識別，能夠為防災減災事業(yè)及政府部門決策提供一種有效的手段。

地表形變是反應當前坡體穩(wěn)定性及運動狀態(tài)最直接的物理量，因此，監(jiān)測地表形變可以為探測滑坡等地質災害的隱患點提供重要信息[1]。傳統(tǒng)的地表形變監(jiān)測方法主要采用精密水準測量、全球定位系統(tǒng) (GPS)等。但這些監(jiān)測方法存在變形監(jiān)測工作量大、費時、花費大、測點難以保存等缺陷，同時，傳統(tǒng)的監(jiān)測方法無法對大區(qū)域滑坡形變進行探測[5?6]。與傳統(tǒng)地表形變相比，合成孔徑雷達干涉測量 (Interferometric Synthetic Aperture Radar,InSAR) 技術是近幾年發(fā)展起來的一種新型大地測量手段。具有覆蓋范圍廣、穿透云層、全天候運作、精度高的特點，理論上可以獲得非常精確的數(shù)字高程模型和毫米量級的地表形變信息，已被成功用于滑坡災害監(jiān)測[7?9]。目前，利用InSAR技術對滑坡早期識別的研究取得了一些成功范例：張路等[2]采用自主研發(fā)的相干散射體時序InSAR 方法，成功識別出了17 處持續(xù)變形中的不穩(wěn)定坡體。Dai 等[10]利用InSAR 對潛在的滑坡進行早期識別，10 個潛在滑坡區(qū)域被早期發(fā)現(xiàn)，結果表明，利用InSAR 可以作為獲取滑坡早期識別的有效手段。韓守富等[11]證明了InSAR技術在黃土高原地區(qū)地質災害早期識別方面的適用性和準確性,可以應用于黃土高原地區(qū)地質災害隱患識別預警。馮文凱等[12]利用SBAS- InSAR技術對金沙江流域沃達村滑坡進行地表形變監(jiān)測，表明SBAS- InSAR技術在復雜山區(qū)地質災害監(jiān)測預警領域有較為廣闊的應用前景，為類似老滑坡監(jiān)測預警提供了新的思路與借鑒。戴可人等[13]利用時間序列InSAR技術對雅礱江流域雅江縣—木里縣段的高山峽谷區(qū)域進行了滑坡災害隱患廣域早期識別，成功探測到8 處隱患區(qū)域。以上方法能夠有效識別滑坡災害，但對于深切割高山峽谷區(qū)的滑坡早期識別，僅利用SAR 單軌道數(shù)據(jù)監(jiān)測會致使SAR 成像幾何畸變造成部分滑坡不能識別，只有通過升降軌數(shù)據(jù)結合的方式才能全面準確的對滑坡災害進行早期識別。

文章利用SBAS-InSAR技術，采用升降軌數(shù)據(jù)結合的方式對深切割高山峽谷區(qū)的滑坡進行早期解譯識別，并根據(jù)識別結果,對不同潛在滑坡類型進行了分析與討論，為高山峽谷地區(qū)防災減災及滑坡災害早期識別提供一種更為全面的方法。

1 SBAS-InSAR技術簡介

小基線集（SBAS-InSAR）是在差分InSAR 基礎上發(fā)展起來的一種新的時間序列分析方法，能夠降低相位噪聲和誤差[14]。

假定在時間t1至ts內 獲取同一地區(qū)的S幅SAR 影像，然后根據(jù)干涉組合條件，在短基線距的條件下形成N幅干涉條紋圖，且滿足：

對tA和tB（tA

式中：r——斜距；

假定tk時刻和tk+1時刻不同干涉圖間的形變速率為vk,k+1，則tA—tB間的累積形變可表示為式：

對N幅干涉條紋圖進行三維時空相位解纏即可求出不同SAR 獲取時間的形變速率。

2 研究方法

文中利用SBAS-InSAR技術對深切割高山峽谷區(qū)滑坡災害早期識別，方法的重點有兩個關鍵：①針對深切割高山峽谷區(qū)的滑坡如何識別；②如何保證識別結果的準確性。

2.1 升降軌數(shù)據(jù)結合的方式對深切割高山峽谷區(qū)的滑坡識別

對于深切割高山峽谷區(qū)的滑坡早期識別，由于地形條件等因素的影響，僅利用SAR 單軌道數(shù)據(jù)監(jiān)測會致使SAR 成像幾何畸變造成部分滑坡不能識別。為此，在技術上采用升降軌數(shù)據(jù)結合的方式，對深切割高山峽谷區(qū)滑坡災害進行早期識別。彌補了SAR 單軌道數(shù)據(jù)識別存在的不全面、不準確等弊端。

2.2 引入高分辨率光學影像等作為輔助識別保證準確性

為保證識別結果的準確性，僅憑形變監(jiān)測結果無法區(qū)分是否為潛在的滑坡導致形變，為此，引入高分辨率光學影像對形變區(qū)域進行輔助識別，結合形變范圍、高程、坡度、植被覆蓋和坡體是否具有滑坡特征等進行識別，避免過度依賴形變結果導致的誤判等問題。同時高分辨率光學影像能起到驗證識別結果的作用，早期的滑坡一般會有一定的滑動痕跡，但不代表所有早期滑坡都有滑動痕跡。綜上，引入高分辨率光學影像等作為輔助識別能有效保證識別結果的準確性。

3 滑坡災害早期識別試驗

3.1 研究區(qū)概況和數(shù)據(jù)來源

東川小江流域為世界典型暴雨泥石流區(qū)，被稱為“泥石流的天然博物館”，該地區(qū)泥石流的發(fā)生往往與滑坡密切相關，是屬于暴雨型滑坡泥石流，即先滑坡，再經(jīng)暴雨沖刷后形成泥石流[15]。文中選取小江沿線兩側高山峽谷作為研究區(qū)（圖1），以小江為界，河谷凹陷，形成“V”字型，屬于典型的深切割高山峽谷，東側為牯牛寨山，最高峰海拔4 017.3 m；西部為拱王山，最高峰海拔4 344.1 m。該區(qū)域地勢陡峭，其獨特的地形和地質構造，致使局部區(qū)域暴雨多、土質松軟、水土流失嚴重，導致該區(qū)域內地質災害頻發(fā)。

圖1 研究區(qū)位置Fig.1 Location of study area

研究數(shù)據(jù)是從歐空局（European Space Agency,ESA）下載的40 景C 波段Sentinel-1A 升降軌影像，時間跨度為2018年4月25日—2019年4月20日，極化方式為VV，成像模式為IW。數(shù)據(jù)參數(shù)如表1所示。為了提高影像軌道精度，引入了POD 精密定軌星歷數(shù)據(jù)，使用JAXA（日本宇宙航空研究開發(fā)機構）提供的ALOS WORLD 3D 30 m 空間分辨率的數(shù)字高程模型（Digital Elevation Model,DEM），用于去除地形相位影響。拼接后的DEM 如圖2所示。

圖2 研究區(qū)DEMFig.2 Digital elevation model of study area

表1 Sentinenl-1A 數(shù)據(jù)參數(shù)Table 1 Sentinenl-1A data parameters

3.2 地表形變監(jiān)測試驗

選取時間跨度2018年4月25日—2019年4月20日的Sentinel-1A 斜距單視復數(shù)（Single Look Complex，SLC）影像，其中升降軌數(shù)據(jù)各20 景。利用SARscape軟件進行處理，升軌和降軌數(shù)據(jù)選取日期為2018-08-11和2018-08-13 的影像作為超級主影像，通過設置臨界基線和時間基線，生成103 和100 對像對，設置多視數(shù)為1∶4 可以較好地抑制斑點噪聲，采用Minimum Cost Flow 解纏方法和Goldstein 濾波方法做干涉工作流，最終生成干涉圖，調整刪除不理想的數(shù)據(jù)，在研究區(qū)生成較為理想的部分升降軌干涉圖如圖3所示。從圖3 可以看出，針對山區(qū)地區(qū)，運用升降軌數(shù)據(jù)處理后得到的干涉圖相干性較為理想。

圖3 部分升降軌干涉圖Fig.3 Part of the lifting rail interference diagram

經(jīng)過軌道精煉和重去平，估算和去除殘余的恒定相位和解纏后還存在的相位坡道，進行第一次反演、第二次反演，最后對序列信息進行地理編碼獲得研究區(qū)域2018年4月25日—2019年4月20日雷達視線方向（LOS）的形變速率圖。如圖4所示，正值表示形變朝著衛(wèi)星的方向，負值表示形變沿著遠離衛(wèi)星的方向。圖4（a）為采用升軌數(shù)據(jù)獲取的地表形變速率圖，圖4(b)表示采用降軌數(shù)據(jù)獲取的地表形變速率圖。從圖4 可以看出，采用升軌獲取的形變主要分布在小江東側，最大形變速率為?162.074 mm/a，降軌獲取的形變主要分布在小江西側西北方向，最大形變速率為?120.425 mm/a。致使不同軌道得到的形變結果不同，是由于升軌數(shù)據(jù)飛行方向大致從南到北，雷達視線（LOS）方向位于右側，能夠很好地將峽谷兩側由西向東的地表形變監(jiān)測出來，相反，降軌數(shù)據(jù)的飛行方向與之相反，能夠將峽谷兩側由東向西的地表形變監(jiān)測出來。為此，利用不同軌道數(shù)據(jù)可以互補，使得監(jiān)測結果更為準確全面，能夠避免單一軌道帶來的幾何畸變等問題。

利用InSAR技術在植被覆蓋地區(qū)進行地形測量非常困難，這主要是因為電磁場和（或）散射體的物理特性隨時間的變化而造成的，植被覆蓋度高地區(qū)，失相干嚴重，導致監(jiān)測地面形變信息的能力較差，對形變監(jiān)測精度的影響較大[16?17]。為考慮監(jiān)測形變結果的有效性，引入歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）對研究區(qū)植被覆蓋數(shù)據(jù)進行分析。如圖5所示，NDVI 結果被限定在[?1,1]，負值表示地面覆蓋為云、水、雪等，對可見光高反射，0 表示有巖石或裸土等，正值表示有植被覆蓋，且隨覆蓋度增大而增大。對比圖4 和圖5 可以發(fā)現(xiàn)，升降軌獲取的部分形變區(qū)處于植被覆蓋較高的地區(qū)，這正是由于植被造成失相干嚴重，致使這些區(qū)域監(jiān)測結果不準確，在分析和識別滑坡時，需把其剔除，避免導致滑坡識別的誤判。

圖4 升降軌獲得的地表形變速率圖Fig.4 Surface deformation rate map obtained by lifting rail

圖5 研究區(qū)NDVIFig.5 NDVI in the study area

3.3 滑坡風險區(qū)的總體識別

深切割高山峽谷地區(qū)同樣存在非滑坡原因的沉降變形，為保證獲取的形變信息為滑坡災害，僅憑形變監(jiān)測結果無法區(qū)分是否為潛在的滑坡導致形變，為此，引入高分辨率光學影像對形變區(qū)域進行輔助識別，結合形變范圍、高程、坡度、植被覆蓋和坡體是否具有滑坡特征等進行識別，避免過度依賴形變結果導致的誤判等問題。由于不同軌道獲取的形變數(shù)據(jù)能夠對研究區(qū)不同方向的滑坡進行早期識別，分別對不同軌道獲取的數(shù)據(jù)結合光學影像解譯，升軌數(shù)據(jù)解譯結果如圖6所示，黑色區(qū)域為滑坡形變區(qū)域，紅色虛線區(qū)域為對應光學影像滑坡災害隱患，黃色箭頭代表滑坡方向,共解譯出11 處潛在滑坡，分別用H01、H02、······、H11進行編號，其中H01、H02、H08 和H09 位于現(xiàn)存泥石流溝處，H01 位于尖山溝老尖山、土巖子一帶，H02 位于蔣家溝一帶，H08、H09 位于大白泥溝和小白泥溝一帶，其他潛在滑坡有一定的滑坡痕跡。識別出來的滑坡區(qū)域平均最大形變速率超過?40 mm/a，最大形變速率?120 mm/a，位于H06 巖子腳村上方。根據(jù)滑坡是否直接威脅附近周圍村落，對11 處潛在滑坡進行危險劃分，共識別出4 處高風險滑坡，分別是H06、H09、H10 和H11，H06 坡體對巖子腳村構成威脅，H09 坡體對大村和魯納窩村構成威脅，H10 對小多紅村和妥托村構成威脅，H11 對小凹子村、廠上村、大麥地、小村子等構成威脅，這些區(qū)域一旦發(fā)生滑坡，極有可能對當?shù)厝罕娫斐缮柏敭a(chǎn)威脅。

圖6 升軌數(shù)據(jù)滑坡災害識別結果Fig.6 Landslide disaster identification results of lifting data

降軌數(shù)據(jù)解譯結果如圖7所示，共解譯出7 處潛在滑坡，分別用H12、H13、······、H18進行編號，滑坡區(qū)域平均形變速率均超過?40 mm/a，最大形變區(qū)域位于H18 大白泥溝和小白泥溝連接處，其最大形變速率為?94.11 mm/a。根據(jù)滑坡是否直接威脅附近周圍村落，對7 處潛在滑坡區(qū)進行危險劃分，共識別出1 處滑坡高風險區(qū)，編號為H16，該坡體對姑莊村和新寨田村構成威脅。對比圖6 和圖7 發(fā)現(xiàn)，不同軌道識別出來的潛在滑坡分布不同，利用升軌獲取的滑坡主要分布在小江右側峽谷中部，降軌獲取的滑坡主要分布在小江左側西北方向，其中不同軌道識別到的潛在滑坡有相同位置，但形變區(qū)域不一致，例如，H01、H13 共處于尖山溝一帶，H02、H15 共處于蔣家溝一帶，H08、H18 共處于大小白泥溝一帶。由此證明，利用升降軌結合的方式能夠有效識別深切割高山峽谷地區(qū)不同坡度方向存在的潛在滑坡，避免了單一軌道存在的識別結果不準確，不全面等問題。

圖7 降軌數(shù)據(jù)滑坡災害識別結果Fig.7 Landslide disaster identification results of rail descent data

3.4 典型災害滑坡分析

經(jīng)過對比光學影像發(fā)現(xiàn)，研究區(qū)識別出來的潛在滑坡可定義為三種類型，分別是處于泥石流區(qū)域的潛在滑坡（H09）、有滑坡痕跡的潛在滑坡（H06）和無滑坡痕跡的潛在滑坡（H16）。為有效了解不同類型潛在滑坡的形變趨勢，分別對三種典型的潛在滑坡進行分析。

3.4.1 H09 典型滑坡

H09 屬于典型的泥石流區(qū)域潛在滑坡，位于小白泥溝，該區(qū)域為暴雨型泥石流多發(fā)區(qū)。利用獲取的形變速率值與三維光學影像疊置分析，其形變速率如圖8所示，存在A、B 和C 三個潛在滑坡區(qū)，呈V 字型分布，滑動方向由坡面向下滑動，在每個潛在滑坡上方能明顯的觀察到斷裂滑動痕跡，最大形變速率為?103.013 mm/a，位于潛在滑坡A 面上方。從圖8 中可以看到，B、C 面分別位于小白泥溝兩側，其潛在滑坡屬于長期滑坡，潛在滑坡上端形變速率大，下端由于滑坡堆積物導致一定的抬升。但A 面潛在滑坡還未形成真正意義的滑坡，一旦A 面滑坡，極有可能威脅大村和魯納窩村民的生命及財產(chǎn)安全。為更有效全面地分析該種類型潛在滑坡的特點，選取A 面進行詳細分析，可以看到A 面呈現(xiàn)三處不均勻形變，有明顯的拉張裂縫且有滑坡痕跡，其中P1 位于整個潛在滑坡頂端，P2 位于潛在滑坡中部，P3 位于潛在滑坡右側。三個點均有可能在后期的降雨等因素下發(fā)展為新的潛在滑坡區(qū)。圖9 為3 處潛在滑坡特征點形變量與當月降雨量的關系圖，可以看到特征點的形變量與降雨量有一定的相關關系。研究區(qū)總體位于小江斷裂帶沿線，該區(qū)域斷裂帶寬5～20 km，呈現(xiàn)弱剪切強擠壓活動特征[18]，以擠壓穹起隆升變形為主[19]，從形變時序圖9 可以看出，三個特征點在2018年5月先抬升后沉降，抬升可能是由于斷裂帶擠壓活動致使，后經(jīng)過雨水沖刷導致后期慢慢滑動沉降。P1、P2 和P3 點的平均形變量分別為?10.571 mm、?22.564 mm、?19.516 mm，年形變量分別為?49.063 mm、?58.740 mm、?62.635 mm。該區(qū)域月均降雨量為102 mm。隨著時間的推移，形變量逐漸增加，到一定臨界值便形成真正意義上的滑坡。

圖8 H09 潛在滑坡形變速率圖Fig.8 Deformation rate diagram of H09 potential landslide

圖9 P1—P3 時間序列曲線與降雨量Fig.9 Time series curve of P1—P3 and rainfall

3.4.2 H06 典型滑坡

H06 屬于典型的有滑坡痕跡的潛在滑坡，位于牯牛山巖子腳上方。圖10 為該潛在滑坡形變速率與三維影像的疊置圖，從圖10 中可以白色虛線內兩處潛在滑坡有明顯的滑動痕跡，滑動方向為坡頂沿著坡腳移動。最大形變速率位于P5 附近，達到?127.093 mm/a。選取特征點P4、P5 形變量與月降雨量構建關系圖（圖11）。隨著時間推移，形變量不斷增大。特征點的形變量與月降雨量有一定的相關關系，從形變時序圖11 可以看出，兩個特征點在2018年5月先抬升后沉降，該潛在滑坡位于小江斷裂帶沿線，可能是由于斷裂帶擠壓活動致使其抬升。當月降雨量增多加劇潛在滑坡的形變，這是由于降雨致使表層松散土體隨著雨水滾落流失。P4、P5點的平均累積形變量分別為?37.562 mm、?42.054 mm，年形變速率分別為?105.903 mm/a、?112.469 mm/a。該潛在滑坡區(qū)最高形變點高程為3 606 m，地勢陡峭，從年形變量可以看出該區(qū)域存在很大安全隱患，一旦發(fā)生滑坡，巖子腳村村民將面臨極大的生命及財產(chǎn)安全。

圖10 H06 潛在滑坡形變速率圖Fig.10 Deformation rate diagram of H06 potential landslide

圖11 P4—P5 時間序列曲線與降雨量Fig.11 P4—P5 time series curve and rainfall

3.4.3 H16 典型滑坡

H16 屬于典型的無滑坡痕跡的潛在滑坡，位于姑莊村和新寨田上方。圖12 為該潛在滑坡區(qū)形變速率與三維影像的疊置圖，從圖12 中可以看到，位于潛在滑坡上方有明顯的形變，最大形變速率為?80.141 mm/a。通過光學影像看到山體表面并無明顯的滑動痕跡。該潛在滑坡區(qū)最高高程為2 304 m，坡度大，滑動方向大致沿著山溝滑動。選取特征點P6、P7 形變量與月降雨量構建關系圖，如圖13所示，形變與降雨有一定關系。從形變時序圖13 可以看出，P6、P7 兩個特征點總是先抬升后沉降再抬升這一反復過程，在2018年5月可能是由于斷裂帶擠壓活動致使其抬升，其他月份抬升微小可能是由于雨水沖刷導致土地堆積抬升，再經(jīng)雨水沖刷致使微小沉降。P6、P7 平均累積形變量分別為?29.219 mm、?21.281 mm，年形變速率分別為?63.788 mm/a、?53.110 mm/a。兩個特征點形變曲線大致一直，從潛在滑坡上方到中部形變速率逐漸減小。該潛在滑坡一旦發(fā)生滑坡，將直接威脅姑莊村和新寨田兩個村。

圖12 H16 潛在滑坡形變速率圖Fig.12 Deformation rate diagram of H16 potential landslide

圖13 P6—P7 時間序列曲線與降雨量Fig.13 P6—P7 time series curve and rainfall

4 結論

文中基于SBAS-InSAR技術，采用Sentinel-1 升降軌數(shù)據(jù)結合互補的方式對東川小江沿線兩側深切割高山峽谷區(qū)滑坡災害進行早期識別實驗，得出以下結論：

（1）升軌數(shù)據(jù)識別出的滑坡主要分布在小江右側，降軌反之，這表明相比單一軌道獲取的識別結果，利用升降軌結合的方式能夠更全面的監(jiān)測和識別高山峽谷滑坡，避免單軌道對高山峽谷區(qū)滑坡進行早期識別存在SAR 成像幾何畸變造等問題。

（2）在高山峽谷地區(qū)植被覆蓋度過高容易致使失相干，同時，有部分地區(qū)并不屬于滑坡導致的形變，為此，引入光學影像，結合形變范圍、高程、坡度、植被覆蓋和坡體是否具有滑坡特征等進行識別，避免過度依賴形變結果導致的誤判等問題。

（3）通過H09、H06、H16 等3 個典型潛在滑坡災害分析，可以看到滑坡的形成和降雨量有一定的關系。

（4）文中共識別出18 處潛在滑坡區(qū)，其中H06、H09、H10、H11、H16 等5 處為滑坡高風險區(qū)，證明利用該方法可作為高山峽谷區(qū)滑坡災害識別的有效手段，但目前還有一定的不足，比如在高山峽谷地區(qū)SAR 數(shù)據(jù)量過少致使部分滑坡并未獲取到形變量等問題。