基于InSAR技術的三峽庫區巫山—奉節段潛在滑坡識別

徐 帥，王尚曉，牛瑞卿

(1.中國地質大學(武漢)地球物理與空間信息學院，湖北 武漢 430074；2.中國地質調查局武漢地質調查中心(中南地質科技創新中心)，湖北 武漢 430205)

滑坡是指斜坡上的巖土體在自然條件和人類工程活動等因素的影響下，受重力作用沿貫通破壞面整體或者分散地順坡向下滑動的地質災害現象，同時也包括處于不穩定狀態，可能演化成滑坡的斜坡體[1]。常規的滑坡調查方法主要包括野外調查、目視解譯航片等，采用這些方法對典型的大型滑坡進行識別與監測，往往會花費大量的人力和財力，且多數情況下不能及時、準確地確定滑坡范圍，而且對滑動不明顯、潛在的滑坡的監測和識別能力不足[2-3]。合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)作為一種新型的空間地表測量技術，不受自然天氣條件的影響，能夠全天候、全天時地對地表進行觀測[4]。目前該技術已日漸成熟，并運用于由地震、火山、地面沉降、滑坡等造成的地表形變的測量，其不僅能夠獲取地表發生的微小形變，精度可達到厘米級甚至是毫米級，而且具有分辨率高、時效性高、覆蓋范圍廣等優勢，已經成為國內外學者研究的熱點[5-9]。

目前，將InSAR技術應用于滑坡災害的研究主要集中在對特定研究區域已知的單體滑坡進行滑坡體的變形監測和形變特征的分析[10-14],也有不少學者嘗試將該技術用于區域滑坡的早期識別以及潛在滑坡范圍的確定[15-16]，但都是基于形變點速率的數學統計分析，未考慮到滑坡的形成是復雜的孕災因子和外界致災因子綜合影響的過程，因此在滑坡的識別精度和潛在滑坡范圍的確定方面仍存在一定的問題。另外，常規的永久散射體干涉測量(Permanent Scatterer InSAR,PS-InSAR)技術和小基線集干涉測量(Small Baseline InSAR，SBAS-InSAR)技術有一定的局限性。PS-InSAR技術需要的影像數據較多，可以不用過多考慮由于時間基線和空間垂直基線過大引起的去相關問題，能有效地提取非線性形變信息并消除大氣相位的影響，但是處理方法較復雜，對數據的時間覆蓋范圍要求高；SBAS-InSAR技術可以極大地增加相干影像對的數目，在估計出平均形變速率之后，利用干涉解纏后的相位，通過空間域的低通濾波和時域的高通濾波去除大氣效應的影響，最大限度地解決了時間、空間失相關和大氣效應對InSAR測量的影響。在軌道精煉步驟中，引入穩定的地面控制點(Ground Control Point，GCP)可以有效地消除平地效應。然而，在選擇GCP點時，很難找到質量好的點，而且也不可能在未知區域建立大量可識別的人工穩定性目標，導致在SBAS處理流程中過于主觀性，沒有定量的標準去評價所選的GCP點，最終干涉結果不理想，生成的形變點密度太低，導致沒有足夠的形變信息。

本文以三峽庫區滑坡災害多發的巫山—奉節段為研究區，充分利用PS-InSAR技術和SBAS-InSAR技術各自的優勢，綜合考慮該區域滑坡的形成機制，通過統計區域歷史滑坡的發育規律，將形成滑坡的孕災因子引入到InSAR技術中，實現了高密度、高相關性的形變點的生成，并成功篩選出概率高、置信度水平較高的潛在滑坡體的形變點，最后結合該形變點區域地形地貌特征和高分辨率光學影像識別出潛在滑坡的范圍，實現了研究區潛在滑坡的早期識別，用以指導滑坡災害的預測預警和野外地質災害的實地調查。

1 研究區概況與SAR數據源

1.1 研究區概況

研究區位于三峽庫區巫山—奉節段(見圖1)，地理位置為109°33′E～110°11′E、30°45′N～23°28′N之間，長江橫貫東西。該地區地處中國地形第二階梯向第三階梯的過渡地帶，是川東褶皺與鄂西山地會合部位，以中、低山侵蝕峽谷地貌為主，地形起伏較大，地質條件復雜，地層巖性以二疊系下統、三疊系中統巴東組、三疊系下統大冶組和嘉陵江組、志留系中統羅惹坪群第三段以及第四系為主，巖層軟硬相間，次級褶皺和斷裂褶皺發育，極易孕育滑坡等地質災害。

1.2 SAR數據源

本文采用的數據源是歐洲空間局新一代C波段的中高分辨率衛星哨兵一號A星(Sentinel-1A)SAR數據，選取2017年3月—2018年11月共48景漸進式掃描對地觀測(Terrain Observation by Progressive Scans SAR，TOPSAR)成像模式下的單視復數干涉寬模式(IW)升軌數據，入射角為33.9°，重訪問周期為12 d，空間分辨率為5 m(R)×20 m(A)，測繪帶幅寬為250 km，極化方式為VV(下載網址為：https://vertex.daac.asf.alaska.edu/),并使用衛星成像21 d后的精密軌道數據進行軌道矯正，定位精度優于5 cm。數字高程模型(DEM)為30 m的SRTM DEM(下載網址為：http://gdex.cr.usgs.gov/gdex/)。

圖1 研究區概況及歷史滑坡位置圖Fig.1 Overview of the study area and historical landslide location map

1.3 地形可視性分析

InSAR技術應用于區域的形變監測時，往往由于地形的影響而產生疊掩和陰影，對落入該區域的形變點造成較大的干擾，因此需要對其進行地形可視性分析。本文采用Cigna等提出的R指數計算方法[17]和Notti等提出的疊掩和陰影計算方法[18]，并結合SAR傳感器采集數據的幾何參數、視線(light of sight，LOS)方向的角度以及地形的坡度和坡向，對研究區SAR數據的地形可視性進行分析，提取出干擾嚴重的區域，用于掩膜最終落入該區域的不可信的點。具體計算公式如下：

R=sin[θ-β·sin(A)]

(1)

Visibility=R·Lv·Sh

(2)

上式中：R為地形效應綜合指數；θ為視線入射角(°);β為地形坡度(°)；當升軌數據時，A=α-ε[其中,α為地形坡向(°)，ε為衛星飛行方位向與北方向的夾角(°)]；Lv為研究區的疊掩；Sh為研究區的陰影；Visibility為SAR數據的地形可視性。

先采用ArcGIS山體陰影模型和DEM，輸入參數為視線(LOS)方位向(γ)和入射角的余角(90°-θ)，計算出研究區陰影(Sh),采用Notti等提出的疊掩計算方法[18]，設置參數為視線(LOS)方位向的補角(γ+180°)和入射角(θ)，計算出研究區的疊掩(Lv);然后利用重分類將計算得到的陰影(Sh)和疊掩(Lv)分別設為0，其他的值設為1。

根據SAR成像特點以及Cigna等在文獻[17]中的大量統計試驗可得：當R≥sin(θ)，且背向雷達衛星時，地形坡度大于視線入射角余角的區域為主動陰影；當0≤R≤sin(θ)時，出現雷達影像的透視收縮現象；當-1≤R<0，且面向雷達衛星時，地形坡度大于視線入射角的區域，出現雷達影像的主動疊掩現象；當R≥sin(θ)，且非主動陰影和疊掩區時，為地形可視性較好的區域，在這個區域內InSAR測量技術能夠獲取有效的地形監測結果。因此，本文定義：R≥sin(θ)(Sentinel-1A升軌視線入射角θ=39.6°，即R≥0.56)且非主動陰影和疊掩區為地形可視性較好的區域，從而提高了干涉形變點的可信度。

2 研究方法

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

穩定的地面控制點(GCP)可以有效地估算殘余相位，進而消除平地效應，提高結果的可靠性。本文首先利用PS-InSAR技術的優勢，依次通過設置振幅離差閾值(不大于0.25)、相干性閾值(不小于0.70)以及空間域的低通濾波和時域的高通濾波，減小干涉相位中的大氣延遲相位誤差，獲取更加可靠的穩定散射點的形變信息；然后，設置形變速率閾值(±1 mm/a)，篩選出相對穩定的PS點；最后，在保證了高相干性和較好穩定性的基礎上，進行地理坐標系到SAR斜距坐標系的轉換，作為SBAS處理流程中的GCP點，避免了選點的盲目性。

振幅離差(DΔA)是依據SAR影像上的同一位置像元在時間序列上強度值的離散特性進行點目標的識別，其計算公式如下：

(3)

式中：σΔA為干涉圖幅度的標準差；μA為干涉圖幅度的均值。

定義相干性系數γ的計算公式為

(4)

式中：M、S分別為構成干涉對的SAR影像，其中M為主影像，S為輔影像；“*”為復數共扼;m、n分別為選擇窗口的大小；i、j則為像元的坐標值。

2.2 研究區滑坡發育的統計規律

滑坡的形成是受復雜的內外因素共同作用的結果，因此僅僅使用InSAR技術獲得地表形變信息來作為滑坡災害識別與監測的依據，是不具有較強說服力的。由于同一區域內滑坡形成的外部因素大致相同，而且受人為干擾嚴重，因此本文嘗試通過統計與分析研究區的歷史滑坡數據，歸納出該區域滑坡的孕災因子變化規律[19]，并確定各孕災因子不同等級的滑坡數及其對滑坡的概率貢獻值，進而對干涉生成的形變點進行潛在滑坡的概率估算。從三峽庫區地質災害防治工作指揮部1∶10 000災害地質圖數據庫中獲取了研究區范圍內103處歷史滑坡的相關數據，該區域以土質滑坡為主，分別統計出滑坡孕災因子各等級的滑坡數及其對滑坡的概率貢獻值，見表1。

2.3 形變點空間域的異常值分析和聚類

通過SBAS-InSAR技術流程生成大量形變干涉點，采用 Anselin Local Moran’sI指數進行形變點空間域的異常值分析和聚類[20]，識別出具有統計

表1 滑坡孕災因子各等級的滑坡數及其對滑坡的概率貢獻值

顯著性的高值(熱點)和低值(冷點)的空間聚類以及形變數據集范圍內的高異常值和低異常值，自動聚合形變點異常數據，確定出適當的空間分析范圍，并糾正數據的空間依賴性，計算出形變聚類點置信水平p值和得分，其相關計算公式如下：

(5)

(6)

(7)

(8)

3 研究結果與分析

3.1 研究區潛在滑坡點篩選

根據研究區形變點空間域的異常值分析和聚類處理結果，研究區內共生成了85 439個高相干性的形變點，沿LOS方向的平均形變速率值的范圍為-29.83～50.53 mm/a[見圖2(a)]。形變速率的正負號表示形變點相對于衛星傳感器視線的運動方向，正值表示形變點沿雷達視線方向向靠近傳感器的方向移動，負值表示形變點沿雷達視線方向向遠離傳感器的方向移動，形變速率的絕對值則代表時間域內平均形變速率的大小。沿LOS方向的形變往往不能真實地反映地表變形情況，因此本文利用視坡夾角及余弦定理將沿雷達視線方向的形變速率轉化為沿斜坡坡度方向的形變速率，進而客觀地反映沿斜坡面的地表形變信息[21]。當雷達視線方向與斜坡面夾角接近90°時，其余弦值接近0，此時沿斜坡面的形變速率趨于無窮大。為了避免出現絕對值極大的異常值，采用Herrera等提出的以余弦值±0.3為合理的閾值，即沿斜坡面的形變速率不能大于沿雷達視線方向形變速率的3.33倍[22]。沿斜坡面的形變速率為正值或者負值分別代表了地表點移動方向沿坡面向上或者向下。依據滑坡體的運動規律，雖然在斜坡坡腳可能出現垂直方向上的正向移動，但水平方向上的位移方向仍應該與沿坡面向下方向保持一致，因此剔除沿斜坡面的形變速率為正值的形變點。

本文選取地形可視性較好范圍內的形變點，結合研究區統計的7個滑坡孕災因子的變化規律以及各孕災因子對滑坡概率貢獻的值，假設這7個孕災因子對該區域滑坡的概率貢獻相同，依據全概率公式(9)，計算出所選形變點的累計概率貢獻值P，即各孕災因子對形變點作為潛在滑坡點的累計概率貢獻值P(B)，并篩選出累計概率貢獻值大于50%的點作為潛在滑坡點，對這些形變點空間域進行異常值分析和聚類，通過設置p值的大小，選擇出置信水平在95%以上的形變點，其空間范圍分布見圖2(b)。

(9)

圖2 沿雷達視線方向上的平均形變速率圖(a)、置信水平在95%以上的作為潛在滑坡點的形變點空間 分布圖(b)、圖2(b)的局部放大圖(c)Fig.2 Average rate image of deformation along the line of sight of the SAR(a)，Spatial distribution map of deformation points as potential landslides with confidence levels above 95%(b)，Partial enlargement map (c)

式中:P(B)為各孕災因子對形變點作為潛在滑坡點的累計概率貢獻值;P(Ai)為每個獨立孕災因子對滑坡的概率貢獻值(即1/n);P(B|Ai)表示每個孕災因子類別在發生滑坡事件下的概率值。

3.2 研究區潛在的滑坡提取

結合潛在滑坡點區域的地形地貌特征以及空間域異常值分析和聚類處理結果，在沿雷達視線方向上的平均形變速率圖[見圖2(a)]上標定出15處概率高、置信水平較好的疑似潛在滑坡范圍[見圖2(b)]，借助多時相的高分辨率光學影像數據和Google Earth影像數據，對這些潛在的滑坡范圍做進一步分析，發現有3處(編號為8、9、13)位于歷史滑坡上或者附近，這些區域斜坡有明顯的地表形變，存在不穩定狀況，出現再次滑動的概率較高；有5處(編號為3、4、10、12、15)是人工修路或者建設房屋造成的地表形變；還有7處(編號為1、2、5、6、7、11)地表有不同程度的破壞，發育成滑坡的概率較大。人類工程活動和地表季節性運動等因素也會影響研究區域內潛在滑坡識別的概率，且由于雷達的可視性有限，并不能完全識別出所有的潛在高概率滑坡，本文利用高分一號影像識別出9處新的潛在滑坡范圍，見圖3。

圖3 利用高分一號影像識別出的潛在滑坡范圍標定圖Fig.3 Calibration map of potential landslides range based on GF-1 image

4 結 論

本文在滑坡災害多發的三峽庫區巫山—奉節段，利用InSAR技術對潛在的滑坡災害進行了早期識別研究，探索了合理的干涉測量技術，并結合滑坡孕災因子和區域地形地貌條件，從滑坡的成因機制方面對識別出的潛在滑坡點進行概率估計，成功識別出概率高、置信水平較好的潛在滑坡范圍，可為該區域滑坡災害的預測預警和野外地質災害的實地調查提供科學依據，得到如下主要結論：

(1) 聯合PS-InSAR和SBAS-InSAR兩種技術的優勢，獲取了高密度、相干性較好的研究區形變結果，共計85 439個形變點，沿LOS方向的平均形變速率值的范圍為-29.83～50.53 mm/a。

(2) 充分考慮了SAR的地形可視性，計算了地形綜合效應指數R，篩選出R≥0.56且非主動陰影和疊掩區為地形可視性較好的區域，從而提高了干涉形變點的可信度。

(3) 引入了滑坡災害的7個孕災因子(地形坡度、地形坡向、植被覆蓋度、地形濕度指數、斜坡形態、距離斷層距離和地層巖性)，統計分析各孕災因子對研究區滑坡概率的貢獻值，采用全概率公式統計分析每個形變點作為潛在滑坡點的概率貢獻值，并通過形變點的Anselin Local Moran’sI指數值對作為潛在滑坡概率貢獻值較高的形變點進行空間域異常值分析和聚類處理，獲得15處概率高、置信水平較好的疑似潛在滑坡范圍，再借助區域地形地貌特征和高分辨率光學影像識別出研究區高概率的9處新的潛在滑坡。

本研究受限于單一的升軌數據，對地形可視性差的區域無法進行有效識別，對于未知區域，由于缺乏地面的監測站點，因此無法完全提取出區域內所有的高概率潛在滑坡體。本文所做的研究并不是追求更高的形變數據精度，而是通過考慮形變數據的異常情況和孕災因子對滑坡的概率貢獻，獲取可能發育成潛在滑坡的概率，進而有針對性地調查和監測高概率潛在的滑坡范圍，以達到防治滑坡災害發生的目的。