應(yīng)用D-InSAR技術(shù)識別分析通渭縣李店鄉(xiāng)滑坡群災(zāi)害

王睿博，甘淑，張薦銘

(1.昆明理工大學(xué)國土資源工程學(xué)院，云南 昆明 650093；2.云南省高校高原山地空間信息測繪技術(shù)應(yīng)用工程研究中心，云南 昆明 650093)

1 引 言

滑坡是中國常見的地質(zhì)現(xiàn)象，每年由于滑坡災(zāi)害造成的損失不計其數(shù)。而滑坡群則是大范圍的滑坡，其災(zāi)害影響更勝于滑坡。2019年9月14日～9月16日，受連續(xù)降雨影響，甘肅定西市通渭縣常家河鎮(zhèn)李家店鄉(xiāng)發(fā)生特大型滑坡群災(zāi)害(圖1)。雖未造成人員傷亡，但是災(zāi)害造成的農(nóng)業(yè)、企業(yè)、供電、通訊、道路、橋梁等經(jīng)濟損失約 1 994.7萬元。

圖1 無人機航拍影像

合成孔徑雷達差分干涉測量(D-InSAR)是指利用同一區(qū)域的兩幅干涉圖像，其中一幅是形變前獲取的干涉圖像，另一幅是形變后獲取的干涉圖像，通過差分處理(除去地球表面，地形起伏)來獲取地表形變的測量技術(shù)。其作為一種先進的地面監(jiān)測技術(shù)，與常規(guī)方法相比，在識別大范圍地區(qū)的地表微小形變方面有著巨大優(yōu)勢，精度最高可達毫米，實現(xiàn)方法主要有3種：二軌法，三軌法以及四軌法。目前，已有許多相關(guān)學(xué)者將此技術(shù)應(yīng)用于形變觀測，文獻[1]利用D-InSAR技術(shù)監(jiān)測京津城際高速鐵路沉降情況，發(fā)現(xiàn)了沉降漏斗，并為鐵路方案的比選提供重要依據(jù)。文獻[2]以青藏鐵路拉薩段為試驗區(qū)，使用D-InSAR技術(shù)監(jiān)測了大范圍凍土形變，發(fā)現(xiàn)其形變情況非常吻合凍土凍脹融沉的物理變化規(guī)律。文獻[3]利用D-InSAR技術(shù)監(jiān)測了九寨溝震后高精度形變情況，快速高效地獲取震損物源的分布與數(shù)量特征。文獻[4]利用D-InSAR技術(shù)監(jiān)測了河北峰峰煤礦地區(qū)地表沉降情況，其利用多模式雷達，獲取了該礦區(qū)的沉降量，精確地反映了礦區(qū)的地理信息。文獻[5]利用D-InSAR技術(shù)監(jiān)測米林滑坡，獲取滑坡形變值，證明了D-InSAR技術(shù)用于滑坡早期識別與動態(tài)監(jiān)測的可行性。但D-InSAR仍存在許多局限，如受到影像個數(shù)、DEM精度、時空間基線等影響，均會造成形變結(jié)果精度不同。本實驗則利用二軌D-InSAR技術(shù)對此滑坡群形變進行觀測，并根據(jù)形變結(jié)果分析不同精度DEM下的滑坡范圍，判定災(zāi)害原因。

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

研究區(qū)(圖2)位于通渭縣李店鄉(xiāng)，其地形山高坡陡，土地貧瘠，自然條件嚴(yán)酷。李店鄉(xiāng)是一個干旱鄉(xiāng)，年降水量 380 mm左右，年蒸發(fā)量則在 1 500 mm以上。根據(jù)文獻[6]研究記載，通渭縣受到了1718年通渭7.8級大地震和1920年海原8.5級大地震的影響，在境內(nèi)誘發(fā)了大量的黃土地震滑坡，現(xiàn)有的地形由大面積的滑坡體組成，結(jié)構(gòu)較為松散。從土質(zhì)方面來看，該地區(qū)與文獻[7]研究區(qū)域地理位置相近，土質(zhì)為黃土，一般呈淺黃色，具大孔隙，垂直節(jié)理發(fā)育，濕陷性強，易崩塌，黃土溶蝕微地貌發(fā)育，存在不穩(wěn)定性[6]。黃土濕陷系數(shù)隨著天然含水率、飽和度、干密度、液限、塑限、塑性指數(shù)的增大而減小，而隨著孔隙比、壓縮系數(shù)的增加呈增大的趨勢[7]。

本次研究選取的數(shù)據(jù)來源于哨兵1A(Sentinel-1)衛(wèi)星，模式為干涉測量寬幅(IW)的SAR影像數(shù)據(jù)，并使用了與SAR影像相匹配的精確軌道數(shù)據(jù)(表1)。

同時，本實驗采用了不同分辨率精度的DEM，分別是 90 m，30 m以及 12.5 m(圖3)，用以比較不同精度下形變量范圍。根據(jù)文獻[8]結(jié)論，90 m與 30 m精度的DEM選用SRTM數(shù)據(jù)——主要是由美國航空航天局(NASA)和國防部國家測繪局(NIMA)聯(lián)合測量，由SRTM數(shù)據(jù)中心(http：//srtm.csi.cgiar.org/srtmdata/)獲取。12.5 m精度的DEM則是來自于日本的ALOS衛(wèi)星，其搭載的PALSAR傳感器所制作的DEM，精度約為 12.5 m，可滿足此次實驗需求，由歐空局網(wǎng)站(https：//search.asf.alaska.edu/#/)獲取。

Sentinel-1干涉對參數(shù)表 表1

圖2 研究區(qū)范圍

圖3 3種精度DEM

3 D-InSAR處理過程

本次D-InSAR技術(shù)應(yīng)用實驗主要采用二軌法，其基本思想是利用已知的外部DEM來去除地形引起的干涉相位，其關(guān)鍵步驟是外部DEM與主影像的配準(zhǔn)[9，10]。因此，本實驗利用SARscape軟件，運用不同精度的DEM來進行D-InSAR配準(zhǔn)，以提高實驗精度，判定滑坡區(qū)域。二軌D-InSAR處理流程如圖4所示。

圖4 D-InSAR處理流程圖

3.1 基線估算與干涉圖生成處理

(1)基線估算

基線估算主要用于評價干涉相對的質(zhì)量。此步驟主要計算空間基線、時間基線、多普勒偏移等系統(tǒng)參數(shù)，用以判斷獲取的SAR數(shù)據(jù)影像是否能夠做干涉，這是D-InSAR技術(shù)進行的基礎(chǔ)?？臻g基線是兩個傳感器之間相對位置的向量，可以將其分解為兩個基線分量：平行基線和垂直基線。其中，垂直基線是判定影像能否做干涉的最重要的一個參數(shù)，若其超過了臨界基線最小或最大值時，兩幅圖像則無相位信息，沒有相干性，即無法做干涉。垂直基線的長短也反映了干涉測量系統(tǒng)對高程的敏感程度[9]。

本次實驗中，將兩景SAR影像導(dǎo)入SARscape軟件后，對VV極化影像裁剪出研究區(qū)域，再對裁剪后的兩景影像進行基線估算，得到空間基線長度為 102.356 m，臨界基線的最小/最大值為 -5 452.484 m與 5 452.484 m，2π模糊位移為 0.028 m。由此可見，空間基線遠小于臨界基線，存在相位信息，有相干性，可以滿足D-InSAR技術(shù)的處理要求。

(2)干涉圖生成處理

裁剪后的兩幅SLC影像共軛相乘，并分別與3幅 90 m，30 m和 12.5 m精度的DEM進行配準(zhǔn)，此過程設(shè)置相同的制圖分辨率，配準(zhǔn)結(jié)束后分別生成了3幅精度不同的干涉條紋圖。此時還無法用目視判定不同精度DEM配準(zhǔn)的影像差別，但經(jīng)過去除平地效應(yīng)之后，根據(jù)去平干涉圖像(圖5)，已經(jīng)可以大致看出，DEM精度不同，去平圖像已經(jīng)有很大差別，但此時去平后的干涉條紋圖的噪聲還異常明顯，必須要進行去噪處理。

圖5 部分去平干涉圖

3.2 濾波相干處理與相位解纏

(1)濾波相干處理

SARscape提供了3種濾波器，分別是Adaptive濾波器、Boxcar濾波器和Goldstein濾波器。本實驗均采用Goldstein濾波器，實驗處理得到濾波后的干涉條紋圖(圖6)以及干涉圖的相干系數(shù)圖(圖7)。相干系數(shù)圖中，相干系數(shù)越大，表示相干性越高，反之，則表示變化較大，相干性較低。

圖6 部分濾波干涉圖

圖7 相干系數(shù)圖

(2)相位解纏

相位的變化是以2π為周期的，當(dāng)相位的變化量超過2π時，相位就會重新開始，并不斷開始循環(huán)。因此，相位解纏是對前幾個步驟(去平地效應(yīng)和濾波)處理后的相位進行解纏處理，即把纏繞相位解纏繞后，恢復(fù)成真實的干涉相位，目的主要是解決2π的模糊問題，使地形上的線性變化的信息與相位一一對應(yīng)。此步驟也是InSAR數(shù)據(jù)處理過程中最重要和最困難的步驟。

SARscape提供了3種解纏方法：Region Growing(區(qū)域增長法)、Minimum Cost Flow(最小費用流法)和Delaunay MCF。其中，Minimum Cost Flow是默認的解纏方法，這種方法是將相位解纏問題轉(zhuǎn)化為最小化問題，通過在全局范圍內(nèi)搜索路徑和最短枝切來求得最小化問題的最優(yōu)解。該方法可以應(yīng)用于規(guī)則網(wǎng)絡(luò)(網(wǎng)格)，也可以用于不規(guī)則網(wǎng)絡(luò)(三角網(wǎng))[11，12]。用此方法解纏時，優(yōu)先采用正方形的格網(wǎng)，考慮所有像元，當(dāng)像元的相干性小于所設(shè)閾值時，可對該像元做掩膜處理。相較而言，當(dāng)解纏相對困難時，Minimum Cost Flow比Region Growing可以取得更好的結(jié)果。

通過實驗對比，此處解纏選擇默認的Minimum Cost Flow方法效果較好。由此得到相位解纏后條紋圖(圖8)。

圖8 部分相位解纏后條紋圖

3.3 相位精確修正與形變圖示表達

(1)軌道精煉與重去平優(yōu)化

軌道精煉和重去平是利用選取的GCP控制點，將干涉后的相位和解纏后的相位精確修正。在選取GCP控制點時，其原則應(yīng)該遠離形變區(qū)域、解纏錯誤的相位躍變區(qū)域和殘余地形相位。因此，平地區(qū)域是比較好的選擇控制點的位置。將GCP控制點作為穩(wěn)定的參考點，通過計算可以得到控制點的相位誤差，從而將這些誤差去除，并重新修改解纏圖像頭文件中的軌道參數(shù)，使得相位更加精確。

軌道精煉主要有軌道優(yōu)化和多項式優(yōu)化兩種方法，如果選擇軌道優(yōu)化的方法，則至少需要7個控制點，否則會自動改為多項式優(yōu)化。多項式優(yōu)化原理是從解纏相位中估算相位斜坡，但要求控制點個數(shù)必須與多項式次數(shù)所需要的控制點個數(shù)相對應(yīng)，否則多項式次數(shù)會自動降低[13]。多項式優(yōu)化雖然條件較為苛刻，但其健壯性更好，即使是在基線小的情況下也可以使用。

本實驗選取Automatic Refinement(自動優(yōu)化法)，讓SARscape自行判定最佳優(yōu)化方式。

(2)相位轉(zhuǎn)形變和地理編碼

此步驟把合成的相位與經(jīng)過絕對校準(zhǔn)的解纏相位相結(jié)合，轉(zhuǎn)換為形變數(shù)據(jù)，再通過地理編碼將形變量顯示到制圖坐標(biāo)系統(tǒng)中，默認得到的是LOS方向的形變，最終得出形變結(jié)果。如圖9所示，圖中形變量的值可能是由于滑坡形變方向和數(shù)據(jù)的LOS方向相差較大，從而導(dǎo)致測得數(shù)據(jù)值與實際形變值相比總體偏小。

圖9 形變區(qū)域圖

4 滑坡識別應(yīng)用分析探討

4.1 不同DEM精度的形變監(jiān)測滑坡識別

將圖9中3種不同精度DEM形變結(jié)果分別通過ArcGIS軟件將形變區(qū)域放大分析，并將3種結(jié)果疊加，提取出主要沉降與抬升區(qū)域(圖10(g))，用以輔助分析滑坡。經(jīng)分析，分別得到三種DEM精度研究區(qū)形變圖(圖10(a)、(b)、(c))。

將ArcGIS中3種形變結(jié)果處理為圖層文件，轉(zhuǎn)換KML格式，疊加至Google Earth經(jīng)緯度相同區(qū)域，形成Google Earth疊加形變圖(圖10(d)、(e)、(f))，此圖是將平面二維形變信息轉(zhuǎn)化為立體三維信息。

圖10 分析圖

現(xiàn)結(jié)合圖10中3種不同精度DEM形變結(jié)果及Google Earth地形情況對滑坡區(qū)域進行詳細分析：

(1)根據(jù)90 m精度DEM形變結(jié)果，結(jié)合Google Earth疊加地形效果，可以看到滑坡群明顯滑坡有5處，形變范圍較為明顯，形變值較其他兩種精度DEM所做結(jié)果偏大。A滑坡處坡度較緩，有明顯的過渡區(qū)域，形變量依次由大到小，山頂滑坡形變最大，約為 0.041 m～0.047 m，到山谷有明顯抬升部分，約為 0.008 m～0.014 m。B滑坡處山體較陡，過渡區(qū)域較少，滑坡形變較為明顯。C、D、E三處滑坡與B滑坡情況較為相似，都屬于山體較陡區(qū)域，形變量變化速度較快。

(2)根據(jù)30 m精度DEM形變結(jié)果，測出6處滑坡，A、B、D、E、F區(qū)域與 90 m精度DEM形變結(jié)果在形變范圍方面基本相似，形變量范圍介于 90 m精度DEM與 12.5 m精度DEM測得結(jié)果之間，故此精度形變值可靠性較高，形變量最大處的值約為 0.035 m～0.041 m，抬升處形變值約為 0.007 m～0.012 m。新測得C處滑坡，屬于小范圍滑坡，其形變量變化范圍不大。

(3)12.5 m精度DEM形變結(jié)果總體相似于 30 m精度DEM形變結(jié)果，但是其形變區(qū)域較為零散，形變值相較于其他兩種精度DEM測得結(jié)果整體偏小，抬升區(qū)域不明顯，總體判斷滑坡較為困難。

4.2 滑坡識別結(jié)果分析探討

通過分析3種不同精度的DEM測出的形變結(jié)果，可以大致判斷滑坡群主體范圍：李店鄉(xiāng)東部及東北部方向，以陽山里村西部區(qū)域形變量最大，滑坡區(qū)域最密集。以 30 m精度DEM顯示效果最好，既能顯示出所有的形變區(qū)域，又能顯示出較為突出的形變量，整體相較于其他兩種精度DEM測得效果更好，此結(jié)果與文獻[14]結(jié)論相同。

現(xiàn)利用30 m精度DEM測得形變結(jié)果，對6處滑坡進行量化分析，分別在坡頂和坡腳各選取一點并連接成線，生成此線上6幅地形剖面形變圖，分布如圖11(a)所示。

從圖11(b)滑坡剖面圖可以看出，6處剖面圖均符合滑坡特征，坡頂?shù)狡碌卓傮w呈現(xiàn)下沉趨勢，中部有滑坡物堆積部分，形變量有明顯抬升。由形變量及像元個數(shù)可以判定出A滑坡的波浪形堆積物大致是從第5像元到第30像元處，約有 0.001 m～0.004 m的抬升，總抬升量約有 0.01 m；B處滑坡堆積物存在于第7像元到第12像元區(qū)域，約有 0.005 m抬升；C滑坡堆積物存在于6～10像元，約抬升了 0.013 m；D滑坡相似于C滑坡，堆積物存在于4～14像元，約有 0.012 m抬升；E滑坡抬升區(qū)域較少，大致是由于滑坡體未完全滑下，屬于整體性糯滑，分布于第5～9像元區(qū)域，抬升了 0.008 m；F滑坡堆積物最多，根據(jù)圖10(b)(e)及圖11(a)可以看出，此處為較大U型山谷，四周滑坡均往中部堆積，故剖面顯示出較高的抬升量，約有 0.021 m，存在于約第9～26像元，跨度約有17個像元。

圖11 剖面分析圖

結(jié)合圖10中(d)(e)(f)，滑坡最大形變處(深藍色區(qū)域)發(fā)生在道路左側(cè)，說明滑坡對道路影響較大，結(jié)合無人機航拍，符合現(xiàn)場情況。Google Earth疊加影像可以明顯看出山體區(qū)域兩側(cè)均有整體性的滑動形變，根據(jù)滑坡的特點，滑坡變形范圍呈現(xiàn)的是連續(xù)的面狀體，通常形狀上陡下緩，近似于圓弧形，由圖11中可以看出形變部分也呈現(xiàn)面狀特點，且滑坡變形部分(藍色區(qū)域)均分布在坡體較陡區(qū)域，其周圍區(qū)域范圍內(nèi)也出現(xiàn)了零散的抬升區(qū)域(黃色區(qū)域)，在經(jīng)過三次不同精度DEM的處理后，每次處理結(jié)果均有形變，形變量及形變主體范圍基本不變。這說明了此處形變在這12天里空間和時間上是連續(xù)變化的，因此實驗結(jié)果符合滑坡特征。

由于滑坡還具有“大雨大滑、小雨小滑、無雨不滑”的特點，滑坡當(dāng)日，結(jié)合天氣情況，通渭縣李店鄉(xiāng)地區(qū)降水量達 72.6 mm，已經(jīng)遠超平均日降水量，隨著降水入滲土地中，致使斜坡土體含水量增大，斜坡體自重也隨之增大，降雨入滲到達滑坡帶，使滑坡帶飽和，使得土體抗剪強度顯著降低，從而引發(fā)斜坡失去穩(wěn)定性并蠕滑變形，最終發(fā)生滑坡群災(zāi)害。

5 結(jié) 論

本此實驗以通渭縣常家河鎮(zhèn)李店鄉(xiāng)滑坡群為觀測目標(biāo)，利用二軌D-InSAR技術(shù)對滑坡群進行差分干涉測量，獲取了不同精度DEM下的形變值、形變范圍。通過比較，最終以 30 m精度DEM測出的形變效果最好，以此判定了6處滑坡群主體部分，均發(fā)生在李店鄉(xiāng)東部及東北部方向，以陽山里村西部區(qū)域形變量最大，滑坡區(qū)域最密集。根據(jù)形變值與形變范圍，結(jié)合滑坡特點、剖面情況、無人機航拍、土質(zhì)狀況及天氣情況，判定滑坡原因是降水量過大引發(fā)的災(zāi)害。

此次實驗只利用D-InSAR技術(shù)觀測了研究區(qū)12天內(nèi)滑坡群的沉降值，并未研究其時序變化。因此，若要持續(xù)監(jiān)測此區(qū)域，獲取時間序列上的變化情況，還可以利用PS-InSAR、SBAS-InSAR等技術(shù)進一步研究。