1. 引言
近年来，全球自然灾害频繁发生，而我国由于其复杂而脆弱的地质环境，已成为世界上受地质灾害影响最为严重的国家之一。滑坡灾害作为自然界最常见的灾害之一，表现出破坏性大，危害性强的特点。据《2022年中国自然资源统计公报》显示，我国全年发生地质灾害5659起，其中滑坡3919起，占比达到69% [1]。大别山地区位于中国安徽省、河南省和湖北省的交界处[2]，大别山地区地形复杂，气候多变，滑坡灾害频发。近年来，随着全球气候变化和人类活动的加剧，大别山地区的滑坡灾害呈现出日益严重的趋势。
随着科技的飞速发展，地质灾害评估已经从基于定性的分析向基于定量的研究转变，并且已经取得了显著的成就。常用的数据模型有：信息量模型[3]、逻辑回归模型[4]、人工神经网络模型[5]、支持向量机模型[6]、随机森林模型[7]等。但是各种模型都有其优势和局限，例如，支持向量机(SVM)模型的基本原理是最小化结构风险，这比传统方法中的经验风险最小化更具优势，因为它是通过解决一个有约束的二次优化问题来保证模型的表现。不过，在小样本情况下，样本的噪声或矛盾信息可能会对SVM的预测效果产生明显的不利影响。因此，在特定的灾害易发环境中，无法找到一个模型能够普遍适用于各种情况。
国内外很多学者在滑坡易发性评价方面进行了卓有成效的研究，并取得众多显著的成果。其中，经验模型[8]、信息量模型[9] [10]、统计预测模型[11] [12]以及机器学习模型[13] [14]在滑坡易发性评估中得到了广泛应用。21世纪初，Saro Lee [15]在2003年开发了一种基于人工神经网络的方法来分析滑坡易发性，通过建立包含地形、土壤类型、植被、地质和土地覆盖等因素的空间数据库，以此训练神经网络模型来预测滑坡的位置。Achu [16]使用了概率均值、概率中位数、概率加权均值和委员会平均值等集合模型，基于较多独立MLT的平均值，得出基于MLTs的微观层次区划可能会提高滑坡危险度图的效率。国内在该领域研究起步较晚，武雪玲等[17]基于地形、地质和遥感影像等多源数据，提取滑坡易发性评价因子，构建卷积神经网络模型，定量预测滑坡易发性。而在最近几年，刘帅等人[18]基于优化随机森林模型，对西秦岭极端降雨事件诱发的群发滑坡易发性进行了评价研究；李泽芝等人[19]则针对秦巴山区堆积层滑坡易发性不同单元进行了性能评价对比研究；王本栋等人[20]采用支持向量机(SVM)、BP神经网络和随机森林(RF)三种典型机器学习算法进行滑坡易发性评价；林琴等人[21]则基于基尼系数的加权随机森林、XGBoost和LightGBM算法在滑坡易发性中的性能进行了研究。
滑坡的发生不仅与其所处的地理位置有关，还与地质环境存在紧密的关联。目前的滑坡易发性模型大多以单个像素作为评估单位，这种方法往往忽略了目标像素与其周围地质环境的相互作用，从而影响了制图的精度[22]。深度学习，特别是以卷积神经网络为代表的方法，因其强大的非线性预测能力而广泛应用于图像分类[23] [24]、迁移学习[25]等领域，能有效提取图像数据中的深层信息。此外，深度学习模型通常使用矩阵图像[26]-[28]作为基础数据，相比于传统机器学习模型，这一方式能更好地考虑到目标样本周边像素的影响，从而提高模型的综合表现。
基于此，本文将收集关于安徽省大别山区的地质、气象和地形数据，利用随机森林、支持向量机、神经网络等多种机器学习算法，对大别山区的滑坡易发性进行综合评估。通过对比分析不同算法的评价结果，探讨各种算法在滑坡易发性评价中的优缺点，比较三个模型的评价效果。此外，通过对不同区域的滑坡易发性进行对比分析，揭示大别山区滑坡灾害的主要影响因素和潜在危险区域。该研究为安徽省大别山区的滑坡防治工作提供科学依据。
2. 研究区域与研究数据
2.1. 研究区域
Figure 1. Location of the study area
图1. 研究区位图
大别山地区的山体构造复杂，是淮阳山字型构造体系的核心部分，并且构成了秦岭褶皱带的延伸，这一地区的山体主要呈西北–东南走向，而东段则转为东北–西南走向。整个山脉长达270千米，一般海拔在500至800米之间，山地的主要部分海拔约为1500米，山脊的海拔通常在1200至1600米之间。窗体顶端大别山区的地形由断层活动切割形成众多菱形断块，在东南侧的黄梅至桐城一带，山脚线挺直，山坡陡峭，坡度超过50˚，形成了显著的断层崖。山区两侧发育了丰富的水系，北侧主要河流如灌河、史河、潢河、竹竿河、洗河等流入淮河，南侧的如靳河、大悟河、淆水、摄水等则流入长江。山间谷地宽广开阔，并有河漫滩和阶地平原，是主要农耕地区。山地多深谷陡坡，地形复杂，坡向多变，坡度多在25˚~50˚。研究区域如图1所示。
2.2. 数据来源
本文用到的数据主要有中国多年度地市行政区划边界数据[29]，滑坡泥石流灾害点数据，DEM数字高程数据，中国城市基本土地利用分类数据，全国降雨量数据等。具体来源见表1。
Table 1. Research data and their sources
表1. 数据来源
数据名称 |
数据来源 |
历史滑坡 |
水工环地质信息服务平台 |
DEM |
地理空间数据云平台 |
行政区划 |
阿里云数据可视化平台 |
多年降雨 |
国家青藏高原科学数据中心 |
土地利用 |
中国科学院资源环境科学数据中心 |
卫星影像 |
地理空间数据云平台 |
3. 研究方法
3.1. 环境因子的选取
影响滑坡发生的因素主要有地形地貌、降雨量、植被覆盖及人类工程活动等，滑坡是多种因素联合作用的结果[29]。现有的滑坡易发性评价过程环境因子选取尚未确定统一标准，滑坡的发生与地形因素发生密切相关，数字高程模型(digital elevation model, DEM)通常被作为提取地形数据的基础数据源[30]。本文相关的地形因素分析涉及高程、坡度、坡向、剖面曲率、地形湿度等5个因子，都是基于DEM数据进行提取而获得。并且本文对NDVI、多年平均降水等因子进行分析，以增强研究成果的合理性和准确性。
3.2. 滑坡易发性模型构建
在滑坡影响因子的权重分析的基础上，本文基于机器学习方法构建滑坡易发性区划模型，通过逻辑回归、人工神经网络和随机森林三种机器学习方法的模型进行计算，对研究区域进行滑坡易发性区划并制图。通过随机采样的方法，选择与滑坡点等量的非滑坡点，这些非滑坡点均匀分布于整个研究区内，将筛选出的影响因子的栅格属性值赋予了这些采样数据(滑坡和非滑坡点)，这些影响因子属性被用作模型的输入变量，其中滑坡(标记为1)和非滑坡(标记为0)作为模型的输出变量。然后将这些带有属性的数据进行划分，其中70%的数据用作训练数据集，剩余的30%则用作测试数据集，将这些影响因子图层作为输入数据，通过训练得到的滑坡易发性模型，生成反映滑坡易发性影响因子权重分布的图。
3.2.1. 样本数据的选择与处理
本研究包含滑坡样本与非滑坡样本，其中滑坡点样本的数量为616，为了构建所需的二分类模型，本文采用1:1的比例进行样本选择，即从非滑坡区域随机抽取616个样本作为非滑坡样本数据，这样的抽样策略有助于确保模型训练的公平性和数据的平衡。
3.2.2. 模型精度评价方法
为了准确评估二分类问题的模型性能，通常采用ROC曲线和AUC值这种评价指标。在这种分类问题中，实例分为两类：正类(如“滑坡”)和负类(如“非滑坡”)。对于分类的结果，存在四种可能的情况：若实际为“滑坡”且预测正确，称为真正例(TP)；若实际为“非滑坡”却错误预测为“滑坡”，称为假正例(FP)；若实际为“非滑坡”且预测正确，则为真负例(TN)；若实际为“滑坡”却被预测为“非滑坡”，称为假负例(FN)。
ROC曲线，全称为受试者工作特性曲线，是一种评估分类模型性能的工具，该曲线通过设置不同的概率阈值，生成一系列的分类决策，然后将这些决策的结果与实际情况进行对比。曲线越靠近左上角，表明模型的性能越优秀，该曲线最靠近左上角的点表示模型的最佳阈值，此时模型的预测误差最小，而AUC值，即ROC曲线下的面积，是用来衡量模型整体预测准确度的指标。
3.2.3. 基于逻辑回归的滑坡易发性区划模型
逻辑回归是一种统计学上的分析方法，广泛用于分类问题。尽管其名称中含有“回归”二字，逻辑回归实际上是用来解决二分类问题的。其基础来自于概率论和统计学，它可以将变量的线性组合映射为一个预测值，这个值表示的是一个事件发生的概率，该模型解决了线性回归输出值可能会超出0和1范围的问题。对数几率转换是通过一个逻辑函数(或称作sigmoid函数)来实现的，逻辑函数的数学表达式为：
(1)
这里的x是输入变量的线性组合，比如：
，其中，
是模型参数，
是特征。
3.2.4. 基于随机森林的滑坡易发性区划模型
随机森林是一种集成学习方法，该模型通过结合多个决策树的预测结果来改善模型的精确度和稳定性的技术，其核心思想是集体智慧的概念。随机森林由许多个决策树组成，每个树都是一个分类器(对于分类问题)或预测器(对于回归问题)。单个决策树可能对数据过度拟合，但当多个树集成时，整体模型的泛化能力通常会增强。训练每个决策树时，随机森林算法从原始数据集中使用自助采样(有放回的抽样)的方式选出一个大小相等的新数据集作为训练集。这意味着一些观测可能会被多次选中，而一些则可能不被选中。在每个决策点，算法会从所有特征中随机选出一部分特征，然后仅使用这些特征来决定最佳的分割方式。这种随机特征选择减少了模型的方差，提升了模型在未知数据上的表现力。
3.2.5. 基于XGBoost的滑坡易发性区划模型
XGBoost (eXtreme Gradient Boosting)是一个优化的分布式梯度提升库，设计用来高效地实现梯度提升算法。该模型属于梯度提升框架，它在梯度提升的基础上进行了系统的优化和加速。梯度提升是一个通过加模型(additive model)的方法逐步在当前的基础上添加新的预测器来改进模型的过程。每个新的预测器(通常是决策树)都根据之前所有预测器预测的残差(实际值与预测值之间的差距)进行训练。与传统的梯度提升决策树(GBDT)相比，XGBoost在目标函数中引入了正则化项。这两个正则化项是用来控制模型的复杂度的，它们分别对应树的叶子节点上的分数的L2范数和树的结构(例如叶子节点的数量)。通过并行化树学习过程中的某些步骤来加速计算。虽然树的构建本身是顺序的，不能并行化，但是对于树中每个节点的分割点的计算可以并行化。XGBoost采用的是深度优先的策略，而不是传统梯度提升使用的广度优先。其在树的分支过程中进行剪枝，当一个分支的分裂导致目标函数的增益(包括正则化项)小于一个阈值时，就停止分裂。
4. 研究结果
4.1. 滑坡影响因子分析结果
4.1.1. 高程
由图2可知，中低高程区域的滑坡密度较高，在100~250米的高程区间的滑坡密度最高，达到了0.080个/km2。这表明在这个高程范围内的区域可能因为特定的地形、土壤和水文条件，导致滑坡风险增加。随着高程的增加，滑坡密度逐渐减少，从250~405米和405~580米的高程区间，滑坡密度分别为0.050个/km2和0.040个/km2，显示出随着高程增加，滑坡密度逐渐减少的趋势。当高程达到580~780米时，滑坡密度进一步降低到0.020个/km2。在低高程区域的滑坡密度最低，在高程小于100米的区域，滑坡密度仅为0.002个/km2，这可能是因为这些低地区域相对平坦，土壤稳定性较高，或者是因为水文条件对滑坡的促进作用较小。这些发现表明，高程是影响滑坡密度的重要因素之一。在中低高程区域，可能由于地形较为陡峭，加之特定的土壤和水文条件，导致滑坡风险增加。而在较高的高程以及低平地区，滑坡密度则相对较低。
Figure 2. Statistical diagram of landslide density at elevation and within elevation range
图2. 高程和高程范围内滑坡密度统计图
4.1.2. 坡度
根据图3所示，低坡度区域的滑坡密度较低，在坡度小于2.9˚的区域，滑坡密度仅为0.003个/km2，这表明在较为平缓的地形中，滑坡事件较少发生。这可能是因为在低坡度的条件下，重力引起的下滑动力相对较小。中等坡度区域的滑坡密度较高，当坡度在2.9˚到17.4˚之间时，滑坡密度显著增加，尤其是在7.8˚到12.7˚的坡度区间，滑坡密度达到最高，为0.080个/km2。这可能是因为这个坡度范围内，土壤、岩石等地表物质在重力作用下更易发生位移。较高坡度区域的滑坡密度有所减少，当坡度超过17.4˚时，滑坡密度开始减少，回落到0.050个/km2。这可能是因为在非常陡峭的坡度条件下，虽然下滑动力增大，但是由于物质较少或者已经在之前的事件中发生滑移，可动物质变少，导致新的滑坡事件相对减少。这些发现表明，坡度是影响滑坡密度的一个重要因素，特别是中等坡度区间内的滑坡风险较高。在坡度2.9˚到28.1˚范围内，滑坡发生的概率较高，这些地区的稳定性不如低坡度地区，斜坡稳定性较差。
Figure 3. Slope distribution and landslide density statistics of each slope
图3. 坡度分布和各坡度的滑坡密度统计图
4.1.3. 坡向
Figure 4. Slope aspect distribution and landslide density statistics of each slope aspect
图4. 坡向分布和各坡向的滑坡密度统计图
北向、西向和西北向坡面的滑坡密度相对较低，这三个方向的坡面，滑坡密度为0.02个/km2，这是所有已列出方向中最低的(图4)。这可能与坡面受到的日照量相对较少，从而导致土壤湿度保持在较高水平，植被覆盖相对较密集有关。较高的土壤湿度和良好的植被覆盖能够增加土壤的稳定性，减少滑坡的发生。其他方向的坡面滑坡密度则较为一致，东北、东、东南和南向的坡面，滑坡密度均为0.03个/km2，表明这些方向上的坡面相对于北向坡面，滑坡风险稍高。这可能与这些方向坡面受到的日照更多，土壤干燥，植被覆盖较少有关，从而减少了土壤的稳定性。
4.1.4. 剖面曲率
图5中剖面曲率小于2.9的范围的滑坡密度为0.005个单位/km2。对于2.9~7.8的曲率范围，密度增加到每平方公里0.04个单位。曲率范围为7.8~17.4时，密度保持在每平方公里0.06个单位。对于17.4~22.3的曲率范围，密度略增到每平方公里0.18个单位。而22.3~28.1的曲率范围，密度持续保持在每平方公里0.06个单位。这些数据说明大部分滑坡集中在剖面曲率较高的区域。剖面曲率较高通常意味着该区域的坡度较大。坡度是影响滑坡发生的重要因素之一，坡度越大，重力作用下的向下力分量越大，材料(如岩石、土壤等)向下滑动的可能性也就越大。其次，剖面曲率较高的区域通常意味着这些区域的地形较为陡峭，雨水和地表水更容易快速流失，不易渗透地表。然而，在某些情况下，这种快速流动的水可以侵蚀坡脚或增加坡面负荷，导致滑坡发生。此外，水分的快速流动也可能导致地下水位的局部上升，增加滑坡的风险。在植被的覆盖方面，由于坡度大，植被的根系很难深入地下，因此土壤的固结作用减弱，容易在重力和其他外力作用下发生移动。此外，这些区域在雨季可能会遇到快速的土壤湿度增加，进一步降低土壤的稳定性，增加滑坡的风险。
Figure 5. Statistical diagram of terrain profile curvature and landslide density of each profile curvature
图5. 地形剖面曲率和各剖面曲率的滑坡密度统计图
4.1.5. 地形湿度指数
在图6中，低湿度地形的滑坡密度较高，当地形湿度小于5.9时，滑坡密度最高，达到0.05个/km2。这可能是因为低湿度条件下，土壤更加干燥，抗剪强度降低，更易发生滑坡。随着地形湿度的增加，滑坡密度并不是简单的线性减少，而是显示出一种复杂的非线性关系。在湿度为5.9~7.5时，滑坡密度降低至0.03个/km2，而在湿度进一步增加至7.5~9.0时，滑坡密度降至最低，为0.01个/km2。然后，当湿度增加到9.0~11.0和11.0~13.2时，滑坡密度略有上升，稳定在0.02个/km2。这种关系提示，地形湿度对滑坡密度有重要影响，但这种影响受多种因素共同作用的影响，包括土壤类型、植被覆盖状况、降雨量等。在低湿度条件下，土壤更容易干裂，减少了土壤的凝聚力，从而可能增加了滑坡的风险。而在湿度适中的区间，可能因为土壤湿润度适宜和植被的稳定作用，滑坡密度降低。但是，当湿度过高时，过多的水分可能导致土壤饱和，增加了滑动面的形成概率，因此滑坡密度有所上升。
Figure 6. Distribution of terrain humidity index and statistical diagram of landslide density in different terrain humidity
图6. 地形湿度指数分布和各地形湿度的滑坡密度统计图
4.1.6. 植被覆盖率
Figure 7. NDVI distribution and landslide density statistics for each NDVI
图7. NDVI分布和各NDVI的滑坡密度统计图
从图7可知，低NDVI值的区域滑坡密度较低，当NDVI值小于0.5时，滑坡密度为0，表明这些区域的植被覆盖度较低，可能是裸露土地或者植被稀疏的地区。这些区域虽然植被覆盖度低，但滑坡密度并不高，可能是因为这些地区的地形、土壤类型等其他因素不利于滑坡的发生。随着NDVI值的增加，滑坡密度开始出现并增加，特别是NDVI值在0.5~0.6时，滑坡密度为0.001个/km2；而当NDVI值在0.6~0.7时，滑坡密度显著增加到0.007个/km2。这表明在植被覆盖度较高的区域，滑坡密度也较高。这可能是因为植被既可以通过其根系加固土壤，降低滑坡风险，但在某些情况下，植被的存在也可能因为雨水截留和土壤湿润度增加，反而增加了滑坡的可能性。这些发现表明，NDVI作为一种反映地表植被覆盖度的指标，在分析滑坡风险时，提供了一种重要的视角。尽管植被覆盖通常被认为可以减少侵蚀和滑坡的风险，但在某些条件下，过密的植被可能因为增加了土壤的湿度，反而增加滑坡的风险。
4.1.7. 多年平均降水量
从图8可以看出，低降雨量区域的滑坡密度非常低，当年平均降雨量小于1020 mm时，滑坡密度仅为0.0002个/km2，这可能是因为降雨不足以渗透深层土壤或造成地表径流，从而减少了滑坡的风险。随着降雨量的增加，滑坡密度显著增加，在降雨量为1020~1140 mm的区间，滑坡密度增加到0.0008个/km2；而当降雨量进一步增加到1140~1270 mm时，滑坡密度急剧上升到0.0100个/km2。这种趋势在1270~1390 mm的降雨区间达到顶峰，滑坡密度达到0.0300个/km2。在到达降雨量阈值之后，滑坡密度有所下降，当年平均降雨量在1390~1490 mm时，滑坡密度减少到0.0200个/km2。这可能是由于在极端的降雨条件下，虽然滑坡的概率增加，但是过量的降雨也可能导致潜在的滑坡区域已经发生滑坡，从而在一定程度上减少了新的滑坡事件的发生。以上分析表明，多年平均降雨量与滑坡密度之间存在明显的相关性。在降雨量适中的区域，随着降雨量的增加，滑坡密度增加，反映了降雨是滑坡发生的重要触发因素。然而，在降雨量非常高的区域，可能因为地表和土壤已经达到饱和状态，再增加的降雨对滑坡密度的影响有所减弱。
Figure 8. Distribution of multi-year average rainfall and statistical graph of landslide density by rainfall period
图8. 多年平均降雨量分布和降雨时段滑坡密度统计图
4.1.8. 土地利用类型
从图9可知，林地的滑坡密度最高：林地的滑坡密度为0.040个/km2，这是所有土地利用类型中最高的。这可能与林地区域的地形、土壤类型、植被覆盖和根系发展有关，虽然植被通常可以稳定土壤，减少侵蚀，但在某些情况下，如极端降雨或不适宜的土壤条件下，林地区域仍然可能面临较高的滑坡风险。耕地和草地的滑坡密度中等：耕地的滑坡密度为0.010个/km2，而草地为0.020个/km2。这表明这些土地利用类型的滑坡风险处于中等水平。耕地可能因土壤被翻动和作物种植的方式影响了土壤的稳定性；草地可能因植被覆盖和根系深度的不同而具有不同的滑坡风险。水域和建设用地的滑坡密度最低：水域和建设用地的滑坡密度分别为0.004个/km2和0.003个/km2，这是所有类型中最低的。这可能是因为这些区域的土地利用方式导致地面更加稳固，或者是因为这些区域通常不位于易发生滑坡的地形上。
Figure 9. Statistical graph of land use types and landslide density of each land use type
图9. 土地利用类型和各土地利用类型的滑坡密度统计图
4.2. 模型评价结果与精度分析
三大模型评价结果
通过比较三种机器学习方法的ROC曲线和AUC值(图10)，可以得出随机森林模型的精度最高(AUC值为0.8462)，其次是XGBoost模型(AUC值为0.8443)，而逻辑回归模型的精度较低(AUC值为0.7400)。这一结果表明随机森林模型在处理复杂的地质数据和非线性关系时表现更优越。由于在随机森林模型中，每棵树的构建都是在随机选取的数据子集和特征上进行，这种方法能有效提高模型在多样化数据上的稳定性和准确性。XGBoost模型作为一种高效的梯度提升框架，也显示出了很高的预测精度，其优化的计算速度和对过拟合有很好的控制，尤其适用于有大量数据的情形。相比之下，逻辑回归模型虽然在数学上较为简单并易于理解，但在处理复杂和高维数据时可能不够灵活，因此在本研究中显示出较低的AUC值，逻辑回归在特征关系较为线性时效果较好，但对于地形、地质等多因素复杂交互的滑坡易发性评价，其性能有限。
图11是基于最优机器学习算法随机森林模型计算得到的大别山滑坡易发性区划图。在诸多的影响因子中，高程的权重系数0.24是最高的，这也是山区的易发性远高于周边地区的主要原因。另外多年平均降水的权重系数也达到0.13，大别山区南坡的降雨量高于北部山区，降雨量越大，其洪涝易发性越强，因此大别山南坡是易发性是最高的。
Figure 10. ROC curves of logistic regression model, random forest model and XGBoost model
图10. 逻辑回归模型、随机森林模型和XGBoost模型的ROC曲线
Figure 11. Landslide susceptibility zoning map based on random forest model
图11. 基于随机森林模型的滑坡易发性区划图
5. 结论与展望
本文通过综合运用机器学习技术，对安徽省大别山区滑坡易发性进行了深入分析。采用随机森林、逻辑回归、XGBoost等多种机器学习算法，构建了滑坡易发性评估模型，有效地评价了区域滑坡的潜在风险。结果表明，大别山区植被覆盖NDVI在0.7~0.8之间其滑坡点分布最密集，且在林地滑坡点分布最多。随机森林是滑坡易发性拟合最优的模型。在诸多因子分析中，高程、坡度、剖面曲率和多年平均降水对滑坡易发性有显著影响。大别山区的滑坡易发性高于周边地区，且南坡的滑坡易发性风险最大。