0.   引言

近年来，基于“一带一路”发展战略需求，对于南亚通道的建设开始逐渐完善。作为中国与尼泊尔密切往来的口岸，吉隆是中尼铁路建设的必经之地。这里地质构造复杂、内外动力作用强烈，具有强地震风险、强地壳变形、强应力集中等特征，是重大地质灾害频发区域^[1]。经沿线实地调查，该地区滑坡灾害发生较频繁，且滑坡具有突发性和不确定性。因此掌握滑坡的易发性对于滑坡预测和防灾减灾具有重要意义。

滑坡易发性评价指某一区域内现有自然条件下已经发生了和诱发因素下滑坡灾害与空间分布的定性或定量关系^[2-3]。目前多种滑坡灾害易发性的评价方法中，大体可以总结为定性和定量两类方法。早期基于专家经验的现场分析和评价因子图层叠加分析的方法过分依赖评价者的主观意识，可比性较低。随着地理信息系统技术、卫星遥感技术和数字地形数据的全面推进，机器学习模型借助GIS平台逐步被引用到地质灾害评价领域^[4-6]。

MaxEnt 模型是以最大熵原理为理论，以JAVA语言为支撑的机器学习模型。在1957年由杰恩斯首先提出最大熵原理，它的中心思想是，在只知道关于未知分布的部分约束条件时，应选取满足这些约束条件且熵值最大的概率分布。Phillips等在2006年提出了最大熵模型，专门用于生态建模和物质分布预测。该模型在2013年被引入滑坡灾害的预测研究中。相比其他概率模型，如频率比、证据权法、多标准决策分析和逻辑回归等，最大熵模型既可以考虑滑坡与评价因子之间的关系，同时还能充分考虑所有因子本身之间的关系。相较于层次分析法，最大熵模型不需要专家经验，在易发性分析时不需要加入主观假设。比较同为机器学习模型的神经网络，最大熵模型可以克服人工神经网络的“黑箱”限制。总之，最大熵模型虽然数学推导较复杂，但模型适应性、稳定性和灵活性更高，预测准确度也更高^[7-16]。

文中使用最大熵模型对研究区滑坡易发性建模，然后通过ROC曲线验证模型的效果。并进一步对各评价因子的贡献度和敏感阈值进行研究。为中尼交通廊道区域的滑坡灾害预警提供参考。

1.   研究区概况

吉隆县地处西藏自治区日喀则地区西南部，南面和西南面与尼泊尔相邻，东临聂拉木县，北面与萨噶县搭界，自古就是西藏通往南亚的交通要道^[17]（图1）。

图  1  研究区地理位置

Figure  1.  Location of the study area

下载: 全尺寸图片 幻灯片

研究区地处青藏高原西南部，全县平均海拔在4000 m以上，地势北高南低。按地貌成因可划分为侵蚀剥蚀山地地貌和河谷侵蚀堆积地貌。位于冈底斯构造带和北喜马拉雅构造带之间，新构造运动活跃，属缓慢抬升区^[18-19]。加之地震多发且强烈，诱发差异性抬升等特殊的第四纪地质环境，使之成为崩塌、滑坡、泥石流多发区。研究区出露地层主要为第四系松散堆积层与前震旦纪江东岩组变质岩和前震旦纪中酸性侵入岩。

研究区调查的G216国道是吉隆县通往吉隆镇的唯一公路，总长86.61 km。范围是以道路为中心向两侧扩展5 km区域。经过对研究区的现场细致排查及室内遥感解译补充，共确定169处滑坡，如图2（a）所示。调查沿途多为高山峡谷路段，地层多为风化破碎的抗侵蚀能力较弱的基岩，重力侵蚀活跃，切坡造成的大临空面、高陡坡脚的斜坡体。沿线性工程两侧地表松散物质较多，坡脚均存在小范围垮塌的可能。较为典型的滑坡如图2（b）（c）。图2（b）滑坡位于G216国道K3470处，滑坡体主体为土石混合体。人工开挖后，坡脚被削，形成牵引式滑坡。坡脚处修筑有45 m的浆砌片石挡墙。在暴雨、强风作用下，坡面上块石易发生滚落。图2（c）本拉错滑坡，位于吉隆镇北侧。坡脚紧靠交通线，剪出口冲过公路与护栏，已严重影响道路通行，目前斜坡处于破坏休止阶段，但在地震、暴雨等因素情况下，会加速边坡的破坏，影响道路安全。

图  2  研究区灾害点分布及典型灾害

Figure  2.  Distribution of disaster points and typical disasters in the study area

下载: 全尺寸图片 幻灯片

2.   研究方法

2.1   评价因子选取

滑坡的发生受到多种因素的影响。文中在充分考虑了研究区尺度和特性后，选取了高程、坡向、坡度、断层密度、河流密度、岩性、地震峰值加速度、植被覆盖指数8个评价因子（表1）。利用ArcGIS对评价因子进行栅格化及重分类处理，生成图3所示8个图层。各图层栅格列数行数保持一致分别为2547、5594，栅格总数6575436个。

图  3  研究区评价因子图层

Figure  3.  Layer of evaluation factors in the study area

下载: 全尺寸图片 幻灯片

表  1  评价因子选取及其来源

Table  1.  Selection of evaluation factors and sources

评价因子	数据来源
高程	地理空间数据云（http://www.gscloud.cn）
坡向	ArcGIS 高程提取
坡度	ArcGIS 高程提取
断层密度	西藏自治区吉隆县地质图（全国地质资料馆http://www.ngac.org.cn/）及现场调查
河流密度	西藏自治区吉隆县地质图（全国地质资料馆http://www.ngac.org.cn/）及现场调查
岩性	西藏自治区吉隆县地质图（全国地质资料馆http://www.ngac.org.cn/）及现场调查
地震加速度	中国地震动参数区划图 (http://www.gb18306.net/)
植被覆盖指数	地理空间数据云 （http://www.gscloud.cn）

下载: 导出CSV 

| 显示表格

（1）高程可以直接反映研究区的地形起伏^[20]。将地理空间数据云DEM数据，通过ArcGIS进行裁剪、拼接等处理，生成研究区DEM图；（2）坡向对斜坡的影响表现为向阳坡受日照时间较长，相应湿度较低，岩石表面风化程度较高，易形成松散堆积物，从而斜坡的危险性变高。通过ArcGIS将DEM图层转化得到坡向图层；（3）坡度大小影响了滑坡体的移动速度和规模程度。通过ArcGIS将DEM图层转化得到坡度图层；（4）断裂构造使断层带及其附近一定范围内的岩土体结构遭到破坏，降低了坡体的完整程度。断裂带越密集岩体也就越破碎，因此选用断层带的线密度作为一个评价因子。以现场调查为基础地质资料为补充，利用ArcGIS线密度计算工具计算得到；（5）研究区地表水系交错，根据地质图利用ArcGIS线密度计算工具生成河流线密度作为滑坡易发性评价因子；（6）地层岩性决定着斜坡岩土体的强度，控制着地质灾害发育和分布规律。根据西藏地区区域性地质图及现场调查将研究区地层岩性划分为四组，分别为：较坚硬变质岩岩组、坚硬岩浆岩岩组、坚硬沉积岩岩组、较软弱沉积岩岩组；（7）地震峰值加速度（PGA）的大小表征地球内动力的强弱^[21]。吉隆县受地震的影响较大，该地区PGA为0.1～0.20 g；（8）植被覆盖指数用归一化差分植被指数（NDVI）来定量表征。通常与植物蒸腾作用、光合作用有关，是指示植被生长状态及植被覆盖度的最佳指示因子^[22]。将Landsat8的数据1～7波段加载到ArcGIS中进行波段合成。然后计算NDVI ，计算公式为：

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

R——红波段处的反射率值。

NDVI取值[1，−1]，负值表示地面覆盖为云、水、雪等；0表示有岩石或裸土等；正值表示有植被覆盖，且随覆盖度增大而增大。研究区的NDVI以0、0.3、0.5为界限分成四类。

2.2   评价因子检验与筛选

为保证各评价因子的独立性，采用皮尔逊相关系数法对8个评价因子进行检验与筛选，计算得到的相关性系数见表2。由表2可知各因子相关性系数较低，绝对值均小于0.5。考虑计算准确性，选取全部8个评价因子进行模型的建立与训练。

表  2  评价因子相关性检验

Table  2.  Correlation test of evaluation factors

	高程	坡向	坡度	断层密度	河流密度	岩性	地震峰值加速度	植被覆盖指数
高程	1
坡向	0.056	1
坡度	−0.307	0.388	1
断层密度	0.414	−0.177	−0.104	1
河流密度	0.251	−0.141	−0.124	0.030	1
岩性	0.459	−0.125	−0.098	0.398	0.200	1
地震峰值加速度	0.491	0.322	−0.441	0.446	−0.182	0.442	1
植被覆盖指数	0.406	-0.466	0.468	−0.156	−0.107	−0.163	−0.467	1

下载: 导出CSV 

| 显示表格

2.3   最大熵模型

最大熵模型是由最大熵原理推导实现的。它的中心思想是，在只知道关于未知分布的部分约束条件时，应选取满足这些约束条件且熵值最大的概率分布。最大熵原理通过熵的最大化来表示其可能性。因此最大熵模型就是在所有可能的概率模型中，熵最大的模型就是最好的模型^[13]。

2.4   模型计算

首先将169个滑坡数据和8个评价因子加载到MaxEnt软件V3.4.1中。然后将灾害点按照7∶3的比例随机分配为训练集（118个）与测试集（51个），具体分布情况见图2。其中训练集用来训练数据和建立模型，测试集用于验证模型的可靠性。设置迭代次数500次，随机计算10次，结果取均值。将结果以ASCII文件输出，导入ArcGIS经过重分类后得到滑坡易发性的分区，具体分布见图4。

图  4  滑坡易发性分区

Figure  4.  Landslide susceptibility zoning

下载: 全尺寸图片 幻灯片

3.   结果及分析

3.1   滑坡易发性分区

结果参考正态分布理论与经验对滑坡易发区进行分类^[23]，共分为五个易发区。P<0.08为极低易发区，0.08≤P<0.25为低易发区、0.25≤P<0.53为中易发区、0.53≤P<0.72为高易发区、P≥0.72为极高易发区。得到的结果见图4。经统计，分区占比依次为11.48%、41.28%、25.21%、10.87%、11.16%。

其中极高易发区主要集中在道路两侧及山谷地带，路东路西分布平均，高度普遍不高。说明人为工程扰动已经严重的影响了斜坡的稳定性，促使了滑坡的发生。此外有三处极高易发区集中地带（图4）。结合各个评价因子图层的分布可知，其中1处岩体较软弱且断层密度较高。2处的岩体较为破碎，河流发育且有一处断裂穿过。3处在研究区南部，该区域的地震加速度较高，灾害极为发育。中、高易发区以极高易发区为中心向外辐射，基本覆盖了道路两侧山体分水岭以下的部分。低易发区和极低易发区的辐射范围更远，主要集中在一些河谷平原地区。

3.2   模型准确性评价

滑坡易发性评价结果是否准确决定了模型的适用性和可靠性，因此有必要对结果精度进行验证。受试者工作特征曲线（ROC）常用来检测模型的预测精度，ROC曲线下面积AUC值在0～1, 曲线越靠近左上角AUC值越大, 模型预测精度越高, 一般认为AUC值>0.7时, MaxEnt模型的预测结果便为可信^[24]。

从预测精度ROC曲线可知运行10次的结果中训练集的AUC最大值为0.864，最小值为0.839。结果偏差不大，且属于较高的数值。可以证明模型具有优秀的空间预测能力和精度，且计算结果十分稳定（图5）。

图  5  ROC曲线和AUC值箱形图

Figure  5.  ROC curve and AUC value box diagram

下载: 全尺寸图片 幻灯片

3.3   评价因子与滑坡的关系

评价因子的贡献率检验结果如图6所示。坡度、高程、植被覆盖指数、坡向是影响滑坡发生的主导因素，4项累计贡献率达77.3%。

图  6  评价因子对滑坡发生的贡献率

Figure  6.  Contribution rate of evaluation factors to landslide

下载: 全尺寸图片 幻灯片

通过分析响应曲线(图7)来判断滑坡与评价因子之间的关系。一般认为当存在概率大于0.5时，其对应的评价因子对滑坡发生产生影响^[25]。滑坡的存在概率随坡度的增加逐渐增加，随着坡度的增加，坡面附近的应力卸荷的范围扩大，坡脚处应力集中程度加大，因此滑坡的发生几率也随之增加^[26-28]。当坡度大于20°时最容易发生滑坡灾害。在本研究区的模拟结果中，高程的敏感范围在4 445 m以下，产生这样的结果主要是由于灾害多发生在沿公路区域，斜坡由于公路切坡，更容易产生滑坡危险，因此公路两侧的边坡要重点防护和清理。当NDVI指数小于0.23时，滑坡表现为敏感。NDVI越小表示斜坡上的裸露比例越大，坡面失去植物根系的保护，在雨水的冲刷下更容易受到破坏，因此较为容易发生滑坡。对于坡向来说，敏感区间在90°～315°。结果表明阳坡比阴坡更能促进滑坡的形成，与现场实际情况相符。

图  7  评价因子响应图

Figure  7.  Response curve of evaluation factors

下载: 全尺寸图片 幻灯片

4.   结论

（1）以现场调查及遥感解译的169个滑坡灾害点作为样本数据，选取8个评价因子，建立最大熵模型进行滑坡易发性分区。研究区滑坡极高和高易发区面积占总面积的11.6%和10.87%。主要集中在公路两侧。同时在断层线密集处、岩性软弱破碎处、地震动加速度峰值高处，均有灾害点集中的现象。

（2）结合模型预测结果的响应曲线可以得到控制滑坡发生的变量。最主要的四项因素分别为坡度、高程、植被覆盖与坡向。其中当坡度大于20°、高程靠近路面、NDVI小于0.23、坡向在90°～315°之间时最容易发生滑坡。

（3）经过10次运算，得到AUC最大值为0.864，最小值为0.839。说明最大熵模型适用于研究区的滑坡易发性研究，且结果准确可靠。可以为中尼交通廊道地区的防灾减灾工作提供预警。