基于机器学习均值化的地质灾害易发性评价

邵大江，叶 辉，王金亮**，周京春，角媛梅，沙晋明

(1.云南师范大学 地理学部，云南 昆明 650500；2.云南省高校资源与环境遥感重点实验室，云南 昆明 650500；3.云南省地理空间信息工程技术研究中心，云南 昆明 650500；4.云南大学 国际河流与生态安全研究院，云南 昆明 650500；5.福建师范大学 地理科学学院，福建 福州 350007)

地质灾害是在自然因素、人文因素或者自然与人为因素共同作用下形成的对生态环境或人类生命财产造成破坏的地质作用或现象[1].常见的如滑坡、泥石流等地质灾害具有突发性、不可控性、破坏性的特征.频发的地质灾害给人类生命财产安全造成了极大地威胁，开展地质灾害易发性评价，可以基于评价结果针对不同的易发性等级提出相应的防灾减灾对策[2].

地质灾害易发性评价是以地质环境条件为基础，参考地质灾害现状的静态因素，能够预测一定区域内发生地质灾害的可能性[3].地质灾害评价模型，可以笼统地分为非确定性模型以及确定性模型[4].更进一步则可分为经验模型、力学模型、因子模型[5]，其中的因子模型是目前地质灾害评估的主流方法，是众多学者的研究对象，如信息量法[6]、地理加权回归模型等.当机器学习在地质灾害预测中发展时，按新的分类法大致可分为以下5 类：机器学习模型、数理统计模型、确定性模型、启发式模型、基于滑坡编录的概率模型[7].机器学习在地质灾害的预测中，存在处理大量数据且对数据相关性和环境因素等的关联性没有过多要求的特点.随着获取数据的便利、计算能力的提升、模型评价算法的日趋完善，机器学习在地质灾害易发性的评价上更加广泛，如梯度提升树、人工神经网络、决策树、随机森林等[8].刘艳辉等[8]利用多种机器学习算法构建了青川县的滑坡灾害预警模型，选取2018 年6 月26 日青川县日常预警业务进行实例校验，结果表明当日的17 处灾害点全部落入预警区；屠水云等[9]基于确定性系数模型（Certain Factors,CF）与基于CF 的逻辑回归模型（Logistic Regression,LR）对沿河县进行评价，其评价结果表明CF-LR 耦合模型AUC 提高了0.096，具有更高的精度；李坤等[10]基于随机森林模型（Random Forest,RF）和支持向量机模型（Support Vector Machine,SVM）对东川区泥石流易发性进行评价，发现2 种机器学习算法结合流域单元建立的山区泥石流易发性评价模型均具有很好的准确性及稳定性，其中RF 优于SVM 模型；Jin 等[11]利用逻辑回归模型从时空角度对西昌市火灾后泥石流的易发性进行预测，预测结果表明随着时间变化高易发区逐渐减少，人工播种加速了土壤的自我修复；Pal 等[12]基于贝叶斯和树模型的集成框架评价泥石流和滑坡灾害，研究表明坡度的变量重要度最高，是泥石流、滑坡的主要诱因，并且RF 模型精度最高；Qiu 等[13]利用一种混合机器学习模型，绘制气候变化下吉隆藏布江流域冰川相关泥石流易发性，通过研究发现研究区61%的面积具有较高易发性，同时确定性因子-遗传算法（Genetic Algorithm,GA）-支持向量机（CF-GA-SVM）模型具有相对较高的鲁棒性.

由于单一模型的局限性以及地质灾害的复杂性，多模型结合的预测方法可以更好地对研究区地质灾害易发性评价做出较准确的预测.本文对南华县地质灾害易发性提出相应预防对策，对高易发区实施重点防控，对中低易发区加强群测群防监控.

1 研究区与方法

1.1 研究区概况南华县，隶属于云南省楚雄彝族自治州，位于100°43′27″～101°22′06″E，24°43′34″～25°22′24″N，其东邻牟定、楚雄市，南接景东、楚雄市，西连弥渡，北毗姚安、祥云，是四川、云南的交通要道，享有“九府通衢”的美称（图1）.同时，南华县内道路河网交错，其中320 国道、217 省道、楚大高速公路和广大铁路穿境而过[14].南华县地处云南“山”字构造的脊柱部分以及滇中高原西部，境内最高点海拔2 861 m，最低点海拔963 m，地势东南低、西北高，东北丘陵连绵，西南群山迭起，东部和中部地形起伏较为平缓，地形复杂，山河沟壑相间陈列.地质构造复杂，断层交错，其中以洱海—红河断裂带居大.全县年均降水量约843.2 mm，年均温约14.8 ℃.全县总面积约为2 343 km2，山区面积居多，约占96%，森林覆盖率高达63.73%，主要的土地利用类型有耕地、灌木林、林地、水域、城镇用地及其他用地.4 条主要河流分别为龙川江、礼社江、兔街河、马龙河.因地质环境脆弱，地质灾害频发，据南华县地质灾害详细调查成果显示，截止2015 年，南华县共发育有地质灾害点188 个，单位面积灾害密度远高于其他县（市）.其中崩塌灾害点4 个、滑坡灾害点155 个、泥石流21 个、潜在不稳定斜坡8 个，其中滑坡占84.2%，为主要地质灾害类型.

图1 研究区地理位置及灾害分布Fig.1 Geographical location and disaster distribution of the study area

因此，本文选取道路、水系、断层、岩体类型、高程、降雨、土地利用、坡度、坡向、土壤类型、植被覆盖度共11 个影响因子作为自变量，采取梯度提升树算法，主要包括LightGBM（Light Gradient Boosting Machine）、XGBoost（eXtreme Gradient Boosting）、CatBoost（Categorical Boosting）3 种算法，以及基于3S 的信息量模型与地理加权回归模型对其进行预测评价，将机器学习输出结果投射为发生概率，取3 种算法均值为其最终结果，并将信息量模型与地理加权模型输出值归一化，引入接收者操作特征曲线（Receiver Operating Characteristic Curve,ROC）与AUC（Area Under ROC Curve）值检验各模型预测结果精度.通过各模型的对比以寻求精度上的提高，以期为提高区域地质灾害预报预警精度提供参考.

1.2 研究方法

1.2.1 梯度提升树算法 梯度提升树算法是由Friedman 等提出的迭代决策树算法的一种.其基本思想是由多棵决策树共同决策，最终把所有的树结论累加得到最终的答案.其中损失函数采用平方差来表示，每一棵树通过学习之前的树的结论和残差拟合得到当前的树，其中残差值为真实值与预测值之差，整个迭代过程所得到的累加树即为提升树[15].梯度提升树算法可以分为回归和分类，回归与分类在思想上并没有太大的区别，在研究中采用了分类算法中的二分类，如公式（1）：

式中：FM(x)为强学习器，F0(x)为初始分类器，M为分类树总数量，Jm为叶子节点个数，Cm,˙j为最佳叶子节点拟合值，Rm，j为第m棵回归树下的叶子结点.

1.2.1.1 基于梯度提升树的算法实现 LightGBM、XGBoost、CatBoost 作为基于梯度提升树算法的3种代表性实现，在各自的算法创新上都具有各自的优势，也曾在各领域获得过较高的荣誉.

（1）与传统的GBDT 模型不同，XGBoost 模型采用泰勒二阶展开对损失函数进行优化，支持CPU 的多线程并行运算，并在损失函数中加入了正则项，使其运算效率和泛化能力相比于其他机器学习算法大大提升[16].该算法具有可分布式处理、精度高、可适应特征维度高、可移植性、可有效防止过拟合等优点[17].

（2）LightGBM 算法的诞生是为了弥补XGBoost的缺陷，解决处理海量数据时存在的问题.其特点是将XGBoost 的level-wise 分裂策略替代为leafwise 分裂策略，使用基于直方图的决策树算法.具有优化计算速度与内存的使用、稀疏优化、准确率优化、网络通讯优化等优势，并且支持3 种模式并行：特征并行、数据并行、投票并行[18].

（3）CatBoost 算法是集成学习算法中一种较为新颖的算法，该算法具有独特的对称树结构，在通过计算叶子结点的值构造决策树的过程中，对特征进行了量化的度量[19].自动采用特殊的方式处理类型特征；考虑到特征之间的联系，采用了组合类别特征，极大地丰富了特征维度；CatBoost 算法的基模型采用了对称树，并且对于叶子结点值的计算方式和传统的Boosting 算法也有区别，对其进行优化，以防止过拟合.

1.2.1.2 参数调节及评价系数 在Sacikit-learn的梯度提升树类库中，梯度提升树算法分为回归类和分类类，参数基本相同，大致可分为两类：一类为Boosting 框架参数；另一类是弱学习器参数即CART 回归树参数.重要参数及取值见表1.

表1 重要参数取值Tab.1 Values of key parameters

判断机器学习模型好坏的评价标准，常用AUC 值衡量.接受者操作特征曲线（ROC）又称为感受性曲线，广泛应用于各种模型精度的检验，通过ROC 曲线可以衡量模型对正负类样本的预测能力及可分性，可以作为衡量一个模型好坏的标准.但是ROC 不能说明哪个分类器的效果更好，而AUC 值定义为ROC 曲线去坐标轴围成的线下面积，该值介于0.5～1.0，AUC 值越高则表明该模型精度越高，预测能力越强[20].

在真实值和预测值的关系中，TP（True Positive）为正样本分为正例，FN（False Negative）为正样本分为负例，TN（True Negative）为负样本分为负例，FP（False Positive）为负样本分为正例.设置不同阈值，计算不同阈值的真假正类率TPR和FPR，然后以TPR和FPR分别为纵轴与横轴，生成的曲线即为ROC 曲线[21]，TPR与FPR计算公式如下：

1.2.2 均值运算法 虽然梯度提升树算法的3种实现有其独特的优势，但依然存在不同程度的缺陷.XGBoost 由于每一轮的迭代都需要反复遍历读取整个训练数据集，会消耗非常多的时间；其次，预排序方法对时间以及空间的消耗都十分巨大；CatBoost 在没有包含分类变量时，效果是最差的；LightGBM 容易导致过拟合.在分类的研究过程中，得出Catboost 对正样本的预测效果最差，但是对负样本的预测效果最好，而XGBoost 对正样本的预测效果最好，对负样本的预测效果最差.由于地质灾害是一种多因素多影响力所引起的效应，单一的模型及算法难以达到完美的预测结果，因此，对于当前的预测模型都有其所不足，应结合各模型相对应的优点进行预测，结合各模型相对应的优点进行预测，以达到预测精度的提升.本文将机器学习的3 种算法预测结果转化为概率，求其均值作为该研究区的地质灾害的发生概率，以达到降低单一模型算法不足的问题.

2 数据及处理

2.1 数据源基于地质灾害统计数据的可获取性，选择2015 年南华县地质灾害详查数据，并结合道路、水系、断层、岩体类型分布、高程、降雨、土地利用、坡度、坡向、土壤类型和植被覆盖度等数据对地质灾害易发性进行评价研究.行政区划、水系、道路分布图等矢量数据来源于全国地理信息资源目录服务系统（https://www.webmap.cn/）.灾害点分布图（1∶5 万）、断层分布图（1∶20 万）、岩体类型分布图（1∶20 万）、土地利用数据（30 m）、NDVI数据（1 km）来源于中国科学院地理科学与资源研究所（https://www.igsnrr.ac.cn/）.高程数据（30 m）从地理空间数据云（https://www.gscloud.cn/#page1）获取，2015 年年降水量数据（1 km）从国家地球系统科学数据中心（http:www.geodata.cn/）获取，土壤类型数据（1 km）从世界土壤数据库（http:www.ncdc.ac.cn/portal/metadata/a948627d-4671-4568-b1b6-fe02e302af09）获取.坡度、坡向数据由高程数据计算得到，植被覆盖度数据由NDVI 计算得到.

2.2 影响因子量化及处理本文基于数据可获取性及前人研究基础，共选取11 个影响因子（图2）.由于数据类型包括连续型、离散型、以及文本类型，需要对数据进行量化处理，便于接下来的模型学习计算.

图2 南华县地质灾害的影响因子Fig.2 Impact factors of geological disasters in Nanhua County

如表2 所示，根据土地利用现状分类准则，将南华县土地利用类型分为6 个类型；根据诱发地质灾害的特点将降水量分为3 个区段，56.30%的地质灾害点发生在900～960 mm 区段；根据南华县岩体类型的物理力学属性，将其分为3 种类型[22]；根据美国农业部所制作的土壤质地三角图，将土壤类型分为3 种类型；根据地质灾害在不同高程的分布情况将高程设置为4 个等级，其中89.30%的地质灾害发生点发生在1 600～2 400 m 区段；利用像元二分法通过归一化植被指数提取NDVI 植被覆盖度数据，并根据45%、60%、75%等3 个阈值将植被覆盖度划分为4 种覆盖类型；利用ArcGIS 软件处理DEM 获取坡度及坡向数据，坡向数据分为9个类型；根据地质灾害的分布特点将坡度划分为4个区段；利用缓冲区分析对道路、水系、断层3 个影响因子制作多环缓冲区，距离设为500、1 000 m和1 500 m[23].为方便模型计算，对数据进行清洗，并将11 个已分类的影响因子的不同类映射为0～8 之间的代码.对于位于裁剪边缘而造成影响因子数据缺失的地质灾害点，在模型预测研究中作删除处理，后续分析阶段加入用做验证.

表2 评价因子代码表Tab.2 Table of evaluation factor codes

3 分析与讨论

3.1 训练结果对南华县188 个地质灾害点进行机器学习数据集制作，根据前人研究按1∶2 的比例随机选取非灾害点，通过数据预处理剔除数据缺失点，将数据集按照3∶7 的比例分为验证集与训练集，训练结果（表3）显示XGBoost、LightGBM、CatBoost等3 种梯度提升树算法准确率分别为0.828、0.822、0.804，AUC 值分别为0.797、0.789、0.760.由表3可知，3 种算法均具有较高的准确率，且3 种模型算法均具有较好的精度，适用于南华县的地质灾害易发性评估，其中XGBoost 在准确性和模型精度上都略高于LightGBM 和CatBoost，但LIghtGBM运算时长远小于XGBoost 和CatBoost.

表3 机器学习训练集结果Tab.3 Results of machine learning training set

3.2 地质灾害易发性预测结果本文将整理好的数据分别代入3 种机器学习模型进行训练，将训练好的模型保存并利用ArcGIS 软件连接研究区内的各个栅格数据点进行地质灾害评估.同时，带入信息量模型以及地理加权回归模型进行易发性评估，利用自然间断法对预测结果划分为易发区和非易发区（图3）.

图3 南华县地质灾害易发性预测结果Fig.3 Prediction results of geological hazard susceptibility in Nanhua County

各种模型算法所计算得到的易发性评价专题图有些许差异，但对于易发区的预测大致相同.总体上，地质灾害易发区几乎遍布南华县全境，呈现2 种分布情况，即4 条由西北—东南走向的易发带与南华县西南部的易发区.统计得到信息量模型预测易发区面积最大，为1 169.07 km2，占总面积的55.5%；地理加权回归模型预测易发区面积最小，为801.97 km2，占总面积的37.3%.根据地质灾害详查数据，得到各乡镇已发生的地质灾害数量，如表4 所示地质灾害多发乡（镇）主要有一街乡、马街镇、兔街镇，以上乡（镇）均位于预测结果的易发区.

表4 南华县各乡（镇）地质灾害数量Tab.4 Number of geological disasters in each township(town) of Nanhua County

基于易发区对预测结果进行二次划分，地质灾害分为非易发区、低易发区、中等易发区、高易发区4 个等级（图4）.通过二次划分，可以更加突出机器学习模型与信息量模型在南华县中部和北部的4 条地质灾害易发带以及该县西南部的1 个易发区.“四带一区”主要由中等易发区和高易发区构成，结合所在地理位置，由东北往西南4 条易发带分别位于龙川江与杭瑞高速公路相交处、羊草河与南大线相交处、一街河与王湛庄至红土坡公路相交处、永景线与兔街河相交处，西南易发区则位于弥河线、永景线、龙街河相交处，并且“四带一区”所处地理位置均有断层带、其中“一区”所处地理位置有5 条断层带.根据图4（d）所示，相比于机器学习与信息量模型，地理加权回归模型对南华县东北部的预测结果并没有出现易发带.其原因是详查数据中，该区域地质灾害发生点较少，地灾样本点的空间分布相对稀疏，而地理加权回归模型更倾向于反映样本点之间的空间关系，其结果则相反：东北部大部分地区属于非易发区，西南部等地属于地质灾害高易发区.根据图4（f）所示，信息量模型预测结果相对于机器学习与地理加权回归模型预测结果，在东北部出现较多由中等易发区和低易发区组成的易发带，其原因是在西南部地质灾害多发区存在多条道路与河流，在信息值的计算过程中河流与道路被赋予较大的信息熵，从而导致其他影响因子所占权重降低，因此在东北部形成围绕道路与河流的易发带.

图4 南华县地质灾害易发性分级预测结果Fig.4 Classification prediction results of geological hazard susceptibility in Nanhua County

根据详查地质灾害数据统计得到不同易发区灾害频数统计表（表5）与不同易发区面积及灾害密度统计图（图5）.在不考虑自然间断法的影响下（表5），地理加权回归模型高易发区已发生地质灾害数量最多，有118 个，占南华县总灾害数量62.7%.而信息量模型高易发区最少，有56 个，占比29.8%.除信息量模型外，其他模型高易发区的灾害分布均达到45%以上，其中均值法模型以及地理加权模型预测结果达到了50%以上.在非易发区中，地理加权回归模型已发生地质灾害数量最少，有28 个，占比14.9%.而CatBoost 模型非易发区地质灾害数量最多，有41 个，占比21.8%.如图5 所示，在不同易发区的面积变化上，地理加权回归模型高易发区与非易发区面积均高于其他模型，而信息量模型高易发区与非易发区面积则均低于其他模型.机器学习算法所预测的各易发区面积无明显差异.不同模型各易发区的灾害密度变化在非易发区与低易发区无明显差异，但在中等易发区与高易发区变化明显.

表5 南华县易发程度分级及灾害频数表Tab.5 Classification of susceptibility and disaster frequency in Nanhua County

图5 南华县不同灾害易发性等级面积及密度统计图Fig.5 Statistical map of area and density of different disaster susceptibility grades in Nanhua County

3.3 卡方检验卡方检验是一种用途很广的计数资料的假设检验方法[24].该方法主要是比较多个样本及分类变量的相关性，其基本思想是提出一个无效假设，即期望频数与观察频数没有差别.其计算结果称为统计量，表示观察值与理论值的偏离程度，统计量越小则不拒绝无效假设，表明实际情况与期望值并无较大差别，统计量越大则表明比较资料之间具有显著差异性.通过卡方检验求得各影响因子统计量，统计量大的影响因子与地质灾害的发生具有较强的相关性，即各因子对灾害发生的贡献率（表6）.通过统计量可知，道路、断层、降水量、水系的数值远高于其他影响因子，说明上述因子与地质灾害发生的相关性极大.

表6 影响因子统计量Tab.6 Statistics of influencing factors

3.4 核密度分析通过核密度（Kernel Density Estimation,KDE）分析可以得到研究对象的密度变化，客观准确地表达其空间分布状况.通过卡方检验得到道路、断层、水系对地质灾害的发生具有较高的相关性，对上述3 个影响因子进行线状要素的核密度分析，通过分析结果以达到对卡方检验结果的侧面支撑.本文将核密度估算结果划分为低密度、中等密度、高密度3 个区段，得到核密度分布图（图6）.根据分布结果可知，南华县已发生的地质灾害中，48 个位于低密度地区、96 个位于中密度地区，44 个位于高密度地区.图6 结果显示，74.5%的地质灾害发生点位于中高密度地区，其结果与卡方检验具有一致性.

图6 南华县核密度分析图Fig.6 Thematic map of kernel density estimation in Nanhua County

3.5 精度分析在剔除位于边缘影响因子缺失的样本后，对南华县177 个地质灾害样本与375 个非地质灾害样本的预测中，Catboost 对地质灾害样本的预测效果最差，但是对非地质灾害样本的预测效果最好；XGBoost 对地质灾害样本的预测效果最好，对非地质灾害样本的预测效果仅高于信息量模型.基于3 种梯度算法的均值所预测的效果取得了较好的提升，对正负样本的预测都得到了较好的增强效果，其预测结果见表7.

表7 分类预测结果Tab.7 Classification prediction results

ROC 曲线[25]（图7）可以作为衡量一个模型好坏的标准，其线下面积定义为AUC 值，可以直观地衡量模型精度.在南华县地质灾害易发性评价的研究中，各模型AUC 值从高至低依次是均值法、地理加权回归模型、XGBoost、LightGBM、Catboost、信息量模型，分别为0.933 7、0.917 4、0.917 1、0.915 4、0.899 9、0.898 1，其中均值法精度较地理加权回归模型、XGBoost、LightGBM、Catboost、信息量模型分别提高了1.7%、1.8%、2.0%、3.8%、4.0%.

图7 ROC 曲线图Fig.7 ROC curve

3.6 讨论（1）在南华县的地质灾害易发性预测研究中，基于机器学习获取海量数据的原则，未考虑各因子之间的相关性.6 种模型预测结果表明南华县地质灾害高易发区具有沿河性、沿路性、沿断层性的特征.基于这一特征，在研究中通过卡方检验对各影响因子进行统计量计算，结果显示道路、断层、降水量、水系的统计量远高于其他影响因子，即相对于其他因子，上述因子对灾害发生的贡献率较大.同时，在此基础上对道路、河流、断层做核密度分析.结果表明74.5%的地质灾害发生在中高密度地区，与卡方检验结果具有一致性.

（2）通过不同模型预测结果对比分析，信息量模型在预测过程中具有负拟合现象.其易发区面积最大，但是高易发区却面积小，因此高易发区地质灾害数量远低于其他模型，并低于中等易发区，不符合客观事实.同时信息量模型忽略了空间差异性，计算过程容易对单一因子赋予高熵值，导致高易发区较少且集中.在南华县的研究中对研究区整体的河流、道路及断层赋予了较高的信息熵值，导致东北部出现“假性易发带”.

地理加权回归模型在研究中倾向于反映样本点的空间关系，导致过拟合现象的发生.其易发区面积最小，但是高易发区面积却远高于其他模型预测结果，其易发区主要集中在高易发区，所以导致低易发区与中等易发区面积较少，这也是高易发区地灾数量过多的原因.非易发区面积也明显高于其他模型预测结果，这一现象说明地理加权回归模型的过拟合现象容易造成预测结果两极化.

基于梯度提升树算法的3 种模型以及均值法预测结果并未出现明显异常现象.其预测算法基于各影响因子的综合特征，进行多次迭代计算，不单独赋予其中任何一个影响因子较高的权重，从而较少出现单一因子权重过高的现象，同时没有受到各样本点空间关系的影响，因此各易发区面积及灾害数量分布合理.

（3）在以梯度提升树为主的南华县地质灾害易发性对比分析中，均值法、地理加权回归模型、XGBoost、LightGBM、CatBoost、信息量模型等6种模型对于研究区的易发性预测具有相似性，高易发区主要分布在南华县南部、中部及北部部分地区，其中一街乡、兔街镇、马街镇为地质灾害多发乡（镇）.但是，在精确程度上均值法的AUC 值达到最高，为0.933 7，精度较其他5 种模型分别提高了1.7%、1.8%、2.0%、3.8%、4.0%.关于地质灾害的正负样本预测，Catboost 对地质灾害样本的预测效果最差，但对非地质灾害样本的预测效果最好，而XGBoost 对地质灾害样本的预测效果最好，均值法则对正负样本的预测效果都取得了较大的提升.

4 结论

（1）南华县地质灾害高易发区主要位于南华县南部、中部及北部部分地区，其空间分布呈现“四带一区”现象，即4 条易发带与1 个易发区.其中一街乡、兔街镇、马街镇为地质灾害多发乡（镇）.

（2）6 种模型对于研究区的易发性预测具有相似性，但在精确程度上，均值法的AUC 值达到最高（0.933 7），精度较其他5 种模型分别提升了1.7%、1.8%、2.0%、3.8%、4.0%.地理加权回归模型则容易出现过拟合现象，信息量模型则容易造成单一因子权重过高现象.在地质灾害正负样本的预测中，Catboost 对地质灾害样本的预测效果最差，但对非地质灾害样本的预测效果最好；XGBoost 对地质灾害样本的预测效果最好；均值法则对正负样本预测精度都取得了较大的提升.

（3）通过卡方检验以及核密度分析可知，道路、断层、降水量、水系是南华县地质灾害发生的主要影响因子，其地质灾害具有沿河性、沿道路性、沿断层性的特征.