徐曉雪,季靈運(yùn),3*,張文婷,劉傳金

(1. 中國地震局第二監(jiān)測(cè)中心，陜西 西安 710054; 2. 國家遙感中心 地質(zhì)災(zāi)害部，北京 100036;3. 防災(zāi)科技學(xué)院 地球科學(xué)學(xué)院，河北 三河 065201)

0 引 言

合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar，InSAR)技術(shù)由于全天時(shí)、全天候、覆蓋范圍廣、監(jiān)測(cè)精度高和空間分辨率高的優(yōu)點(diǎn)，在山區(qū)監(jiān)測(cè)活躍滑坡方面顯示出巨大的潛力。InSAR技術(shù)可以提供滑坡位置、邊界和運(yùn)動(dòng)特征的關(guān)鍵信息，為滑坡的災(zāi)前時(shí)空形變追溯提供了新的技術(shù)途徑。自Achache等首次使用差分合成孔徑雷達(dá)干涉測(cè)量(Different InSAR,D-InSAR)技術(shù)對(duì)法國阿爾卑斯地區(qū)La Clapiere 滑坡進(jìn)行監(jiān)測(cè)以來，為了提高形變測(cè)量的精度和可靠性，在隨后的研究中諸如永久散射體InSAR(PSI)方法、InSAR小基線子集(Small Baseline Subset,SBAS)方法和混合時(shí)序InSAR方法等多種時(shí)間序列InSAR方法被提出來。其中，SBAS方法具有獲取微小形變信息和長時(shí)間序列緩慢地表變形的優(yōu)勢(shì)，能夠在一定程度上克服時(shí)間失相干和空間失相干的不利影響，在復(fù)雜山區(qū)滑坡形變監(jiān)測(cè)方面具有較為廣闊的應(yīng)用前景。很多學(xué)者采用SBAS方法進(jìn)行滑坡的早期識(shí)別與形變監(jiān)測(cè)，如Manunta等使用SBAS方法對(duì)意大利中部翁布里亞地區(qū)進(jìn)行大范圍滑坡識(shí)別與監(jiān)測(cè)；Zhao等使用SBAS方法對(duì)美國西部森林覆蓋區(qū)域的滑坡進(jìn)行調(diào)查，獲得了滑坡的時(shí)間序列形變特征；Liu等采用SBAS方法對(duì)美國加利福尼亞州三只熊滑坡的滑前形變進(jìn)行監(jiān)測(cè)，清晰地厘定了森林山區(qū)緩慢移動(dòng)滑坡的運(yùn)動(dòng)特征。隨著SAR衛(wèi)星的相繼發(fā)射和大量SAR數(shù)據(jù)的存儲(chǔ)，基于雷達(dá)遙感數(shù)據(jù)的SBAS方法在未來將會(huì)發(fā)揮更加重要的作用。

滑坡是中國西南山區(qū)最主要的地質(zhì)災(zāi)害之一，由于地形復(fù)雜多變且滑坡破壞區(qū)域較大，往往難以防范，滑坡的高精度監(jiān)測(cè)對(duì)于了解其運(yùn)動(dòng)學(xué)特征和減少滑坡災(zāi)害具有重要意義。2020年8月21日，四川省雅安市漢源縣富泉鎮(zhèn)中海村6組發(fā)生山體滑坡(簡(jiǎn)稱“漢源滑坡”)，滑坡體總方量達(dá)到500×10m，屬特大型滑坡。此次滑坡造成群眾房屋受損及部分人員傷亡，省道435交通中斷，給人民生命財(cái)產(chǎn)安全帶來威脅。追溯滑坡體災(zāi)前形變狀態(tài)，對(duì)后續(xù)開展針對(duì)性的工程地質(zhì)調(diào)查及滑坡誘因分析具有重要作用。因此，本文基于升軌、降軌哨兵一號(hào)(Sentinel-1)衛(wèi)星影像，求解滑坡體的InSAR時(shí)間序列，研究四川雅安地區(qū)漢源滑坡的災(zāi)前演變情況。由于漢源滑坡發(fā)生在大渡河流域高山河谷地區(qū)，坡體表面作物覆蓋，極大地影響了InSAR數(shù)據(jù)處理中干涉圖的相干性，所以采用基于相干性的小基線集時(shí)間序列算法，以適應(yīng)低相干分析，更好地獲得形變場(chǎng)，從而厘定其時(shí)空變化規(guī)律，為災(zāi)后穩(wěn)定性評(píng)估和未來滑坡災(zāi)害的監(jiān)測(cè)預(yù)警提供思路。

1 研究區(qū)概況

漢源縣位于四川省西南部，隸屬雅安市管轄(圖1中矩形框代表Sentinel-1影像覆蓋范圍)，為四川盆地與青藏高原之間的過渡地帶，省道306經(jīng)縣境西北順著大渡河穿過，與成昆線相接。從地形來看，漢源縣位于大渡河流域高山河谷地區(qū)，河道下切，山高坡陡，是滑坡易發(fā)區(qū)。此次滑坡所在的區(qū)域?qū)俚湫椭械蜕叫逼碌孛?，為臺(tái)階狀地形，斜坡前緣較陡。坡面上既有陡坡又有緩坡，滑坡體整體表面坡度約25°，最陡地方的坡度達(dá)40°，為利于滑坡失穩(wěn)的地形地貌條件。滑坡體整體呈長條形舌狀，寬約280 m，斜長約880 m。

紅色五角星及紅色圓圈為滑坡位置；藍(lán)色方框?yàn)镾AR數(shù)據(jù)覆蓋范圍圖1 四川雅安地區(qū)漢源滑坡位置和SAR數(shù)據(jù)覆蓋Fig.1 Location of Hanyuan Landslide in Ya’an Area of Sichuan and Coverage of the SAR Datasets

漢源縣區(qū)域巖土體可分為松散巖、碎屑巖、碳酸鹽巖、變質(zhì)巖和巖漿巖5種巖類，易滑巖組的存在在降雨誘發(fā)下易于發(fā)生區(qū)域滑坡事件。漢源滑坡發(fā)生后，多部門開展了現(xiàn)場(chǎng)成因調(diào)查。從地質(zhì)結(jié)構(gòu)來看，滑坡體表層的土壤為耕作土和砂性土，結(jié)構(gòu)十分松散；土層下面的巖石是風(fēng)化很強(qiáng)烈的砂泥巖，容易失穩(wěn)變形。漢源縣由于所處的獨(dú)特地理位置，具有年降雨豐沛、降雨集中和連續(xù)多日強(qiáng)降雨的特點(diǎn)。從2020年8月10日起，漢源縣經(jīng)歷了兩輪強(qiáng)降雨天氣，巖土體含水率接近飽和，斜坡整體穩(wěn)定性降低，容易發(fā)生地質(zhì)災(zāi)害。降雨通過表層土壤下滲到巖石，慢慢浸潤巖土體，超過穩(wěn)定臨界值之后就發(fā)生了滑坡。強(qiáng)降雨可能是此次滑坡發(fā)生的誘發(fā)因素。

2 數(shù)據(jù)來源與分析方法

2.1 研究區(qū)SAR數(shù)據(jù)與干涉處理

本文所采用的SAR數(shù)據(jù)為Sentinel-1衛(wèi)星影像。Sentinel-1衛(wèi)星影像不僅可以實(shí)時(shí)免費(fèi)下載，而且覆蓋范圍廣，重訪周期短，軌道精度高，可以進(jìn)行全天候、全天時(shí)的高分辨率監(jiān)測(cè)。由于SAR采用側(cè)視成像方式，而且西南山區(qū)地形復(fù)雜多變，致使覆蓋漢源滑坡區(qū)域的降軌Sentinel-1 SAR影像在滑坡區(qū)域形成陰影，雷達(dá)接收不到后向散射信息，嚴(yán)重影響滑坡監(jiān)測(cè)的準(zhǔn)確度。因此，本文只針對(duì)覆蓋漢源滑坡區(qū)域的升軌Sentinal-1 SAR影像進(jìn)行討論，影像覆蓋范圍如圖1中矩形框所示。本研究所采用的數(shù)據(jù)分別為59景、56景兩個(gè)升軌的Sentinel-1影像數(shù)據(jù)，不同軌道SAR數(shù)據(jù)的時(shí)間跨度統(tǒng)一為滑坡發(fā)生之前的2018年10月至2020年8月，時(shí)間間隔為12 d。使用美國航空航天局(NASA)發(fā)布的30 m空間分辨率SRTM1 DEM地形數(shù)據(jù)模擬并消除地形對(duì)干涉相位的貢獻(xiàn)，與SRTM3數(shù)據(jù)相比具有更高的空間分辨率。

數(shù)據(jù)處理使用開源InSAR處理軟件GMTSAR進(jìn)行。由于此次滑坡的空間范圍較小，在雷達(dá)坐標(biāo)距離向和方位向不進(jìn)行降采樣，多視比設(shè)置為1∶1；經(jīng)過地理編碼后，單個(gè)像素的空間分辨率約為25 m。采用較小的多視比得到的差分干涉圖相干性較低，存在較多斑點(diǎn)噪聲，因此，對(duì)原始差分干涉圖進(jìn)行濾波，以減少噪聲，提高相干性。對(duì)于軌道誤差改正，通過建立二次多項(xiàng)式模型來消除。由于滑坡位于復(fù)雜地貌山區(qū)，茂密的植被使SAR影像難以保持較好的相干性，也很難在植被茂密的滑坡上提供有用信號(hào)，干涉對(duì)在很大程度上受到時(shí)間失相干的影響，特別是臨滑前夏季影像生成的干涉對(duì)失相干現(xiàn)象明顯，導(dǎo)致了長時(shí)間跨度干涉圖的時(shí)間失相干。在去除低相干干涉對(duì)后，最終選擇P128軌道的54個(gè)干涉對(duì)和P26軌道的41個(gè)干涉對(duì)進(jìn)行后續(xù)InSAR時(shí)間序列分析。這些干涉對(duì)的基線通常小于100 m，SAR數(shù)據(jù)的具體獲取日期與時(shí)空基線關(guān)系如圖2所示。

20181216是指2018年12月16日，其他依此類推圖2 SAR數(shù)據(jù)時(shí)空基線圖Fig.2 Spatial and Temporal Baselines of SAR Datasets

2.2 InSAR時(shí)間序列分析

SBAS方法采用較小的時(shí)空基線，在提取滑坡沿雷達(dá)視線方向(Line of Sight,LOS)的時(shí)間序列形變方面較為常用。SBAS方法的主要思路是基于一定的時(shí)空基線閾值，選取多個(gè)主影像組成自由組合的干涉對(duì)，若干涉構(gòu)網(wǎng)形成一個(gè)唯一子集并充分連接，可采用最小二乘法求解時(shí)間序列信息。當(dāng)存在多個(gè)子集時(shí)，依據(jù)各相干像元相位與觀測(cè)時(shí)間的關(guān)系，利用奇異值分解(Singular Value Decomposition,SVD)方法將多個(gè)小基線集合數(shù)據(jù)聯(lián)合起來求解，有效地解決了各數(shù)據(jù)集之間的時(shí)間不連續(xù)問題，從而求得時(shí)間序列形變、平均形變速率結(jié)果。然而，由于研究區(qū)植被覆蓋茂密，傳統(tǒng)SBAS方法容易受到失相干的影響，低相干像元通常使用相干性閾值進(jìn)行掩膜，降低了時(shí)間序列結(jié)果的點(diǎn)位密度。

本文采用基于相干性的SBAS方法進(jìn)行InSAR時(shí)間序列分析，將干涉圖的相位相干性引入到SBAS時(shí)間序列處理算法中，不丟棄干涉圖中存在的部分噪聲數(shù)據(jù)，盡可能多地保留干涉圖中的像素，通過相干性對(duì)觀測(cè)相位進(jìn)行加權(quán)以減弱失相干的影響，提高了時(shí)間序列產(chǎn)品的空間分辨率。在時(shí)間序列處理過程中，通過SNAPHU軟件進(jìn)行相位解纏，選擇較低的相干性閾值0.16，以獲得盡可能多的相位信息，這可能會(huì)帶來相位解纏誤差，但是在時(shí)間序列分析中失相干像元點(diǎn)會(huì)被降權(quán)處理。這種基于相干性的SBAS方法使用加權(quán)矩陣根據(jù)每個(gè)差分干涉圖中每個(gè)像素的相干性對(duì)觀測(cè)相位進(jìn)行加權(quán)，對(duì)差分干涉圖中的噪聲不太敏感，降低了失相干相位的權(quán)重，從而可以得到相干信號(hào)的密集空間覆蓋。每個(gè)像素的加權(quán)最小二乘反演可以表示為

(1)

=diag{,,,,}

(2)

2.3 坡向形變提取

采用InSAR技術(shù)只能得到三維地表變形沿雷達(dá)視線方向的投影，然而大多數(shù)滑坡運(yùn)動(dòng)通常發(fā)生在平行斜坡方向。在一定的假設(shè)條件下，雷達(dá)視線方向變形可以轉(zhuǎn)移到由特定的坡度角和方位角所確定的斜坡方向。本文將雷達(dá)視線方向測(cè)量結(jié)果投影到滑坡的斜坡方向，獲得滑坡沿坡向的運(yùn)動(dòng)情況，平行斜坡運(yùn)動(dòng)可以更好地識(shí)別滑動(dòng)體，提高對(duì)滑坡實(shí)際運(yùn)動(dòng)的量化分析。

經(jīng)過投影轉(zhuǎn)換后，坡向矢量大小總是不小于雷達(dá)視線方向。從雷達(dá)視線方向到坡向的轉(zhuǎn)換可以很容易地通過放大比例因子來實(shí)現(xiàn)。放大比例因子()表達(dá)式為

(3)

式中：=[sinsin-sincoscos]，為雷達(dá)視線方向的形變單位向量；=[-coscos-cos·sinsin]，為坡向的單位向量；和分別表示雷達(dá)飛行方位角和入射角；和分別為坡角和坡向。

3 災(zāi)前形變結(jié)果分析

3.1 雷達(dá)視線方向形變速率與時(shí)間序列

采用基于相干性的InSAR時(shí)間序列方法成功監(jiān)測(cè)到四川雅安地區(qū)漢源滑坡的災(zāi)前運(yùn)動(dòng)特征，各獨(dú)立軌道的雷達(dá)視線方向平均形變速率結(jié)果如圖3所示。其中，正值表示滑坡靠近衛(wèi)星運(yùn)動(dòng)，負(fù)值表示滑坡遠(yuǎn)離衛(wèi)星運(yùn)動(dòng)。升軌P128和P26軌道在2018～2020年變形分布相似，但變形幅度略有不同(圖3)。這主要是因?yàn)榛挛挥赑128軌道的近距離范圍，而位于P26軌道的遠(yuǎn)距離范圍，相對(duì)于同一目標(biāo)而言，其雷達(dá)視線方向存在9.9°的偏差。兩個(gè)升軌的平均形變速率結(jié)果顯示，滑坡存在明顯的滑前形變信號(hào)。在滑坡發(fā)生前的近兩年時(shí)間，P128軌道的最大平均形變速率可達(dá)到-42 mm·年，P26軌道達(dá)到-42.1 mm·年，非形變區(qū)大多數(shù)監(jiān)測(cè)點(diǎn)的平均形變速率位于-10～10 mm·年，兩個(gè)升軌的形變值為負(fù)值表明滑坡遠(yuǎn)離衛(wèi)星運(yùn)動(dòng)。從兩個(gè)軌道的平均形變速率結(jié)果可以看出，位移變化較大處主要發(fā)生在河流沿岸，這與現(xiàn)場(chǎng)識(shí)別的活動(dòng)滑坡位置一致。

為了更好地分析此次滑坡的滑前時(shí)空運(yùn)動(dòng)特征，對(duì)2018～2020年兩個(gè)升軌的形變時(shí)間序列進(jìn)行研究。圖4為P128和P26兩個(gè)軌道在雷達(dá)視線方向上的時(shí)間序列變形。從圖4可以看出，不同軌道滑坡運(yùn)動(dòng)的時(shí)空變化規(guī)律基本一致，滑坡在整個(gè)監(jiān)測(cè)期間存在持續(xù)變形。總的來說，滑坡在2018年底至2019年4月移動(dòng)迅速，而在之后移動(dòng)緩慢，基本處于穩(wěn)定狀態(tài)?；逻\(yùn)動(dòng)在空間上并不均勻，隨著時(shí)間的積累，滑坡形變區(qū)域逐漸形成形變中心，滑坡中心的累積形變可達(dá)到97 mm，而位于周邊區(qū)域的累積形變相對(duì)較小，最大累積形變僅為70 mm。

圖3 雷達(dá)視線方向平均形變速率結(jié)果Fig.3 Images of Average LOS Deformation Rate

兩個(gè)軌道的時(shí)間序列形變結(jié)果表明，漢源滑坡的災(zāi)前形變過程經(jīng)歷了形變快速積累和穩(wěn)定兩個(gè)階段。由于缺少可用的夏季干涉對(duì)，時(shí)間序列變形分析未獲取2019年7月至9月及臨滑前的時(shí)間序列信息，災(zāi)前形變分析并沒有捕捉到臨滑前的形變加速信號(hào)。

20181221是指2018年12月21日，其他依此類推圖4 雷達(dá)視線方向累積形變結(jié)果Fig.4 Images of LOS Cumulation Deformation

圖5 雷達(dá)視線方向單點(diǎn)累積形變統(tǒng)計(jì)結(jié)果Fig.5 Statistical Results of LOS Cumulation Deformation for Single Point

本文提取了漢源滑坡形變中心區(qū)域P點(diǎn)[圖3(a)]的雷達(dá)視線方向累積形變時(shí)間序列，結(jié)果如圖5所示。單點(diǎn)形變時(shí)間序列結(jié)果顯示，兩個(gè)升軌表現(xiàn)出明顯的滑前形變信號(hào)?；逻\(yùn)動(dòng)經(jīng)歷了兩個(gè)時(shí)段的形變：Ⅰ階段為2018年10月至2019年4月，形變加速積累，其中P128軌道累積形變量達(dá)97 mm，P26軌道累積形變量達(dá)86 mm；Ⅱ階段為2019年4月以后，滑坡運(yùn)動(dòng)速率幾乎為0，并在2020年8月21日發(fā)生滑坡。由于失相干現(xiàn)象，缺少夏季的干涉對(duì)，在2019年7月至2019年9月和臨滑前(即2020年6月至8月)并沒有獲得有效形變信號(hào)。2019年4月以后滑坡運(yùn)動(dòng)緩慢的具體原因還需結(jié)合后續(xù)工程地質(zhì)調(diào)查展開。

3.2 坡向平均形變提取與結(jié)果驗(yàn)證

為了分析漢源滑坡滑前沿坡向的運(yùn)動(dòng)特征，首先通過研究區(qū)DEM數(shù)據(jù)提取滑坡的坡度和坡向。雖然不同軌道數(shù)據(jù)集測(cè)得的雷達(dá)視線方向的形變量級(jí)不同，但是平行坡向的分量是相同的。因此，在提取坡度和坡向之后，結(jié)合式(3)分別計(jì)算了P128軌道和P26軌道的放大比例因子。為了控制坡向的偏差，剔除了絕對(duì)放大值大于一定閾值(本文取值為2)的孤立點(diǎn)，便于更好地分析滑坡沿坡向的運(yùn)動(dòng)規(guī)律。

圖6展示了P128軌道和P26軌道經(jīng)放大比例因子逐像素放大校正后得到的平行坡向的平均形變速率結(jié)果。從圖6可以看出，滑坡發(fā)生前在斜坡方向上存在明顯的災(zāi)前形變信號(hào)，滑坡方向的最大平均形變速率可以達(dá)到43.9 mm·年，顯著大于雷達(dá)視線方向?；麦w中心為滑前形變最劇烈的區(qū)域，整體呈現(xiàn)出明顯的漏斗狀特征。滑坡沿坡向的形變結(jié)果顯示，形變中心區(qū)域長約250 m，寬約250 m，形變中心海拔高度1 017 m，該區(qū)域是在滑坡發(fā)生前發(fā)生顯著形變的范圍。

由于沒有獲取研究區(qū)地面監(jiān)測(cè)數(shù)據(jù)，本文使用相同時(shí)間跨度的獨(dú)立SAR數(shù)據(jù)集來驗(yàn)證InSAR監(jiān)測(cè)結(jié)果的精度。通過統(tǒng)計(jì)2個(gè)軌道獨(dú)立觀測(cè)的SAR數(shù)據(jù)結(jié)果重疊點(diǎn)的差值，驗(yàn)證InSAR結(jié)果之間的一致性。統(tǒng)計(jì)結(jié)果顯示，2組相鄰軌道滑坡區(qū)域重疊點(diǎn)的坡向速率差值符合正態(tài)分布，速率差值的標(biāo)準(zhǔn)差為2.9 mm·年，表明InSAR監(jiān)測(cè)結(jié)果是可靠的。

綜合此次滑坡發(fā)生前InSAR時(shí)間序列結(jié)果，漢源滑坡存在明顯的災(zāi)前形變積累。在滑坡發(fā)生地為高山河谷的地形和砂泥巖地質(zhì)結(jié)構(gòu)背景下，經(jīng)過強(qiáng)降雨的加載作用，可能打破邊坡的穩(wěn)定階段。足夠的降水浸透滑坡體基底面，降低有效抗剪強(qiáng)度；當(dāng)邊坡剪切應(yīng)力超過剪切強(qiáng)度時(shí)，就可能觸發(fā)滑坡。結(jié)合滑坡發(fā)生前該區(qū)域經(jīng)歷的強(qiáng)降雨過程，推測(cè)在災(zāi)前形變積累情況下，強(qiáng)降雨可能是此次滑坡發(fā)生的誘發(fā)因素。

4 結(jié) 語

本文利用多軌道Sentinel-1衛(wèi)星影像，采用基于相干性的InSAR時(shí)間序列方法，對(duì)四川雅安地區(qū)漢源滑坡的災(zāi)前形變情況(2018年10月至2020年8月)進(jìn)行追溯。首先，獲取了滑坡沿雷達(dá)視線方向的平均速率與時(shí)間序列，結(jié)果顯示滑坡雷達(dá)視線方向的最大平均形變速率為42.1 mm·年；漢源滑坡從2018年底便發(fā)生了變形，至滑坡發(fā)生前最大累積形變量達(dá)97 mm，共經(jīng)歷了形變加速和穩(wěn)定兩個(gè)階段，形變?cè)诳臻g分布不均勻。然后，提取了滑坡沿坡向的形變，結(jié)果顯示在坡向上同樣存在明顯形變信號(hào)，最大坡向平均形變速率達(dá)43.9 mm·年。本文研究結(jié)果為此次漢源滑坡提供了第一手、相對(duì)完整的形變測(cè)量資料，所揭示的滑坡時(shí)空運(yùn)動(dòng)特征對(duì)于評(píng)估滑坡危害及進(jìn)一步實(shí)地調(diào)查至關(guān)重要。上述基于InSAR技術(shù)的災(zāi)前形變追溯方法也可以用來識(shí)別其他滑坡可能發(fā)生點(diǎn)的前兆運(yùn)動(dòng)。

圖6 漢源滑坡平行坡向的平均形變速率結(jié)果Fig.6 Images of Average Deformation Rate of Hanyuan Landslide Parallel to Aspect

值得注意的是，由于研究使用的Sentinel-1 SAR影像數(shù)據(jù)植被穿透能力較弱，在研究區(qū)植被茂密的情況下相干性會(huì)受影響。本文采用的基于相干性的InSAR時(shí)間序列方法雖然在一定程度上可以提高形變結(jié)果的點(diǎn)位密度，但由于夏季植被非常茂密，臨滑前干涉對(duì)相干性很差，所以本文形變追溯結(jié)果并未捕捉到臨滑前明顯的加速變形信號(hào)。這也說明InSAR技術(shù)在滑坡監(jiān)測(cè)預(yù)警中容易受到觀測(cè)環(huán)境影響，可能存在時(shí)間采樣率不足等缺點(diǎn)。

歐洲航天局提供了Sentinel-1衛(wèi)星影像，在此表示感謝!