雷達數據輔助下的滑坡易發性評價

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

中國地質大學地球物理與空間信息學院，武漢 430074

0 引言

三峽工程作為人類重大工程活動對庫區地質環境影響巨大，滑坡作為典型的地質災害在庫區屢屢發生，對當地居民生命財產安全和三峽大壩正常運營構成了嚴重危害。因此，對三峽庫區范圍進行滑坡災害的易發性評價意義重大。

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

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

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

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

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

1.1巖土體含水量反演

1.1.1 水-云模型

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

(1)

(2)

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

(3)

1.1.2 巖土體含水量反演模型

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

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

(4)

(5)

1.2 邏輯回歸

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

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

(6)

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

2 研究區概況及數據源

2.1 研究區概況

研究區位于長江三峽庫區秭歸縣境內，地理坐標為110°30′00″E—110°51′00″E， 30°55′00″N—31°04′00″N，面積約為293 km2，其中水系面積約41 km2。研究區位于湖北省西部、三峽庫區庫首位置(圖1)，地勢西南高東北低，海拔80～1 500 m，沿江兩岸高山林立，為峽谷地貌；氣候屬亞熱帶大陸季風氣候，雨量充沛。研究區內滑坡災害頻發。

2.2 數據源

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

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

3 巖土體含水量提取

3.1影像預處理

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

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

a. Sentinel-1A雷達影像；b.Lansat8影像。圖2 研究區遙感數據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 水-云模型應用

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

Cveg=1.44INDW2+1.36INDW+0.34；

(7)

(8)

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

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

4 滑坡易發性評價

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

4.1 評價因子選取

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

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

4.2 因子歸一化處理

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

表2 各因子Pearson相關系數

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

4.3 建立滑坡易發性評價模型

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

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

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

4) 滑坡易發性評價。利用建立的模型對研究區全部柵格單元進行計算，結果輸出即為滑坡易發性指數。

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

5 結果分析和驗證

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

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

易發性指數易發性分區所占面積比/%滑坡所占面積比/%0.8～1.00.5～0.80.2～0.50.0～0.2高易發區中易發區低易發區不易發區9.714.916.359.153.635.19.61.7

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

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

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

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

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

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

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

6 結論與展望

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

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

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