基于最大熵模型的中尼交通廊道滑坡易發(fā)性分析

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

（1.中國(guó)科學(xué)院地質(zhì)與地球物理研究所/中國(guó)科學(xué)院頁(yè)巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室,北京 100029；2.中國(guó)科學(xué)院地球科學(xué)研究院,北京 100029；3.中國(guó)科學(xué)院大學(xué),北京 100049）

0 引言

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

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

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

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

1 研究區(qū)概況

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

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

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

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

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

2 研究方法

2.1 評(píng)價(jià)因子選取

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

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

表1 評(píng)價(jià)因子選取及其來源Table 1 Selection of evaluation factors and sources

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

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

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

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

2.2 評(píng)價(jià)因子檢驗(yàn)與篩選

為保證各評(píng)價(jià)因子的獨(dú)立性，采用皮爾遜相關(guān)系數(shù)法對(duì)8 個(gè)評(píng)價(jià)因子進(jìn)行檢驗(yàn)與篩選，計(jì)算得到的相關(guān)性系數(shù)見表2。由表2 可知各因子相關(guān)性系數(shù)較低，絕對(duì)值均小于0.5?？紤]計(jì)算準(zhǔn)確性，選取全部8 個(gè)評(píng)價(jià)因子進(jìn)行模型的建立與訓(xùn)練。

表2 評(píng)價(jià)因子相關(guān)性檢驗(yàn)Table 2 Correlation test of evaluation factors

2.3 最大熵模型

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

2.4 模型計(jì)算

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

3 結(jié)果及分析

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

結(jié)果參考正態(tài)分布理論與經(jīng)驗(yàn)對(duì)滑坡易發(fā)區(qū)進(jìn)行分類[23]，共分為五個(gè)易發(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ū)。得到的結(jié)果見圖4。經(jīng)統(tǒng)計(jì)，分區(qū)占比依次為11.48%、41.28%、25.21%、10.87%、11.16%。

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

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

3.2 模型準(zhǔn)確性評(píng)價(jià)

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

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

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

3.3 評(píng)價(jià)因子與滑坡的關(guān)系

評(píng)價(jià)因子的貢獻(xiàn)率檢驗(yàn)結(jié)果如圖6所示。坡度、高程、植被覆蓋指數(shù)、坡向是影響滑坡發(fā)生的主導(dǎo)因素，4 項(xiàng)累計(jì)貢獻(xiàn)率達(dá)77.3%。

圖6 評(píng)價(jià)因子對(duì)滑坡發(fā)生的貢獻(xiàn)率Fig.6 Contribution rate of evaluation factors to landslide

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

圖7 評(píng)價(jià)因子響應(yīng)圖Fig.7 Response curve of evaluation factors

4 結(jié)論

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

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

（3）經(jīng)過10 次運(yùn)算，得到AUC最大值為0.864，最小值為0.839。說明最大熵模型適用于研究區(qū)的滑坡易發(fā)性研究，且結(jié)果準(zhǔn)確可靠。可以為中尼交通廊道地區(qū)的防災(zāi)減災(zāi)工作提供預(yù)警。