王睿博，甘淑，張薦銘

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

1 引 言

滑坡是中國(guó)常見(jiàn)的地質(zhì)現(xiàn)象，每年由于滑坡災(zāi)害造成的損失不計(jì)其數(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)濟(jì)損失約 1 994.7萬(wàn)元。

圖1 無(wú)人機(jī)航拍影像

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

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

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

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

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

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

圖2 研究區(qū)范圍

圖3 3種精度DEM

3 D-InSAR處理過(guò)程

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

圖4 D-InSAR處理流程圖

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

(1)基線估算

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

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

(2)干涉圖生成處理

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

圖5 部分去平干涉圖

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

(1)濾波相干處理

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

圖6 部分濾波干涉圖

圖7 相干系數(shù)圖

(2)相位解纏

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

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

通過(guò)實(shí)驗(yàn)對(duì)比，此處解纏選擇默認(rèn)的Minimum Cost Flow方法效果較好。由此得到相位解纏后條紋圖(圖8)。

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

3.3 相位精確修正與形變圖示表達(dá)

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

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

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

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

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

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

圖9 形變區(qū)域圖

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

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

將圖9中3種不同精度DEM形變結(jié)果分別通過(guò)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地形情況對(duì)滑坡區(qū)域進(jìn)行詳細(xì)分析：

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

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

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

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

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

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

從圖11(b)滑坡剖面圖可以看出，6處剖面圖均符合滑坡特征，坡頂?shù)狡碌卓傮w呈現(xiàn)下沉趨勢(shì)，中部有滑坡物堆積部分，形變量有明顯抬升。由形變量及像元個(gè)數(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個(gè)像元。

圖11 剖面分析圖

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

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

5 結(jié) 論

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

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