基于InSAR技術(shù)的三峽庫區(qū)巫山—奉節(jié)段潛在滑坡識別

徐 帥，王尚曉，牛瑞卿

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

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

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

本文以三峽庫區(qū)滑坡災(zāi)害多發(fā)的巫山—奉節(jié)段為研究區(qū)，充分利用PS-InSAR技術(shù)和SBAS-InSAR技術(shù)各自的優(yōu)勢，綜合考慮該區(qū)域滑坡的形成機(jī)制，通過統(tǒng)計(jì)區(qū)域歷史滑坡的發(fā)育規(guī)律，將形成滑坡的孕災(zāi)因子引入到InSAR技術(shù)中，實(shí)現(xiàn)了高密度、高相關(guān)性的形變點(diǎn)的生成，并成功篩選出概率高、置信度水平較高的潛在滑坡體的形變點(diǎn)，最后結(jié)合該形變點(diǎn)區(qū)域地形地貌特征和高分辨率光學(xué)影像識別出潛在滑坡的范圍，實(shí)現(xiàn)了研究區(qū)潛在滑坡的早期識別，用以指導(dǎo)滑坡災(zāi)害的預(yù)測預(yù)警和野外地質(zhì)災(zāi)害的實(shí)地調(diào)查。

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

1.1 研究區(qū)概況

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

1.2 SAR數(shù)據(jù)源

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

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

1.3 地形可視性分析

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

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

(1)

Visibility=R·Lv·Sh

(2)

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

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

根據(jù)SAR成像特點(diǎn)以及Cigna等在文獻(xiàn)[17]中的大量統(tǒng)計(jì)試驗(yàn)可得：當(dāng)R≥sin(θ)，且背向雷達(dá)衛(wèi)星時(shí)，地形坡度大于視線入射角余角的區(qū)域?yàn)橹鲃雨幱埃划?dāng)0≤R≤sin(θ)時(shí)，出現(xiàn)雷達(dá)影像的透視收縮現(xiàn)象；當(dāng)-1≤R<0，且面向雷達(dá)衛(wèi)星時(shí)，地形坡度大于視線入射角的區(qū)域，出現(xiàn)雷達(dá)影像的主動疊掩現(xiàn)象；當(dāng)R≥sin(θ)，且非主動陰影和疊掩區(qū)時(shí)，為地形可視性較好的區(qū)域，在這個(gè)區(qū)域內(nèi)InSAR測量技術(shù)能夠獲取有效的地形監(jiān)測結(jié)果。因此，本文定義：R≥sin(θ)(Sentinel-1A升軌視線入射角θ=39.6°，即R≥0.56)且非主動陰影和疊掩區(qū)為地形可視性較好的區(qū)域，從而提高了干涉形變點(diǎn)的可信度。

2 研究方法

2.1 結(jié)合PS-InSAR技術(shù)的GCP點(diǎn)選擇

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

振幅離差(DΔA)是依據(jù)SAR影像上的同一位置像元在時(shí)間序列上強(qiáng)度值的離散特性進(jìn)行點(diǎn)目標(biāo)的識別，其計(jì)算公式如下：

(3)

式中：σΔA為干涉圖幅度的標(biāo)準(zhǔn)差；μA為干涉圖幅度的均值。

定義相干性系數(shù)γ的計(jì)算公式為

(4)

式中：M、S分別為構(gòu)成干涉對的SAR影像，其中M為主影像，S為輔影像；“*”為復(fù)數(shù)共扼;m、n分別為選擇窗口的大??；i、j則為像元的坐標(biāo)值。

2.2 研究區(qū)滑坡發(fā)育的統(tǒng)計(jì)規(guī)律

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

2.3 形變點(diǎn)空間域的異常值分析和聚類

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

表1 滑坡孕災(zāi)因子各等級的滑坡數(shù)及其對滑坡的概率貢獻(xiàn)值

顯著性的高值(熱點(diǎn))和低值(冷點(diǎn))的空間聚類以及形變數(shù)據(jù)集范圍內(nèi)的高異常值和低異常值，自動聚合形變點(diǎn)異常數(shù)據(jù)，確定出適當(dāng)?shù)目臻g分析范圍，并糾正數(shù)據(jù)的空間依賴性，計(jì)算出形變聚類點(diǎn)置信水平p值和得分，其相關(guān)計(jì)算公式如下：

(5)

(6)

(7)

(8)

3 研究結(jié)果與分析

3.1 研究區(qū)潛在滑坡點(diǎn)篩選

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

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

(9)

圖2 沿雷達(dá)視線方向上的平均形變速率圖(a)、置信水平在95%以上的作為潛在滑坡點(diǎn)的形變點(diǎn)空間 分布圖(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)為各孕災(zāi)因子對形變點(diǎn)作為潛在滑坡點(diǎn)的累計(jì)概率貢獻(xiàn)值;P(Ai)為每個(gè)獨(dú)立孕災(zāi)因子對滑坡的概率貢獻(xiàn)值(即1/n);P(B|Ai)表示每個(gè)孕災(zāi)因子類別在發(fā)生滑坡事件下的概率值。

3.2 研究區(qū)潛在的滑坡提取

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

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

4 結(jié) 論

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

(1) 聯(lián)合PS-InSAR和SBAS-InSAR兩種技術(shù)的優(yōu)勢，獲取了高密度、相干性較好的研究區(qū)形變結(jié)果，共計(jì)85 439個(gè)形變點(diǎn)，沿LOS方向的平均形變速率值的范圍為-29.83～50.53 mm/a。

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

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

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