雷達(dá)數(shù)據(jù)輔助下的滑坡易發(fā)性評價

趙金童,牛瑞卿,姚 琦,武雪玲

中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院，武漢 430074

0 引言

三峽工程作為人類重大工程活動對庫區(qū)地質(zhì)環(huán)境影響巨大，滑坡作為典型的地質(zhì)災(zāi)害在庫區(qū)屢屢發(fā)生，對當(dāng)?shù)鼐用裆敭a(chǎn)安全和三峽大壩正常運(yùn)營構(gòu)成了嚴(yán)重危害。因此，對三峽庫區(qū)范圍進(jìn)行滑坡災(zāi)害的易發(fā)性評價意義重大。

隨著遙感、地理信息系統(tǒng)和計算機(jī)技術(shù)的發(fā)展，對于滑坡的易發(fā)性評價，國內(nèi)外學(xué)者已經(jīng)開展了多方面的研究。例如：Carrara[1]利用回歸分析對意大利山區(qū)滑坡進(jìn)行易發(fā)性評價；Lee[2]基于GIS技術(shù)通過邏輯回歸模型進(jìn)行滑坡易發(fā)性預(yù)測；Ballabio 等[3]采用支持向量機(jī)模型對意大利北部區(qū)域進(jìn)行了滑坡易發(fā)性制圖；晏同珍等[4]、殷坤龍[5]多次對三峽庫區(qū)進(jìn)行滑坡災(zāi)害空間預(yù)測；牛瑞卿等[6]利用粗糙集-支持向量機(jī)模型實現(xiàn)了三峽庫區(qū)秭歸到巴東段滑坡易發(fā)性分區(qū)；彭令等[7]、武雪玲等[8]分別在多源數(shù)據(jù)支持下得到了三峽庫區(qū)滑坡穩(wěn)定性分區(qū)圖。雖然上述研究已經(jīng)使滑坡空間預(yù)測達(dá)到了一定的精度，但其中選取的濕度指數(shù)因子都只是間接地反應(yīng)了巖土體含水量，再通過巖土體含水量間接反應(yīng)地下水信息，而缺少因子來直接表示地表水文狀況。

隨著微波遙感技術(shù)的持續(xù)發(fā)展，各種適用模型不斷改進(jìn)，大面積反演巖土體含水量已經(jīng)成為了可能。Paloscia等[9]研究表明，基于神經(jīng)網(wǎng)絡(luò)算法利用Sentinel-1影像反演巖土體含水量(soil moisture content，CSM)是可行的；何連等[10]利用多時相Sentinel-1數(shù)據(jù)實現(xiàn)了農(nóng)田地表巖土體水分的反演；Oh等[11]根據(jù)不同極化方式下粗糙度和介電常數(shù)與后向散射系數(shù)的相關(guān)關(guān)系提出了反演巖土體含水量的經(jīng)驗?zāi)Ｐ?；Fung等[12]提出積分方程模型(integrated equation model，IEM)，該模型基于輻射傳輸方程建立地表散射模型；Fung等[13]對IEM進(jìn)行了改進(jìn)，得到一種高級積分方程模型(advanced integrated equation model，AIEM)； Ulaby等[14]提出了基于微波輻射傳輸方程一階解的微波散射模型(michigan microwave canopy scattering model，MIMICSM)，主要針對植被覆蓋的粗糙地面，是目前應(yīng)用最多的植被覆蓋地區(qū)散射模型，但其較為復(fù)雜；Attema等[15]提出的水-云模型將植被看作各項散射均質(zhì)體，模型簡單實用。

Hassaballa等[16]分別利用RADARSAT-1雷達(dá)數(shù)據(jù)和Landsat-7 ETM+光學(xué)數(shù)據(jù)反演吉隆坡Bukit Antarabangs滑坡地區(qū)地表巖土體含水量，研究其與該地區(qū)歷史性滑坡的空間相關(guān)性，發(fā)現(xiàn)利用雷達(dá)數(shù)據(jù)反演得到的含水量與滑坡空間相關(guān)性的受試者工作特征曲線下的面積(area under curve，AUC)為0.7且高于光學(xué)數(shù)據(jù)，證明地表巖土體含水量高的地區(qū)更易發(fā)生滑坡且雷達(dá)數(shù)據(jù)相對光學(xué)數(shù)據(jù)能更好地提取巖土體含水量；Brocca等[17]引入散射計(ASCAT)數(shù)據(jù)衍生得到的巖土體含水量作為水文因子來預(yù)測意大利中部地區(qū)Torgiovannetto滑坡的張力裂縫，結(jié)果表明其相關(guān)系數(shù)高于采用降雨量作為水文因子的預(yù)測結(jié)果，也就意味著巖土體含水量與滑坡的發(fā)生具有較強(qiáng)的相關(guān)性。

本文在前人研究的基礎(chǔ)上，選取三峽庫區(qū)秭歸段作為研究區(qū)，利用Sentinel-1雷達(dá)數(shù)據(jù)通過水-云模型去除植被影響得到裸露地表后向散射系數(shù)，采取同時相同范圍不同極化方式的Sentinel-1雷達(dá)數(shù)據(jù)獲取地表粗糙度參數(shù)，再利用前人在AIEM的基礎(chǔ)上建立的后向散射模型反演得到巖土體含水量；以此作為反映地表水及地下水信息的因子替代其他濕度指數(shù)，以實現(xiàn)研究區(qū)滑坡易發(fā)性評價。

1 巖土體含水量反演和邏輯回歸

1.1巖土體含水量反演

1.1.1 水-云模型

水-云模型(water-cloud model)由Attma等[15]于1978年提出，主要為了簡化植被覆蓋區(qū)的后向散射機(jī)制，提高模型適用性。該模型將植被層看作一個各項均質(zhì)散射體，忽略了植被層及地表之間的相互多次散射，將植被覆蓋地區(qū)總的后向散射簡單描述為兩部分：由植被層直接反射回的體散射項；經(jīng)植被層二次衰減后地表的后向散射項。模型表述如下：

(1)

(2)

γ2(θ)=exp[-2BCvegsecθ]。

(3)

1.1.2 巖土體含水量反演模型

決定裸露地表后向散射系數(shù)的參數(shù)主要是CSM和地表粗糙度，而地表粗糙度通常由均方根高度(s)在垂直方向相關(guān)長度(l)在水平方向上的組合表示。為了簡化模型參數(shù)， Zribi等[18]提出利用地表粗糙度參數(shù)Zs來表示粗糙度，Zs=s2/l。

李震等[19]通過高級積分方程模型的模擬數(shù)據(jù)，分別建立了同極化雷達(dá)后向散射模型和雙極化雷達(dá)數(shù)據(jù)粗糙度計算模型，然后將兩種模型結(jié)合即可求得CSM：

(4)

(5)

1.2 邏輯回歸

邏輯回歸模型(logistic regression model)是一種常見的機(jī)器學(xué)習(xí)模型，因其便捷高效，得以普遍應(yīng)用。因為滑坡只有發(fā)生或不發(fā)生兩種情況(一般1代表發(fā)生，0代表未發(fā)生)，本質(zhì)上它是一個二分類問題，所以二元邏輯回歸模型可以完美地解決滑坡易發(fā)性評價中的二分類問題[21]。

我們假設(shè)滑坡發(fā)生的概率為P，不發(fā)生的概率則為1-P。以P作為因變量，xi為自變量，建立回歸方程：

(6)

式中：α為常數(shù)；βm為滑坡各影響因子的邏輯回歸系數(shù)；xm即代表滑坡各影響因子的相對權(quán)重；m為影響因子數(shù)目。

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

2.1 研究區(qū)概況

研究區(qū)位于長江三峽庫區(qū)秭歸縣境內(nèi)，地理坐標(biāo)為110°30′00″E—110°51′00″E， 30°55′00″N—31°04′00″N，面積約為293 km2，其中水系面積約41 km2。研究區(qū)位于湖北省西部、三峽庫區(qū)庫首位置(圖1)，地勢西南高東北低，海拔80～1 500 m，沿江兩岸高山林立，為峽谷地貌；氣候?qū)賮啛釒Т箨懠撅L(fēng)氣候，雨量充沛。研究區(qū)內(nèi)滑坡災(zāi)害頻發(fā)。

2.2 數(shù)據(jù)源

本文主要采用的數(shù)據(jù)源有：1)圖2a 為Sentinel-1A IW(interferometric wide swath)模式下 Level-1 地距影像(ground range detected，GRD)數(shù)據(jù)，編號為S1A_IW_GRDH_1SDV_20150826T103524_20150826T103557_007433_00A3CA_0028，用于提取地表巖土體含水量；2)圖2b 為Landsat-8 OLI 衛(wèi)星影像一景，編號為LC81250392015216LGN00，用來提取地表覆被信息；3)1∶5萬比例尺地形圖，用于提取地形、地貌等信息；4)1∶5萬及1∶20萬比例尺地質(zhì)圖，用于提取地層巖性、構(gòu)造等信息；5)其他歷史滑坡資料和野外調(diào)查資料，用于已發(fā)生滑坡的解譯。

Sentinel-1A作為歐空局(European space agency ，ESA)哥白尼全球?qū)Φ赜^測計劃的首顆衛(wèi)星于2014年4月發(fā)射，主要用于取代Envisat衛(wèi)星提供全球范圍的高分辨率C波段雷達(dá)數(shù)據(jù)。本文選取的數(shù)據(jù)為Level-1 GRD產(chǎn)品，它經(jīng)過多視處理并經(jīng)WGS84橢球投影，只有幅度信息而沒有相位信息。Landsat-8 由美國國家航空航天局(national aeronautics and space administration, NASA)于2013年2月發(fā)射，其攜帶的OLI(operational land imager)陸地成像儀除了ETM+傳感器包含的7個波段外還增加了藍(lán)色波段和短波紅外波段。本文所使用的遙感數(shù)據(jù)參數(shù)信息見表1。

3 巖土體含水量提取

3.1影像預(yù)處理

對Landsat-8影像處理首先進(jìn)行輻射定標(biāo)、大氣校正等預(yù)處理，然后獲取歸一化差分水體指數(shù)(normalized difference water index，INDW)、歸一化植被指數(shù)(normalized difference vegetation index,INDV)和土地覆被等信息，主要通過ENVI軟件實現(xiàn)；對于Sentinel-1雷達(dá)數(shù)據(jù)的處理主要包括輻射校正、斑點濾波和幾何校正，目的是獲取植被覆蓋地表總的雷達(dá)后向散射系數(shù)，主要通過ESA開源軟件Sentinel application platform 實現(xiàn)。

圖1 研究區(qū)位置示意圖Fig.1 Location of the study area

a. Sentinel-1A雷達(dá)影像；b.Lansat8影像。圖2 研究區(qū)遙感數(shù)據(jù)Fig.2 Remote sensing data of the study area

影像類型成像時間分辨率/m模式入射角/(°)頻率/GHz極化方式Sentinel-12015-08-2620×20IW GRD32.95.405VH+VVLandsat-82015-08-0430×30OLI

注：IW GRD.干涉寬幅模式。

3.2 水-云模型應(yīng)用

植被含水量即為單位面積植被中所含水的質(zhì)量。作為水-云模型中一個重要的參數(shù)，其大面積提取比較困難，通常利用Landsat等陸地影像進(jìn)行反演。2004年Jackson[22]發(fā)現(xiàn)植被含水量Cveg和Gao[23]提出的歸一化差分水體指數(shù)INDW有很強(qiáng)的聯(lián)系，即：

Cveg=1.44INDW2+1.36INDW+0.34；

(7)

(8)

其中，RNIR、RSWIR分別代表Landsat-8影像的5、6波段的灰度值。

圖3反演得到的含水量只是淺層地表的巖土體含水量，但Uhlemann等[25]通過連續(xù)3年的地電監(jiān)測數(shù)據(jù)發(fā)現(xiàn)滑坡重新激活前淺層巖土體含水量顯著提升，同時鄧孺孺等[26]的研究也發(fā)現(xiàn)深層巖土體含水量與淺層巖土體含水量正相關(guān)；兩者研究表明淺層地表含水量也顯著影響著滑坡的發(fā)生且與深層巖土體含水量成正比。雖然深層巖土體含水量對滑坡穩(wěn)定性具有更重要的意義，但在現(xiàn)有技術(shù)手段下反演地表巖土體含水量可能是更現(xiàn)實的選擇。

4 滑坡易發(fā)性評價

由于本文數(shù)據(jù)來源多樣，坐標(biāo)系、投影不盡相同，首先需要對數(shù)據(jù)統(tǒng)一至相同投影坐標(biāo)系下。本文選取柵格單元作為評價單元，需將地質(zhì)、地表覆被等矢量數(shù)據(jù)轉(zhuǎn)為柵格數(shù)據(jù)，并將柵格單元大小統(tǒng)一至30 m×30 m。研究發(fā)現(xiàn)，土質(zhì)滑坡較巖質(zhì)滑坡更易受水的影響[27]，所以選取土質(zhì)滑坡作為研究對象。

4.1 評價因子選取

滑坡的發(fā)生主要受控制因素和影響因素兩大類因素的影響。控制因素指的是研究區(qū)的地形、地質(zhì)構(gòu)造等因素；影響因素指的是人類工程活動、地表覆被等因素。本文在前人研究的基礎(chǔ)上結(jié)合研究區(qū)的特點選取地質(zhì)(工程巖組指數(shù)(IPT)、斷層距離指數(shù)(IDF)、斜坡結(jié)構(gòu)指數(shù)(ISA))、地形地貌(數(shù)字高程(D)、地形粗糙指數(shù)(ITR)、地形表面紋理(ITST)、坡度(i))、地表覆被(歸一化植被指數(shù)(INDV)、土地利用指數(shù)(ILU))和水文(水系距離指數(shù)(IDR)、地形濕度指數(shù)(ITW)、巖土體含水量(CSM))四大類共12個評價因子。為了確保各因子間的相對獨立性，利用Pearson積差相關(guān)系數(shù)法對因子間相關(guān)性進(jìn)行評價(表2)，發(fā)現(xiàn)相關(guān)系數(shù)的絕對值均小于0.75，可以作為因子評價滑坡易發(fā)性。

圖3 巖土體含水量分布圖Fig.3 Map of soil moisture content

4.2 因子歸一化處理

因所選的評價因子處于不同類型、范圍和數(shù)量級下，歸一化處理可以提高數(shù)據(jù)處理效率，改善模型精度，所以對于連續(xù)型因子，我們采用極差標(biāo)準(zhǔn)化方法對其歸一化；對于離散型因子(包括斜坡結(jié)構(gòu)、土地利用和工程巖組)，我們根據(jù)已解譯出的滑坡柵格單元中各因子不同屬性占總滑坡柵格單元的比例作為各屬性的值，再采用極差標(biāo)準(zhǔn)化方法對各屬性值進(jìn)行歸一化，結(jié)果見表3。

表2 各因子Pearson相關(guān)系數(shù)

表3 離散型因子歸一化處理

4.3 建立滑坡易發(fā)性評價模型

基于邏輯回歸模型的滑坡預(yù)測主要思路就是選取合適的訓(xùn)練和測試樣本建立模型并驗證模型精度，進(jìn)而將模型應(yīng)用到整個研究區(qū)，得到研究區(qū)滑坡易發(fā)性分區(qū)圖。本文分別使用巖土體含水量，以及前人常用的地形濕度指數(shù)為地表濕度指數(shù)因子進(jìn)行易發(fā)性研究，并對不同濕度指數(shù)因子下易發(fā)性評價的精度進(jìn)行分析，驗證采用巖土體含水量因子進(jìn)行評價的可行性。主要步驟如下：

1) 劃分柵格單元。本文采用的數(shù)據(jù)基本上都是以柵格或矢量形式存在的影像數(shù)據(jù)，同時結(jié)合遙感影像的空間分辨率，將所有評價因子重采樣為30 m×30 m大小的柵格單元。經(jīng)統(tǒng)計可得，研究區(qū)共劃分為280 670個柵格單元，其中歷史土質(zhì)滑坡單元6 830個，未發(fā)生滑坡單元273 840個。

2) 選擇建模樣本。我們從所有的滑坡柵格單元中隨機(jī)抽取80%(5 464個柵格單元)和相同數(shù)量的未滑坡柵格單元作為訓(xùn)練樣本，剩下20%(1 366個柵格單元)滑坡柵格單元作為測試樣本。

4) 滑坡易發(fā)性評價。利用建立的模型對研究區(qū)全部柵格單元進(jìn)行計算，結(jié)果輸出即為滑坡易發(fā)性指數(shù)。

5) 結(jié)果分析。通過測試樣本對滑坡預(yù)測結(jié)果進(jìn)行評價，并對比分析采用巖土體含水量、地形濕度指數(shù)因子時的預(yù)測精度。

5 結(jié)果分析和驗證

對模型輸入整個研究區(qū)范圍的因子，輸出即為介于0.0～1.0之間的滑坡易發(fā)性指數(shù)，指數(shù)越大，滑坡越容易發(fā)生，反之滑坡越不容易發(fā)生。本文采用自然斷點法對滑坡易發(fā)性指數(shù)進(jìn)行分區(qū)，劃分為不易發(fā)區(qū)(0.0～0.2)、低易發(fā)區(qū)(0.2～0.5)、中易發(fā)區(qū)(0.5～0.8)、高易發(fā)區(qū)(0.8～1.0)，得到巖土體含水量因子下滑坡易發(fā)性分區(qū)圖(圖4)，并對各分區(qū)所占面積比(表4)進(jìn)行分析。

()內(nèi)為滑坡易發(fā)性指數(shù)。圖4 巖土體含水量因子下滑坡易發(fā)性分區(qū)圖Fig.4 Land slide susceptibility zoning map based on the CSM

易發(fā)性指數(shù)易發(fā)性分區(qū)所占面積比/%滑坡所占面積比/%0.8～1.00.5～0.80.2～0.50.0～0.2高易發(fā)區(qū)中易發(fā)區(qū)低易發(fā)區(qū)不易發(fā)區(qū)9.714.916.359.153.635.19.61.7

預(yù)測結(jié)果顯示，土質(zhì)滑坡多發(fā)生在河流沿岸，尤其在長江支流童莊河、青干河、吒溪河沿岸易發(fā)性指數(shù)較高，表明土質(zhì)滑坡受水的影響較大，基本符合現(xiàn)實滑坡分布狀況。

由表4可知，中高易發(fā)區(qū)面積僅占研究區(qū)總面積的24.6%，而區(qū)內(nèi)滑坡面積卻占到了總滑坡面積的88.7%；這也間接表明了我們采用巖土體含水量因子進(jìn)行滑坡易發(fā)性評價是可行有效的。

選擇研究區(qū)內(nèi)典型土質(zhì)滑坡新灘滑坡進(jìn)行預(yù)測結(jié)果驗證。根據(jù)已有歷史資料和實際調(diào)查資料，新灘滑坡位于屈原鎮(zhèn)西北2 km處長江北岸，與鏈子崖危巖體隔江相望，為古滑坡。最近一次滑動發(fā)生于1985年6月12日，滑坡體積約3 000萬m3，南北長約2 km，滑動面積約為1.1km2；由于預(yù)報及時，新灘鎮(zhèn)居民無一傷亡。在滑坡易發(fā)性分區(qū)圖中，該地區(qū)被預(yù)測為高易發(fā)區(qū)，如圖5所示。

為了檢驗2個不同濕度指數(shù)下滑坡的預(yù)測能力并進(jìn)行對比分析，我們利用測試樣本采用Chung等[28]提出的成功率驗證法測試2種模型的預(yù)測精度(圖6)。圖6中：橫軸代表滑坡易發(fā)性從大到小分成的100個區(qū)間，如橫軸10.0%處即表示滑坡易發(fā)性指數(shù)為90.0%～100.0%的區(qū)間；縱軸即為對應(yīng)易發(fā)性區(qū)間下實際滑坡累計發(fā)生的頻率，曲線下面積越大，證明預(yù)測效果越好。

由圖6可知，對于本文預(yù)測的滑坡易發(fā)性指數(shù)為80.0%～100.0%的區(qū)間(橫軸20.0%處)，即對應(yīng)易發(fā)性分區(qū)中的高易發(fā)區(qū)，采用巖土體含水量、地形濕度指數(shù)分別預(yù)測了64.6%、59.4%的實際滑坡落在高易發(fā)區(qū)中，采用巖土體含水量、地形濕度指數(shù)因子時曲線下面積即模型預(yù)測精度分別為80.2%、77.2%。按照滑坡易發(fā)性指數(shù)和滑坡累計發(fā)生頻率將地形濕度指數(shù)由數(shù)字高程模型中提取的坡度、流域面積計算得到，它主要考慮了地形對地表水的再分配作用來描述巖土體水分的空間分布。但巖土體水分的空間分布還受植被、巖土體特性和人類活動等的干擾，這也決定了地形濕度指數(shù)因子只能間接反應(yīng)巖土體含水量這一重要水文因子對研究區(qū)斜坡穩(wěn)定性的影響。而通過雷達(dá)數(shù)據(jù)可以直接反演出去除植被影響后的裸露地表巖土體含水量，它也與區(qū)域范圍內(nèi)地下水的分布狀況有很強(qiáng)的聯(lián)系。同時滑坡預(yù)測結(jié)果也證明：利用巖土體含水量因子比地形濕度指數(shù)因子有著更高的預(yù)測準(zhǔn)確率，意味著它作為直接反應(yīng)巖土體含水量的水文因子比其他濕度指數(shù)能更好地反應(yīng)地表及地下水信息，可以作為滑坡易發(fā)性評價的重要水文因子。

圖5 典型滑坡驗證Fig.5 Typical landslide validation

圖6 滑坡累計發(fā)生頻率與易發(fā)性指數(shù)曲線Fig.6 Relation between landslide accumulative frequency and prediction index

6 結(jié)論與展望

1)通過C波段Sentinel-1雷達(dá)數(shù)據(jù)反演獲取巖土體含水量，將其作為一種新引入的水文因子代替地形濕度指數(shù)來表征地表及地下水信息。在此基礎(chǔ)之上，再充分利用多源空間數(shù)據(jù)提取其他因子并對數(shù)據(jù)進(jìn)行歸一化處理，分別建立邏輯回歸模型，定量獲取滑坡易發(fā)性指數(shù)，得到研究區(qū)易發(fā)性分區(qū)圖。

2)本文利用雷達(dá)數(shù)據(jù)獲取巖土體含水量因子進(jìn)行滑坡預(yù)測的精度要高于傳統(tǒng)方法采用的地形濕度指數(shù)因子，證明巖土體含水量作為重要的水文因子引入到滑坡預(yù)測中是可行的并且精度較好。但由于巖土體含水量的反演受地形及地表覆被的影響較大，現(xiàn)有模型方法普適性較差，大面積反演水文因子仍存在較大困難。

3)由于條件限制，本文并未測定研究區(qū)實地巖土體含水量，如何相對準(zhǔn)確的大面積獲取巖土體含水量必然是后續(xù)研究的重點。隨著遙感技術(shù)的發(fā)展，滑坡空間預(yù)測的突破點就在于如何結(jié)合利用多源多類型的遙感數(shù)據(jù)。