廖 軍，吳彩燕，王立娟，譚秋焰，朱新婷

（1.西南科技大學環(huán)境與資源學院，四川綿陽 621010；2.四川省安全科學技術研究院，四川成都 610045）

引言

滑坡是指巖體、碎屑、土壤沿斜坡自上而下滑動的現(xiàn)象，也指斜坡在重力作用下沿著坡體內(nèi)部一個或者多個面做剪切運動的現(xiàn)象［1～3］，在我國的山區(qū)和峽谷地帶廣泛分布，尤其是西南高山區(qū)和三峽庫區(qū)，發(fā)生頻率最高，造成的損失最嚴重［4～5］。常規(guī)的滑坡識別手段主要有野外實地勘查和遙感影像識別，野外實地勘查能準確圈定滑坡要素，但其工作既費時又低效，并且勘測區(qū)域受到地形的極大限制［6］；遙感影像識別滑坡大多是利用專家經(jīng)驗知識通過光學遙感影像的光譜、空間、紋理等信息進行的［7～8］，盡管已經(jīng)取得了相當不錯的成就，但光學影像容易受到天氣因素的限制如大霧、降雨等，并且在監(jiān)測和識別滑動不明顯的潛在滑坡能力不足［9～10］。而新出現(xiàn)的合成孔徑雷達干涉測量（InSAR）技術不僅能精確獲取地表微小形變信息，且具有全天候、全天時、時效性高等優(yōu)點，成為了國內(nèi)外學者研究的熱點［11～14］。

近年來，InSAR技術逐漸應用于滑坡災害識別，并取得了一些成功的案例。廖明生等以三峽庫區(qū)滑坡為對象，使用高分辨SAR數(shù)據(jù)對其進行監(jiān)測，成功獲得了滑坡發(fā)生的時間、位置及形變情況，證明了高分辨率雷達數(shù)據(jù)有效應用在滑坡監(jiān)測中［15］；Tofani V等［16］采用永久散射體干涉測量（Persistent Scatterer，PS-InSAR）技術對意大利北部的Santo Stefano d′Aveto滑坡進行分類與監(jiān)測，結果表明PS技術結合地面數(shù)據(jù)能夠準確判別滑坡的邊界和運動狀態(tài)；Bianchini等［17］以意大利的Calabria地區(qū)為研究區(qū)，利用PS-InSAR技術對該地區(qū)的滑坡進行了識別及滑坡編目的更新；趙超英等［18］利用多源SAR數(shù)據(jù)來識別黑方臺地區(qū)的潛在滑坡，結合實地調(diào)查和光學影像驗證了識別結果的可靠性，并分析典型滑坡的失穩(wěn)模式；Dong等［19］使用相干散射體（Coherent Scatterer，CS-InSAR）技術更新了丹巴縣甲居滑坡的活動邊界，繪制了更詳細的滑坡位移圖?，F(xiàn)有的研究雖然在一定程度上可以對潛在滑坡進行識別，但是這些方法大多是從雷達視線方向（LOS）角度出發(fā)進行識別研究，雖然可以識別滑坡發(fā)生的時間和地理位置，但LOS方向與地表和山坡之間存在一定的夾角，不能直接反映滑坡的變形程度，也不能得到滑坡變形的真實時間序列，難以識別滑坡前后的具體變形趨勢。

本研究以滑坡災害多發(fā)的水城縣為研究區(qū)，通過改進的SBAS-InSAR技術，獲得水城縣坡度方向形變信息，結合水城縣雷達可視分區(qū)及谷歌影像剔除無效觀測值，利用Anselin Local Moran′s I對獲取的監(jiān)測值進行異常值分析和聚類處理，獲取低值聚集區(qū)域，最后結合谷歌影像的地表特征，實現(xiàn)水城縣潛在滑坡的識別，用于指導該區(qū)域滑坡災害的預測預警。

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

1.1 研究區(qū)概況

水城縣位于貴州省六盤水市，鄰接云南省，在水系和構造運動作用下，西北地勢較東南高。該區(qū)屬于亞熱帶濕潤季風氣候地區(qū)，降雨集中在5～10月份，其中6～8月降雨量占全年53.03%［20］，縣內(nèi)河網(wǎng)密閉，水量集中，是造成境內(nèi)地質(zhì)災害的主要誘因之一。水城縣地質(zhì)構造復雜，存在多個構造變形區(qū)，地層巖性以二疊系、三疊系以及第四系為主。

水城縣境內(nèi)滑坡廣泛分布在飛仙關組、峨眉山組及龍?zhí)督M等地層中，該地層中節(jié)理發(fā)育，巖層破碎，極易在人類活動或降水等因素下發(fā)生滑坡等地質(zhì)災害［20～21］。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

1.2 數(shù)據(jù)源介紹

文中采用的數(shù)據(jù)源包括歐空局新一代C波段的中高分辨率衛(wèi)星Sentinel－1A SAR數(shù)據(jù)，選取了2018年1月～2018年12月共20景降軌數(shù)據(jù)，入射角為39.5°，訪問周期為12 d，拍攝模式為IW，空間分辨率為5×20m（Rg×Az），極化方式為VV，并使用POD精密定軌星歷數(shù)據(jù)（Precise Orbit Ephemerides，POD）進行軌道矯正，定位精度優(yōu)于5 cm；數(shù)字高程模型（DEM）為30 m分辨率的SRTM DEM（Shuttle Radar Topography Mission，SRTM）；從Google Earth中獲取的研究區(qū)光學遙感影像數(shù)據(jù)；野外調(diào)查資料。

2 研究方法

本研究選取2018/1～2018/12共20景Sentinel－1A SAR降軌數(shù)據(jù)，利用永久散射體干涉測量（PS-InSAR）技術提取獲得相干性高且相對穩(wěn)定的散射體用做地面控制點（Ground Control Point，GCP），引入到小基線集干涉測量（SBAS-In-SAR）技術處理流程中，獲取長時間序列的地表形變信息（LOS向），并轉換為坡度方向形變信息；結合水城縣地形可視性分區(qū)及谷歌影像，剔除不可視區(qū)（陰影、透視收縮、疊掩）、高密度林地及灌木林、河流等不可靠的形變信息，保留其他區(qū)域可靠形變信息；對獲取的可靠形變信息采用Anselin Local Moran′s I指數(shù)進行空間異常值分析和聚類處理，獲取低值聚集區(qū)域，并結合谷歌影像的地表特征，從而對潛在滑坡進行識別?；伦R別路線如圖2。

圖2 潛在滑坡識別流程Fig.2 Landslide identification process

2.1 結合PS-InSAR技術的GCP點選擇

在SBAS-InSAR處理流程中，選取穩(wěn)定的地面控制點（GCP）進行軌道精煉和相位偏移的計算，可以有效地估算和去除殘余相位及可能存在的斜坡相位，進而消除平地效應，對相位轉形變的可靠性起著重要作用［16－18］。

由于人工選點會對形變速率引入誤差，文中利用PS-InSAR技術的優(yōu)勢［22］，依次設置振幅離差閾值（≤0.25）、相干性閾值（≥0.8）以及空間域的低通濾波和時域的高通濾波，減小干涉相位中的大氣延遲相位誤差，得到可靠的穩(wěn)定散射點的形變信息，并對該結果設置速率閾值選取在－0.1～0.1 mm的穩(wěn)定點作為GCP點，通過地理坐標到SAR斜距坐標系的轉換，即可應用于在SBAS處理流程中。

振幅離差指數(shù)（D）是依據(jù)SAR圖像上同一位置像元的強度值在時間序列上的離散特性進行目標的識別，其計算公式如下［23］：

式中：μ是時序上的某個點的振幅平均值；σ是時序上的振幅的標準差。

2.2 SBAS-InSAR技術原理及方法

小基線集干涉測量（SBAS-InSAR）是在差分干涉測量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）和永久散射體干涉測量（PS-InSAR）發(fā)展起來的一種時序監(jiān)測方法，該方法能有效消弱常規(guī)D-InSAR中大氣效應，時空失相干等因素對監(jiān)測結果的影響，同時彌補了PS-InSAR技術不適用于相干點稀疏地區(qū)的不足［10，16～17］。在2002年，Berardino等［24］首次提出了小基線集時序監(jiān)測方法，SBAS-InSAR能夠獲取非線性的形變結果，在滑坡形變監(jiān)測中具有良好的適應性。

小基線集的原理是通過設置時間、空間基線的閾值，將多期SAR影像進行自由配對組合生成多個子集，每個子集由超主圖像連接，從而形成一個大集合，保證在每個小集合中SAR影像間的基線小，集合間的SAR影像基線大，從而提高數(shù)據(jù)的采樣率，克服空間基線過長導致的空間失相干影響。配合外部高精度的DEM消除地形相位，生成一系列差分干涉圖，對差分干涉圖進行篩選得到高相干性像元，然后利用地面控制點進行軌道精煉，剔除低相干性和解纏錯誤的干涉像對，采用奇異值分解法和最小二乘法對多個小基線集聯(lián)合求解，從而得到研究區(qū)的時間形變序列及地表平均形變速率［25－26］。

大多數(shù)滑坡都是沿坡面滑動的，而小基線集干涉測量法只能監(jiān)測雷達視線的形變，不足以反映滑坡的真實變形。因此，為了準確地描述滑坡真實的形變情況，將LOS向形變速率轉換為沿坡度方向形變速率［27－28］：

式中：V s l o p e為沿坡度方向的形變速率；V los為LOS向的形變速率；β為雷達視線方向與斜坡坡面的夾角，即視坡夾角；φ為斜坡的坡度/（°）；α為斜坡的坡向/（°）；θ為雷達入射角/（°）；αs為衛(wèi)星軌道方向與正北方向的夾角，即衛(wèi)星飛行的方向，升軌數(shù)據(jù)為負，降軌數(shù)據(jù)為正。

2.3 聚類和異常值分析原理

采用Anselin Local Moran′s I指數(shù)對SBAS-InSAR技術獲取的形變干涉點進行空間域的異常值分析和聚類分析［29～30］，通過計算要素與其他所有要素作為相鄰要素的數(shù)值關系，從而獲得數(shù)據(jù)集中有統(tǒng)計顯著性的高值（熱點）和低值（冷點），自動對形變點中的異常數(shù)據(jù)進行聚合，確定適當?shù)目臻g分析范圍，并校正數(shù)據(jù)的空間相關性，得到形變聚類點置信水平p值和得分，其相關的計算公式如下：

統(tǒng)計量z I i得分的計算方法如下：

式中：I i為形變點i的Anselin Local Moran′s I指數(shù)值；x i、x j分別為點i，j的屬性；n為選定的窗口內(nèi)要素的總數(shù)；是對應屬性的平均值；w i,j是要素i和j之間的空間權重；S2和S分別為為選定的窗口內(nèi)像元的方差、標準差；μ0為總體均值；Z I i為形變聚類點i的得分；p i為形變聚類點i的顯著性置信水平。

3 數(shù)據(jù)處理

3.1 雷達可視性分析

雷達遙感通過向地面發(fā)送脈沖信號，接收地面反射的回波信號，實現(xiàn)對地物的觀測。盡管具有全天時、全天候、大面積等無法比擬的優(yōu)勢，由于SAR采用側視成像方式，其入射角（θ）與實際地形存在夾角，使得SAR影像會出現(xiàn)陰影、透視收縮、疊掩3類常見的幾何畸變［27］。根據(jù)Sentinel－1A衛(wèi)星降軌條件下的幾何參數(shù)及水城縣DEM提取的地形參數(shù)（坡度、坡向）可獲得水城縣可視性分區(qū)方式，并獲得水城縣可視性區(qū)域（見表2）。

表2 SAR影像可見性分區(qū)Table 2 SAR image visibility zoning

Sentinel－1A降軌影像在水城縣的可視性區(qū)域分布如圖3所示，可見區(qū)域約占35.76%，低敏感區(qū)域約占24.08%，不可見區(qū)域約占40.16%。其中，可見區(qū)域以主要分布于東南、西南區(qū)域；不可見區(qū)域主要分布在中心區(qū)域。

圖3 研究區(qū)雷達可視性分區(qū)Fig.3 Radar visibility zoning of in study area

3.2 基于改進SBAS技術的數(shù)據(jù)處理

文中使用精密軌道數(shù)據(jù)（POD）和30m分辨率的SRTM DEM進行影像軌道誤差和地形相位校正。為了提高時間、空間相干性，保證干涉圖質(zhì)量，將時間基線設置為120 d，空間基線閾值為極限基線的5%（5 162.12×0.05=258.106 m），使用3D解纏。超級主影像為2018年7月3日，共生成了138個干涉像對，如圖4所示。通過查看生成的相干系數(shù)圖、去平地干涉圖、相位解纏圖，剔除相干性較差的干涉對，最終使用了115個干涉相對進行后續(xù)處理。

圖4 數(shù)據(jù)集時空基線分布Fig.4 Spatiotemporal baseline distribution of the dataset

圖5（a）是去平地后的相位干涉圖，水城縣由于植被覆蓋密集，地形起伏較大，導致地形相位去除效果一般，為了減小植被等因素的影響，選取了Delaunay MCF方法進行相位解纏；為了獲得較高的地面分辨率，將方位向與距離向的視數(shù)比設置為1：4；為了獲得更為清晰干涉條紋，選取了Goldstein濾波進行噪聲去除，濾波后的干涉圖結果如圖5（b），可看出濾波后噪聲誤差減小，干涉條紋的清晰度提高。

圖5 研究區(qū)部分干涉處理結果圖Fig.5 Results of partial interference processing in the study area

3.3 坡度方向可靠形變信息獲取

基于改進的SBAS-InSAR技術，獲取了2018年水城縣LOS方向的地表變形速率，對于西南山區(qū)，LOS向的形變速率不足以真實反映斜坡的真正形變情況，因此將LOS向形變速率轉換到坡度方向形變速率。當β接近90°時，cosβ近乎為0，導致V slop e趨于無窮大。為了對產(chǎn)生的異常值進行剔除，Herrera等以cosβ=±0.3為閾值，即V sl ope≤3.33V los。設定0＜C＜0.3時，C=0.3；-0.3＜C＜0時，C=-0.3［27－28，31－32］。V slope表示沿坡度方向的位移速率，塊體移動方向只能沿斜坡向下，故剔除V slope為非負的形變點。

基于水城縣可視性分區(qū)結果，先剔除了由透視收縮、疊掩和陰影不可信區(qū)域的形變信息，保留了好可視區(qū)域和低敏感區(qū)域中的形變信息。由于水城縣南部地區(qū)常年覆蓋著茂密的植被，盡管處于可見區(qū)域內(nèi)，但C波段地表穿透能力有限，即便SAR數(shù)據(jù)時間間隔短，也會導致失相干，得到有效監(jiān)測值較少且監(jiān)測精度有限。此外，滑坡在河流及平坦的地區(qū)發(fā)生的可能性較低。因此，在可視性分區(qū)結果的基礎上，結合谷歌影像，將密度較高的林地、灌木林、河流以及坡度小于5°不可靠區(qū)域的形變信息進行剔除，保留其它區(qū)域可靠形變信息，得到水城縣坡度方向年均形變速率（圖6）。

圖6 研究區(qū)2018年坡度方向年均形變速率Fig.6 Average annual deformation rate of the slope direction in study area in 2018

4 結果與分析

4.1 水城縣變形監(jiān)測結果分析

從圖6中可以看出整個水城縣北部的形變速率較為突出，各地形變速率不均，存在多個形變中心，集中分布在陡坡上。區(qū)域1形變較為嚴重，年平均速率在－75～0 mm/a，礦區(qū)開采處的變速率最明顯，形變速率達到了－142.3 mm/a；區(qū)域2整體地形起伏較大，形變程度也較為突出，其年平均形變速率在－35～0 mm/a之間；區(qū)域3受到雷達監(jiān)測制約以及植被覆蓋影響，獲取的監(jiān)測結果有限，形變速率大多在－10～0mm/a，少部分區(qū)域形變速率達到了－20 mm/a左右，整體較為穩(wěn)定。

由上述分析可知，研究區(qū)內(nèi)的變形區(qū)具有分布不均勻的特點，主要是人類活動因素（礦區(qū)開采、修筑道路）及差異化的地理環(huán)境因素（坡度、坡向、植被覆蓋率等）造成的。根據(jù)已有資料分析，區(qū)域1分布著較多的礦山，受礦區(qū)開采的影響，導致該區(qū)域整體形變量較大，易導致滑坡的發(fā)生。區(qū)域2地形起伏大、山體陡峭、植被覆蓋率低，受地形地貌等地理環(huán)境因素的影響，其山體邊坡的穩(wěn)定性會明顯低于區(qū)域3，發(fā)生滑坡的可能性相對較高。

對獲取的形變信息進行統(tǒng)計分析，形變速率在－15～0 mm/a的形變點占據(jù)總數(shù)的40.5%，形變速率在－35～－15 mm/a的形變點占據(jù)總數(shù)的31.5%，形變速率在－55～－35 mm/a的形變點占據(jù)總數(shù)的13.5%，形變速率大于－55 mm/a的形變點占總數(shù)14.5%。從統(tǒng)計分析結果來看，形變速率大部分集中在－35～0 mm/a，其中形變速率小于－30 mm/a所占比例最大，在一定程度上能夠表明水城縣地表整體較為穩(wěn)定，形變主要集中在水城縣的北部區(qū)域。。

根據(jù)上述對水城縣坡度方向形變速率分析，結合國內(nèi)外眾多學者針對滑坡事故做出的研究成果［13－18，22－23］，文中選取V slop e=15 mm/a為速率閾值，即V slope＞15 mm/a的區(qū)域易發(fā)生滑坡事故。

4.2 水城縣滑坡識別及驗證

4.2.1 滑坡識別

通過Anselin Local Moran′s I對坡度方向形變點進行空間域異常值分析和聚類處理，會得到HH、HL、LH、LL四種聚類類型，其中HH表示高值聚集，HL表示低值包含高值異常，LH表示高值包含低值異常，LL表示低值聚集。通過獲得低值聚集區(qū)域，結合滑坡在谷歌影像上的形狀、分布、微地貌特征，確定滑坡位置及邊界。通過本文的方法最終共識別出了14處潛在滑坡，滑坡詳細信息見表3，滑坡識別結果分布見圖7，局部滑坡識別結果見圖8。

圖8 部分滑坡識別結果圖Fig.8 Some landslide recognition results

表3 研究區(qū)滑坡識別結果Table 3 Landslide identification in the study area

根據(jù)表3、圖7分析，識別出的滑坡有一半以上是由于劇烈的人類工程活動誘發(fā)導致，85%的滑坡類型為土質(zhì)滑坡，在強降雨的情況下更容易出現(xiàn)失穩(wěn)，結合H1、H5、H9的時序形變分析，在6～9月份形變趨勢處于急劇增加的態(tài)勢，是因為該期間迎來降雨高峰期，雨水富集在松散土層和相對防水的基巖接觸面之間，增加了土壤的自重，同時軟化了滑帶土，使其內(nèi)聚力和內(nèi)摩擦角減小，使得坡體的穩(wěn)定性降低，為土層滑坡的形成奠定了基礎。

圖7 研究區(qū)滑坡識別分布圖Fig.7 Distribution map of landslide identification in study area

4.2.2 驗證分析

根據(jù)表3中潛在滑坡的高程值和坡度值進行分析，95%以上的滑坡分布在高程大于1 000 m，坡度大于20°的區(qū)域，借助前人對水城縣滑坡分布規(guī)律特征的研究成果進行對比驗證［20］，表明識別的潛在滑坡與水城縣滑坡分布規(guī)律特征相吻合，在一定程度上說明了識別結果的可靠性。

此外，通過貴州省地質(zhì)災害排查項目及野外實地調(diào)查獲得了水城縣的22處滑坡，用于對本文識別的潛在滑坡進行驗證。驗證結果：識別獲得的14處潛在滑坡，有11處潛在滑坡位置和范圍與野外調(diào)查的滑坡較為吻合（占識別結果的78.6%），有3處潛在滑坡范圍與野外調(diào)查的結果存在一定偏差（占識別結果的21.4%），部分識別結果見圖8，另有8處滑坡未能識別。

根據(jù)圖7，圖8可知，對于成功識別的滑坡點位置，大多是位于坡度和形變程度較大的區(qū)域。未能成功識別的8處滑坡點，但也均監(jiān)測到了部分形變區(qū)域，其未能識別的主要原因是：這些區(qū)域大多處于雷達不可見和不敏感區(qū)域，獲得有效形變點較少；部分地區(qū)植被覆蓋嚴重，為避免監(jiān)測結果受大面積植被因素的影響，剔除了該區(qū)域的形變點，減少了該區(qū)域的有效時序形變點位提取量，導致無法準確識別。

5 結論與討論

5.1 結論

本文以滑坡多發(fā)的水城縣為研究區(qū)，探討了一種基于InSAR技術對潛在滑坡識別中的應用。使用歐空局提供的Sentinel－1 A SAR數(shù)據(jù)，利用改進的SBAS-InSAR方法獲取地表坡度方向形變信息，通過Anselin Local Moran′s I進行聚類分析，獲得高異常值區(qū)域，并結合滑坡的形狀、分布、微地貌特征成功識別出了14處潛在滑坡。將識別的潛在滑坡與野外勘探資料相核實，驗證了該時序InSAR用于滑坡隱患識別的有效性及優(yōu)勢。

本研究受限于單一的降軌數(shù)據(jù)，對地形可視性差的區(qū)域無法進行有效識別，因此無法完全提取研究區(qū)內(nèi)潛在滑坡體。文中是對整個縣城的進行時序監(jiān)測，在處理數(shù)據(jù)上所花費時長較大，為了更快速的獲取識別結果，可利用升降軌數(shù)據(jù)直接對滑坡高易發(fā)區(qū)進行監(jiān)測識別，從而有針對性地調(diào)查和監(jiān)測高概率潛在的滑坡范圍，達到防止滑坡災害發(fā)生的目的。

5.2 討論

由于受到SAR系統(tǒng)監(jiān)測方式的制約，在對西南高山區(qū)域的滑坡識別應用中，導致可探測的形變區(qū)域面積較少，難以全方位對潛在的地質(zhì)災害隱患進行監(jiān)測識別，限制了通過InSAR技術對潛在地質(zhì)災害識別的應用。本研究獲取數(shù)據(jù)類型單一、有效觀測區(qū)覆蓋范圍較小，導致部分滑坡點未能進行監(jiān)測，因此可以選用聯(lián)合升軌和降軌觀測減少觀測盲區(qū)，提升滑坡的有效探測率。

現(xiàn)有的民用雷達衛(wèi)星通常工作于X波段（波長3.1cm）、C波段（波長5.6cm）或L波段（波長2.3.6cm），但是能大量獲取的免費SAR數(shù)據(jù)為C波段的Sentinel－1 A數(shù)據(jù)。對于InSAR形變測量來說，雷達波長決定測量靈敏度，即波長越短，監(jiān)測精度越高。但是在高山峽谷區(qū)域，X波段和C波段的波長較短，難以穿透植被冠層，容易導致時間失相干；使用波段更長、穿透力更強的L波段進行觀測時，數(shù)據(jù)往往具有更好的相干性，在西南山區(qū)滑坡識別應用中更具有優(yōu)勢。