萬 洋，郭 捷，馬鳳山，劉 佳，宋燁煒

（1.中國科學院地質與地球物理研究所/中國科學院頁巖氣與地質工程重點實驗室,北京 100029；2.中國科學院地球科學研究院,北京 100029；3.中國科學院大學,北京 100049）

0 引言

近年來，基于“一帶一路”發(fā)展戰(zhàn)略需求，對于南亞通道的建設開始逐漸完善。作為中國與尼泊爾密切往來的口岸，吉隆是中尼鐵路建設的必經之地。這里地質構造復雜、內外動力作用強烈，具有強地震風險、強地殼變形、強應力集中等特征，是重大地質災害頻發(fā)區(qū)域[1]。經沿線實地調查，該地區(qū)滑坡災害發(fā)生較頻繁，且滑坡具有突發(fā)性和不確定性。因此掌握滑坡的易發(fā)性對于滑坡預測和防災減災具有重要意義。

滑坡易發(fā)性評價指某一區(qū)域內現(xiàn)有自然條件下已經發(fā)生了和誘發(fā)因素下滑坡災害與空間分布的定性或定量關系[2?3]。目前多種滑坡災害易發(fā)性的評價方法中，大體可以總結為定性和定量兩類方法。早期基于專家經驗的現(xiàn)場分析和評價因子圖層疊加分析的方法過分依賴評價者的主觀意識，可比性較低。隨著地理信息系統(tǒng)技術、衛(wèi)星遙感技術和數(shù)字地形數(shù)據(jù)的全面推進，機器學習模型借助GIS 平臺逐步被引用到地質災害評價領域[4?6]。

MaxEnt 模型是以最大熵原理為理論，以JAVA 語言為支撐的機器學習模型。在1957年由杰恩斯首先提出最大熵原理，它的中心思想是，在只知道關于未知分布的部分約束條件時，應選取滿足這些約束條件且熵值最大的概率分布。Phillips 等在2006年提出了最大熵模型，專門用于生態(tài)建模和物質分布預測。該模型在2013年被引入滑坡災害的預測研究中。相比其他概率模型，如頻率比、證據(jù)權法、多標準決策分析和邏輯回歸等，最大熵模型既可以考慮滑坡與評價因子之間的關系，同時還能充分考慮所有因子本身之間的關系。相較于層次分析法，最大熵模型不需要專家經驗，在易發(fā)性分析時不需要加入主觀假設。比較同為機器學習模型的神經網絡，最大熵模型可以克服人工神經網絡的“黑箱”限制?？傊?，最大熵模型雖然數(shù)學推導較復雜，但模型適應性、穩(wěn)定性和靈活性更高，預測準確度也更高[7?16]。

文中使用最大熵模型對研究區(qū)滑坡易發(fā)性建模，然后通過ROC 曲線驗證模型的效果。并進一步對各評價因子的貢獻度和敏感閾值進行研究。為中尼交通廊道區(qū) 域的滑坡災害預警提供參考。

1 研究區(qū)概況

吉隆縣地處西藏自治區(qū)日喀則地區(qū)西南部，南面和西南面與尼泊爾相鄰，東臨聶拉木縣，北面與薩噶縣搭界，自古就是西藏通往南亞的交通要道[17]（圖1）。

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

研究區(qū)地處青藏高原西南部，全縣平均海拔在4 000 m 以上，地勢北高南低。按地貌成因可劃分為侵蝕剝蝕山地地貌和河谷侵蝕堆積地貌。位于岡底斯構造帶和北喜馬拉雅構造帶之間，新構造運動活躍，屬緩慢抬升區(qū)[18?19]。加之地震多發(fā)且強烈，誘發(fā)差異性抬升等特殊的第四紀地質環(huán)境，使之成為崩塌、滑坡、泥石流多發(fā)區(qū)。研究區(qū)出露地層主要為第四系松散堆積層與前震旦紀江東巖組變質巖和前震旦紀中酸性侵入巖。

研究區(qū)調查的G216 國道是吉隆縣通往吉隆鎮(zhèn)的唯一公路，總長86.61 km。范圍是以道路為中心向兩側擴展5 km 區(qū)域。經過對研究區(qū)的現(xiàn)場細致排查及室內遙感解譯補充，共確定169 處滑坡，如圖2（a）所示。調查沿途多為高山峽谷路段，地層多為風化破碎的抗侵蝕能力較弱的基巖，重力侵蝕活躍，切坡造成的大臨空面、高陡坡腳的斜坡體。沿線性工程兩側地表松散物質較多，坡腳均存在小范圍垮塌的可能。較為典型的滑坡如圖2（b）（c）。圖2（b）滑坡位于G216 國道K3470 處，滑坡體主體為土石混合體。人工開挖后，坡腳被削，形成牽引式滑坡。坡腳處修筑有45 m 的漿砌片石擋墻。在暴雨、強風作用下，坡面上塊石易發(fā)生滾落。圖2（c）本拉錯滑坡，位于吉隆鎮(zhèn)北側。坡腳緊靠交通線，剪出口沖過公路與護欄，已嚴重影響道路通行，目前斜坡處于破壞休止階段，但在地震、暴雨等因素情況下，會加速邊坡的破壞，影響道路安全。

圖2 研究區(qū)災害點分布及典型災害Fig.2 Distribution of disaster points and typical disasters in the study area

2 研究方法

2.1 評價因子選取

滑坡的發(fā)生受到多種因素的影響。文中在充分考慮了研究區(qū)尺度和特性后，選取了高程、坡向、坡度、斷層密度、河流密度、巖性、地震峰值加速度、植被覆蓋指數(shù)8 個評價因子（表1）。利用ArcGIS 對評價因子進行柵格化及重分類處理，生成圖3所示8 個圖層。各圖層柵格列數(shù)行數(shù)保持一致分別為2 547、5 594，柵格總數(shù)6 575 436 個。

圖3 研究區(qū)評價因子圖層Fig.3 Layer of evaluation factors in the study area

表1 評價因子選取及其來源Table 1 Selection of evaluation factors and sources

（1）高程可以直接反映研究區(qū)的地形起伏[20]。將地理空間數(shù)據(jù)云DEM 數(shù)據(jù)，通過ArcGIS 進行裁剪、拼接等處理，生成研究區(qū)DEM 圖；（2）坡向對斜坡的影響表現(xiàn)為向陽坡受日照時間較長，相應濕度較低，巖石表面風化程度較高，易形成松散堆積物，從而斜坡的危險性變高。通過ArcGIS 將DEM 圖層轉化得到坡向圖層；（3）坡度大小影響了滑坡體的移動速度和規(guī)模程度。通過ArcGIS 將DEM 圖層轉化得到坡度圖層；（4）斷裂構造使斷層帶及其附近一定范圍內的巖土體結構遭到破壞，降低了坡體的完整程度。斷裂帶越密集巖體也就越破碎，因此選用斷層帶的線密度作為一個評價因子。以現(xiàn)場調查為基礎地質資料為補充，利用ArcGIS 線密度計算工具計算得到；（5）研究區(qū)地表水系交錯，根據(jù)地質圖利用ArcGIS 線密度計算工具生成河流線密度作為滑坡易發(fā)性評價因子；（6）地層巖性決定著斜坡巖土體的強度，控制著地質災害發(fā)育和分布規(guī)律。根據(jù)西藏地區(qū)區(qū)域性地質圖及現(xiàn)場調查將研究區(qū)地層巖性劃分為四組，分別為：較堅硬變質巖巖組、堅硬巖漿巖巖組、堅硬沉積巖巖組、較軟弱沉積巖巖組；（7）地震峰值加速度（PGA）的大小表征地球內動力的強弱[21]。吉隆縣受地震的影響較大，該地區(qū)PGA為0.1～0.20g；（8）植被覆蓋指數(shù)用歸一化差分植被指數(shù)（NDVI）來定量表征。通常與植物蒸騰作用、光合作用有關，是指示植被生長狀態(tài)及植被覆蓋度的最佳指示因子[22]。將Landsat8 的數(shù)據(jù)1～7 波段加載到ArcGIS 中進行波段合成。然后計算NDVI，計算公式為：

式中：NIR——近紅外波段；

R——紅波段處的反射率值。

NDVI取值[1，?1]，負值表示地面覆蓋為云、水、雪等；0 表示有巖石或裸土等；正值表示有植被覆蓋，且隨覆蓋度增大而增大。研究區(qū)的NDVI以0、0.3、0.5 為界限分成四類。

2.2 評價因子檢驗與篩選

為保證各評價因子的獨立性，采用皮爾遜相關系數(shù)法對8 個評價因子進行檢驗與篩選，計算得到的相關性系數(shù)見表2。由表2 可知各因子相關性系數(shù)較低，絕對值均小于0.5?？紤]計算準確性，選取全部8 個評價因子進行模型的建立與訓練。

表2 評價因子相關性檢驗Table 2 Correlation test of evaluation factors

2.3 最大熵模型

最大熵模型是由最大熵原理推導實現(xiàn)的。它的中心思想是，在只知道關于未知分布的部分約束條件時，應選取滿足這些約束條件且熵值最大的概率分布。最大熵原理通過熵的最大化來表示其可能性。因此最大熵模型就是在所有可能的概率模型中，熵最大的模型就是最好的模型[13]。

2.4 模型計算

首先將169 個滑坡數(shù)據(jù)和8 個評價因子加載到MaxEnt 軟件V3.4.1 中。然后將災害點按照7∶3 的比例隨機分配為訓練集（118 個）與測試集（51 個），具體分布情況見圖2。其中訓練集用來訓練數(shù)據(jù)和建立模型，測試集用于驗證模型的可靠性。設置迭代次數(shù)500 次，隨機計算10 次，結果取均值。將結果以ASCII 文件輸出，導入ArcGIS 經過重分類后得到滑坡易發(fā)性的分區(qū)，具體分布見圖4。

3 結果及分析

3.1 滑坡易發(fā)性分區(qū)

結果參考正態(tài)分布理論與經驗對滑坡易發(fā)區(qū)進行分類[23]，共分為五個易發(fā)區(qū)。P<0.08 為極低易發(fā)區(qū)，0.08≤P<0.25 為低易發(fā)區(qū)、0.25≤P<0.53 為中易發(fā)區(qū)、0.53≤P<0.72 為高易發(fā)區(qū)、P≥0.72 為極高易發(fā)區(qū)。得到的結果見圖4。經統(tǒng)計，分區(qū)占比依次為11.48%、41.28%、25.21%、10.87%、11.16%。

圖4 滑坡易發(fā)性分區(qū)Fig.4 Landslide susceptibility zoning

其中極高易發(fā)區(qū)主要集中在道路兩側及山谷地帶，路東路西分布平均，高度普遍不高。說明人為工程擾動已經嚴重的影響了斜坡的穩(wěn)定性，促使了滑坡的發(fā)生。此外有三處極高易發(fā)區(qū)集中地帶（圖4）。結合各個評價因子圖層的分布可知，其中1 處巖體較軟弱且斷層密度較高。2 處的巖體較為破碎，河流發(fā)育且有一處斷裂穿過。3 處在研究區(qū)南部，該區(qū)域的地震加速度較高，災害極為發(fā)育。中、高易發(fā)區(qū)以極高易發(fā)區(qū)為中心向外輻射，基本覆蓋了道路兩側山體分水嶺以下的部分。低易發(fā)區(qū)和極低易發(fā)區(qū)的輻射范圍更遠，主要集中在一些河谷平原地區(qū)。

3.2 模型準確性評價

滑坡易發(fā)性評價結果是否準確決定了模型的適用性和可靠性，因此有必要對結果精度進行驗證。受試者工作特征曲線（ROC）常用來檢測模型的預測精度，ROC 曲線下面積AUC值在0～1,曲線越靠近左上角AUC值越大,模型預測精度越高,一般認為AUC值>0.7 時,MaxEnt 模型的預測結果便為可信[24]。

從預測精度ROC 曲線可知運行10 次的結果中訓練集的AUC最大值為0.864，最小值為0.839。結果偏差不大，且屬于較高的數(shù)值?？梢宰C明模型具有優(yōu)秀的空間預測能力和精度，且計算結果十分穩(wěn)定（圖5）。

圖5 ROC 曲線和AUC 值箱形圖Fig.5 ROC curve and AUCvalue box diagram

3.3 評價因子與滑坡的關系

評價因子的貢獻率檢驗結果如圖6所示。坡度、高程、植被覆蓋指數(shù)、坡向是影響滑坡發(fā)生的主導因素，4 項累計貢獻率達77.3%。

圖6 評價因子對滑坡發(fā)生的貢獻率Fig.6 Contribution rate of evaluation factors to landslide

通過分析響應曲線(圖7)來判斷滑坡與評價因子之間的關系。一般認為當存在概率大于0.5 時，其對應的評價因子對滑坡發(fā)生產生影響[25]?；碌拇嬖诟怕孰S坡度的增加逐漸增加，隨著坡度的增加，坡面附近的應力卸荷的范圍擴大，坡腳處應力集中程度加大，因此滑坡的發(fā)生幾率也隨之增加[26?28]。當坡度大于20°時最容易發(fā)生滑坡災害。在本研究區(qū)的模擬結果中，高程的敏感范圍在4 445 m 以下，產生這樣的結果主要是由于災害多發(fā)生在沿公路區(qū)域，斜坡由于公路切坡，更容易產生滑坡危險，因此公路兩側的邊坡要重點防護和清理。當NDVI指數(shù)小于0.23 時，滑坡表現(xiàn)為敏感。NDVI越小表示斜坡上的裸露比例越大，坡面失去植物根系的保護，在雨水的沖刷下更容易受到破壞，因此較為容易發(fā)生滑坡。對于坡向來說，敏感區(qū)間在90°～315°。結果表明陽坡比陰坡更能促進滑坡的形成，與現(xiàn)場實際情況相符。

圖7 評價因子響應圖Fig.7 Response curve of evaluation factors

4 結論

（1）以現(xiàn)場調查及遙感解譯的169 個滑坡災害點作為樣本數(shù)據(jù)，選取8 個評價因子，建立最大熵模型進行滑坡易發(fā)性分區(qū)。研究區(qū)滑坡極高和高易發(fā)區(qū)面積占總面積的11.6%和10.87%。主要集中在公路兩側。同時在斷層線密集處、巖性軟弱破碎處、地震動加速度峰值高處，均有災害點集中的現(xiàn)象。

（2）結合模型預測結果的響應曲線可以得到控制滑坡發(fā)生的變量。最主要的四項因素分別為坡度、高程、植被覆蓋與坡向。其中當坡度大于20°、高程靠近路面、NDVI小于0.23、坡向在90°～315°之間時最容易發(fā)生滑坡。

（3）經過10 次運算，得到AUC最大值為0.864，最小值為0.839。說明最大熵模型適用于研究區(qū)的滑坡易發(fā)性研究，且結果準確可靠?？梢詾橹心峤煌ɡ鹊赖貐^(qū)的防災減災工作提供預警。