0.   引言

滑坡是中国分布最多、威胁最大的地质灾害并且具有分布范围广、发生概率高以及危害程度深等特点^[1]，是人民生命安全和财产安全的巨大威胁。根据2022年中华人民共和国自然资源部发布的中国自然资源统计公报^[2]显示：2 022年共发生地质灾害5 659起，共造成90人死亡、16人失踪、34人受伤，直接经济损失15亿元。地质灾害对区域内居民的生命和财产安全有着很大的威胁，对地质灾害如滑坡进行易发性评价以及易发程度分区，可以供政府进行政策参考，制定相关监测或巡防方案，优化防灾减灾工作，减少地质灾害发生所引起的各类损失^[3]。

目前，滑坡易发性研究主要分为定性研究和定量研究两类方法^[4]。其中，定性研究方法主要包括专家打分法、层次分析法等^[5]，定性方法依赖主观分析和实践经验，主观性过强。定量方法主要以评价因子和滑坡发生之间的关系为基础^[6]，在初始的研究中一般使用统计学模型，应用较广泛的统计模型包括信息量模型^[7]、逻辑回归模型^[8]等，如陈立华等^[9]将GIS与信息量模型相结合开展滑坡易发性评价并取得较优的预测结果。然而，由于统计模型结构简单，多类非线性数据的输入和滑坡的不同特性都会导致模型的预测精度降低^[10]。对比传统的统计方法模型，机器学习对统计函数的依赖性低，机器学习的机制来源于输入数据，通过训练数据构建模型规则，从而反映滑坡的发生可能性与影响因子之间的关系。机器学习对于非线性数据的处理能力强，处理速度快，能够有效处理特征之间的关联性^[11]，如神经网络模型^[12]、支持向量机^[13]、随机森林^[14]等模型在滑坡易发性评价领域有着广泛的应用，如李成林等^[15]使用神经网络模型对研究区滑坡易发性进行评价，又如黄发明等^[16]学者将支持向量机模型与聚类分析相结合并构建SOM-SVM模型，在非滑坡栅格单元选取方面取得了较好的效果。近年来，随着深度学习的快速发展，更多的模型被使用在滑坡易发性领域，如卷积神经网络^[17]、长短时记忆网络^[18]等新模型，这些模型在滑坡多维特征提取、模型强训练等方面有着优秀的表现。

然而机器学习模型的精度常常依赖于样本的数据量及数据质量和对模型的参数设置，选择合适的参数对机器学习模型至关重要，对滑坡易发性评价的精度有着很大的影响。对于该问题前人的研究成果主要有两种解决方案，一方面可以通过优化特征选择提升滑坡的预测精度^[19]，常用的特征优化方法如遗传算法等，如高轩宇等^[20]建立Bi-GRU模型对滑坡影响因子进行特征重要度排序并选择；另一方面可以通过优化算法去确定模型的超参数，如白鲨优化^[21]、粒子群优化^[22]等进化、群智能算法，智能优化算法在避免模型陷入局部最优解的困境方面有着很好的效果，使用算法优化机器学习模型后能提高全局优化速度从而快速搜索最优超参数解,提高模型的精确度和运行速度，如杨灿等^[23]将贝叶斯优化算法与四种常见的机器学习模型相结合用于探究贝叶斯优化算法的可行性，陈芯宇等^[24]使用麻雀搜索算法对两类机器学习模型进行优化并检验其优化效果。现在的优化算法与机器学习模型的结合已经很成熟，但是不同算法对模型优化效果探究方面研究较少，随着优化算法的不断更新，选择合适的算法对模型进行优化显得尤为重要，因此应当试验不同算法对常见机器学习模型的优化效果并加以检验。另一方面，在对模型效果检验方面，一般多使用统计学方法对模型性能进行评估，因此拟采用SHAP分析，提高模型的可解释性和可信度。基于以上理论，本文旨在：（1）选择3种常见优化算法（贝叶斯优化算法、麻雀搜索算法、大猩猩优化算法），构建BP神经网络模型并将3种优化算法与BP模型相结合，采用AUC值、F1分数等参数分别对BP单模型和三类优化模型进行精度检验，使用K折验证的方法评估模型性能；（2）为不同优化算法提供相同模型的对比平台，增加不同算法之间优化性能对比的可信度，为后续研究提供理论支撑；（3）利用SHAP（shapley additive explanation）分析的方法构建可解释的机器学习结合滑坡易发性评价的模型，绘制秭归地区的滑坡易发性分区图，为预防秭归地区滑坡地质灾害提供理论支持。

1.   研究区域及数据来源

1.1   研究区域

本研究选取湖北省宜昌市秭归县作为研究区域。研究区地处湖北省西部，三峡工程坝上库首，东与宜昌市夷陵区相交界，南与长阳土家族自治县交界，西临巴东县，北接兴山县^[25]。秭归地势西南高东北低，东段为黄陵背斜，西段为秭归向斜，属长江三峡山地地貌^[26]。研究区由长江及其支流、沟溪共同构成其水文网，流域面积达1952.5 km²，占县域总面积的80.4%，香溪河是研究区中最大的一条长江支流，根据赋存条件的不同，将该研究区地下水分为松散岩类孔隙水、碎屑岩类孔隙裂隙水和碳酸盐岩类岩溶水3类。共统计研究区内滑坡532处，研究区概况及滑坡点分布如图1所示。

图  1  研究区位置示意图

Figure  1.  Overview map of the study area

下载: 全尺寸图片 幻灯片

1.2   研究数据来源

研究数据包括地形地貌数据、水文数据、地层数据及卫星遥感数据，其主要数据来源及数据类型见表1所示，其中滑坡灾害点数据来自湖北省地质局水文地质工程地质大队提供的1∶5万秭归县地质灾害调查风险普查数据；高程、坡度、坡向以及水系数据来源于中国科学院计算机网络信息中心地理空间数据云平台 （http://www.gscloud.cn）；工程地质岩组数据来源于湖北省地质局水文地质工程地质大队提供的秭归县1∶5万区域地质图；归一化植被指数（NDVI）数据来源于国家科技基础条件平台—国家地球系统科学数据中心（http://www.geodata.cn）；道路数据来源于OpenStreetMap （https://www.openstreetmap.org）；平面曲率与剖面曲率数据由坡度、坡向数据提取得到。

表  1  数据来源及数据类型

Table  1.  Data sources and types

数据名称	数据来源	数据类型
滑坡灾害点	1∶5万秭归县地质灾害调查 风险普查数据	点数据
高程、坡度、 坡向、水系	地理空间数据云DEM数据	30 m × 30 m栅格数据
工程地质岩组	秭归县1∶5万区域地质图	1∶50000栅格数据
归一化植被指数	Landsat8 OLI	30 m × 30 m栅格数据
距道路距离	OpenStreetMap	矢量数据
平面曲率、剖面曲率	坡度、坡向数据	30 m × 30 m栅格数据

下载: 导出CSV 

| 显示表格

2.   研究方法

2.1   评价因子选取

滑坡的发生取决于多类因素，目前的各项研究对滑坡的评价因子选取尚未达成统一，本文主要研究优化算法改进模型的模型性能，因此借鉴前人经验及其他研究已有的研究成果，采取最为常见的9个评价因子：高程、坡度、坡向、平面曲率、剖面曲率、归一化植被指数、地层岩性、距道路距离以及距水系距离进行研究^[27]。使用ArcGIS软件提取各因子并进行重分类如图2所示：

图  2  评价因子图层

Figure  2.  Evaluation factor layers

下载: 全尺寸图片 幻灯片

2.2   评价因子特征选择

特征选择（feature selection）是指从原始数据中选择最具有代表性和相关性的特征子集以降低数据维度、减少冗余信息和噪声并提高模型的泛化能力^[28]。本文采用Pearson特征选择法对数据进行处理并检验各评价因子之间的相关性。

Pearson线性相关系数是最常用的线性相关系数^[29]，最适用数据的形式：线性数据、连续且符合正态分布、数据间差异不能太大。设有$m$个对象，$n$列指标，可以构成数据矩阵$X={\left(x_{\dot{i}j}\right)}_{m\times n}$，现在研究数据矩阵中第$a$列$x_a$和第$b$列$x_b$的相关性，设$a,b$列的相关系数为$rho(a,b)$，则

其中$m$为每列的长度，相关系数的值的范围从−1到+1，值−1表示完全负相关，值+1表示完全正相关，值越接近0代表相关性越低。

2.3   样本选取

将研究区划分为30 m × 30 m的栅格，以栅格单元作为评价单元^[30]，使用ArcGIS进行栅格转点，提取各评价因子的属性值。采用研究区内532处滑坡点作为灾害样本，并使用ArcGIS对滑坡点周围设置400 m缓冲区，按照间隔200 m生成等量非灾害样本532处，获得样本点共计1064处，建立样本库并导入MATLAB软件。按照8∶1∶1的比例划分为训练集、验证集、测试集，在训练集上采用优化算法对BP模型进行超参数优化，从而获取最优模型，在验证集上进行巡回验证，最后在测试集上测试模型性能。使用MATLAB工具对模型进行SHAP分析并导出相关参数。

2.4   模型建立与优化

2.4.1   神经网络模型（BP）

BP神经网络（backpropagation，BP）是一种常用的人工神经网络模型^[31]，其流程包括：数据准备、网络结构设计、前向传播、计算误差、反向传播、重复迭代和输出结果。

2.4.2   贝叶斯优化

贝叶斯优化（Bayesian optimization）是一种用于优化模型超参数的方法^[32]。

假设训练出的模型为$F_{\mathrm{B}}$，超参数为$\theta$，模型训练是由训练数据$(X_{\mathrm{train}}^{\prime },y_{\mathrm{train}}^{\prime })$以及超参数训练而成$F_{\mathrm{B}}=f(X_{\mathrm{train}}^{\prime },y_{\mathrm{train}}^{\prime },\theta )$，$(X_{\mathrm{train}}^{\prime },y_{\mathrm{train}}^{\prime })$分别指训练数据的两列，适应度函数为准确率最大，分类错误率最小，因此将错误率作为适应度函数，由下列计算得到，设经过训练的模型预测得到值为${\breve{y}}_{\mathrm{vaild}}^{\prime }$，则${\breve{y}}_{\mathrm{vaild}}^{\prime }=F_{\mathrm{B}}{X}_{\mathrm{vaild}}^{\prime }$，$fitness=1-{\displaystyle \frac{1}{N}}{\displaystyle \sum |{\breve{y}}_{\mathrm{vaild}}^{\prime }=={\breve{y}}_{\mathrm{vaild}}^{\prime }|}$。使用贝叶斯优化BP模型工作流程图如图3（a）所示。

图  3  算法优化流程图

Figure  3.  Algorithm optimization flowchart

下载: 全尺寸图片 幻灯片

2.4.3   麻雀搜索算法（SSA）

麻雀搜索算法（sparrow search algorithm，SSA）是一种新型的群智能优化算法，主要是受麻雀的觅食行为和反捕食行为的启发^[33]。在SSA中，具有较好适应度值的发现者在搜索过程中会优先获取食物^[34]。此外，因为发现者负责为整个麻雀种群寻找食物并为所有加入者提供觅食的方向。在每次迭代的过程中，发现者的位置更新描述如下:

式中：t——当前迭代数；

${iter}_{max}$——一个常数，表示最大的迭代次数；

$x_{\dot{i}j}$——第$i$个麻雀在第$j$维中的位置信息；

$\alpha \in (0,1]$、$ST(ST\in (0.5,1])$——预警值和安全值；

$Q$——服从正态分布的随机数；

$L$——一个$1\times d$的矩阵，其中该矩阵内每个元素全部为1。

加入者（追随者）的位置更新描述如下：

式中：$X_{\mathrm{worst}}$——当前所发现的最差的位置；

$A$——一个$1\times d$的矩阵，其中每个元素随机赋值为1或者－1，并且$A^{+}=A^T{\left(A{A}^T\right)}^{-1}$。

当意识到危险时，麻雀种群会做出反捕食行为，其数学表达式如下：

式中：$X_{\mathrm{best}}$——当前的全局最优位置；

$\beta$——步长控制参数，是服从均值为0，方差为1的正态分布的随机数；

$K\in [-1,1]$——一个随机数；

$f_i$——当前麻雀个体的适应度值；

$f_g$、$f_w$——当前全局最佳和最差的适应度值；

$\epsilon$——最小的常数，避免出现分母为0；

$K$——麻雀移动的方向同时也是步长控制参数。

2.4.4   大猩猩优化算法（GTO）

大猩猩优化算法（gorilla troops optimizer，GTO）是一种基于大猩猩部落行为的优化算法，它模拟了大猩猩部落中的社交和协作行为以寻找最优解^[35]。在GTO算法中，使用参数（$P$）来选择迁移到未知位置，当$rand<p$时迁移到未知位置，当$rand⩾0.5$则选择移动向其他大猩猩，当$rand<0.5$时选择迁移到已知位置。这三种机制为GTO算法提供了基础，第一种机制可以使算法控制问题所在的全部范围，第二种机制增强了GTO算法的全局搜索能力，第3个机制可以避免GTO算法陷入局部优化困境，3种机制的描述如下：

式中：$GX(t+1)$——下一个t次迭代的大猩猩候选位置 向量；

$X\left(t\right)$——大猩猩的当前位置向量；

$r_1、r_{2、}{r}_3、rand$——每次迭代中更新的从0到1 的随机值；

$p$——优化前必须给定的参数，取值范围为0～1，该参数决定了迁移到未知位置的概率；

$UB、LB$——变量的上界和下界；

$X_r$——整个种群中随机选择的一个大猩猩；

$G{X}_r$——随机选择的一个大猩猩，其中一个大猩 猩的候选位置向量是随机选择的，并且 包含在每个阶段更新的位置。

式（5）中的C、L和H则用下式计算：

式中：$It$——当前迭代值；

$MaxIt$——优化操作中的迭代总数；

$\cos$——余弦函数；

$r_4$——每次迭代更新的0～1的随机值。

麻雀搜索算法与大猩猩优化算法均为智能优化算法，其对BP模型的优化流程如图3（b）所示。

2.5   模型精度预测与性能评价

2.5.1   K折验证

K 折验证（K-fold cross-validation）是一种用于评估机器学习模型性能的技术，特别是在样本量较小的情况下。它通过将数据集多次划分为训练集和验证集，来获取模型的平均性能，从而更可靠地估计模型的泛化能力^[36]。

2.5.2   分类衡量指标

（1）准确率（accuracy）：范围从0到1，值越接近1表示分类器性能越好。但在不平衡类别分布的情况下，高准确率可能并不意味着好的分类器。

（2）精确率（precision）：范围从0到1，值越接近1表示分类器预测的正例中真正例的比例越高。当假正例的成本较高时，精确率是一个重要的指标。

（3）召回率（recall）：范围从0到1，值越接近1表示分类器对正例的识别能力越高。当假负例的成本较高时，召回率是一个重要的指标。

（4）F1分数（F1 score）：范围从0到1，是精确率和召回率的调和平均值。F1分数是一个综合性能指标，其表明了模型在精确率和召回率之间的平衡性能。

（5）特异性（specificity）：范围从0到1，值越接近1表示分类器对负例的识别能力越好。

（6）AUC值：范围从0到1，AUC值越接近1表示分类器对正例和负例的区分能力越好，即更能正确地对样本进行排序。

2.6   SHAP分析

SHAP方法通过观察模型在有和没有该特征时的行为来估计特征的影响^[37]。为此，将功能归属问题转化为合作博弈论问题，SHAP值是在缺失属性（缺失特征的归属为零）下计算每个特征对最终输出的贡献的唯一一致的局部精确解。SHAP值提供了一个数值，表明当考虑其他所有特征的可能组合时，特征对模型预测的平均贡献。正值表示特征对预测结果正向影响，负值表示对预测结果负向贡献^[38]。对于特征集合S中的特征$i$，shapely 值计算如下：

式中：$\text{boldsymbol}N$——所有特征的集合；

$\text{boldsymbol}S$——不包含特征$i$的任何一个特征子集；

$\left|\text{boldsymbol}S\right|$——$\text{boldsymbol}S$中的特征数；

$v\left(\text{boldsymbol}S\right)$——特征集$\text{boldsymbol}S$对模型预测输出的贡献；

$v\left(\text{boldsymbol}S\cup \left\{i\right\})\right)$——包含特征$i$的特征集$\text{boldsymbol}S\cup \left\{i\right\}$对预测 输出的贡献。

3.   研究结果

3.1   评价因子特征选择

各评价因子的数据符合该特征选择所需要的数据形式，因此使用Pearson特征选择分析后得到的图像如图4所示；其中对各评价因子相关性绝对值进行对比，最大值为道路-高程：0.472，低于0.5，满足独立性要求，可以作为评价因子进行使用。根据Pearson特征选择法对各评价因子进行数据处理，得到特征重要性值排序如表2所示，该特征选择法选择高程、道路、岩性、归一化植被指数、水系、坡度以及平面曲率作为模型评价因子。

图  4  因子相关性分析图

Figure  4.  Factor correlation analysis chart

下载: 全尺寸图片 幻灯片

表  2  因子重要性特征值表

Table  2.  Table of Feature importance values of evaluation factors

特征	高程	距道路距离	岩性	归一化植被指数	距水系距离	坡度	平面曲率	剖面曲率	坡向
特征值	0.4124	0.3264	0.3079	0.2971	0.2496	0.09382	0.09342	0.08059	0.05456

下载: 导出CSV 

| 显示表格

3.2   滑坡易发性区划

将研究区内所有栅格中每个评价因子的值输入贝叶斯-BP、SSA-BP和GTO-BP这3个优化模型中，预测每个栅格可能发生的概率，并将概率值作为对应栅格的易发性指数，按照自然段点法将其划分为：极低易发区、低易发区、中易发区、高易发区和极高易发区如图5所示^[39]。并统计已知的滑坡灾害点所处易发区的位置以及各个易发区的面积占比，如表3所示：

图  5  易发性区划图

Figure  5.  Landslide susceptibility zoning maps

下载: 全尺寸图片 幻灯片

表  3  易发区分区统计

Table  3.  Susceptibility zoning statistics

模型	面积及滑坡数量	极低	低	中	高	极高
贝叶斯-BP	面积占比/%	27.91	24.24	19.44	16.76	11.64
	滑坡数量/个	21	43	86	153	229
	滑坡数占比/%	3.95	8.09	16.16	28.76	43.04
SSA-BP	面积占比/%	33.61	22.07	16.33	14.35	13.63
	滑坡数量/个	30	39	61	124	278
	滑坡数占比/%	5.64	7.33	11.46	23.31	52.26
GTO-BP	面积占比/%	29.38	22.69	17.21	14.77	15.94
	滑坡数量/个	23	38	46	136	289
	滑坡数占比/%	4.33	7.14	8.64	25.56	54.33

下载: 导出CSV 

| 显示表格

由表3可知，在贝叶斯-BP、SSA-BP及GTO-BP模型的预测下，高易发区及极高易发区的滑坡灾害点占总滑坡灾害点的71.8%、75.57%和79.89%。

3.3   K折检验结果对比

取K=10，将整个数据集划分为10个子集，再进行10次训练和验证，每次使用9个子集为训练集，剩下一个子集为验证集，得到每次验证值如图6所示：

图  6  K折验证对比图

Figure  6.  K-fold verification comparison chart

下载: 全尺寸图片 幻灯片

由图6可知BP模型、贝叶斯-BP模型、SSA-BP模型及GTO-BP模型的10折均值分别为：0.64218、0.6751、0.73431、0.74359，3个优化模型较BP单模型的提升比例分别为：5.126%、14.346%、15.792%。

3.4   模型精度评价

模型精度评价包括AUC值、F1分数、准确率、精确率、召回率和特异性^[40]，根据各项指标表现，反应模型预测性能以及精度评价。模型精度评价如表4：

表  4  模型精度检验表

Table  4.  Model performance evaluation table

模型	AUC	提升比例/%	准确率	提升比例/%	精确率	提升比例/%	召回率	提升比例/%	F1值	提升比例/%	特异性	提升比例/%
BP	0.7608	—	0.6981	—	0.6962	—	0.7050	—	0.7006	—	0.7050	—
贝叶斯-BP	0.8653	13.7	0.7924	13.5	0.7852	12.8	0.7954	12.8	0.7903	12.8	0.7954	12.8
SSA-BP	0.8452	11.1	0.7924	13.5	0.7837	12.5	0.7913	12.2	0.7875	12.4	0.7913	12.2
GTO-BP	0.8720	14.6	0.8113	16.2	0.8024	15.3	0.8069	14.5	0.8047	14.9	0.8069	14.5

下载: 导出CSV 

| 显示表格

3.5   SHAP分析结果

利用SHAP算法对研究区滑坡易发性结果进行分析，得到评价因子重要性结果如图7所示：

图  7  SHAP分析图

Figure  7.  SHAP analysis chart

下载: 全尺寸图片 幻灯片

3.6   SHAP分析的模型局部解释

选择研究区内灾害样本和非灾害样本作为研究对象，输出模型的局部解释图并验证易发性评价结果准确度。灾害样本选择为邓家屋场滑坡（第76号灾害样本）、潘家湾滑坡（第175号灾害样本）、大块田滑坡（第363号灾害样本）及水田湾滑坡（第500号灾害样本），非灾害样本选择为第18号非灾害样本、第101号非灾害样本、第200号非灾害样本及第400号非灾害样本。灾害样本和非灾害样本分别如图8和图9所示：

图  8  灾害样本图

Figure  8.  Disaster sample map

下载: 全尺寸图片 幻灯片

图  9  非灾害样本图

Figure  9.  Non-disaster sample map

下载: 全尺寸图片 幻灯片

邓家屋场滑坡地处秭归县梅家河乡龟坪村三组。该滑坡面积2×10⁴ m²，体积10⁵ m³，为中型滑坡；主要诱发因素为降雨作用，变形迹象以后缘地表裂缝、中部土体下座、前缘滑移、地面隆起为主，现状较稳定，近期变形不明显，趋势较稳定。未造成直接经济损失，潜在威胁3户11人，威胁经济损失87万元，险情等级小型。模型预测值为0.539，属于中易发区域。利用SHAP算法对该滑坡进行区域分析，可得结果如图10（a），邓家屋场滑坡分布高程为590～660 m，地层岩性为混合沉积岩，归一化植被指数值为0.47，距水系距离为大于2000 m。图10（a）表明：道路与归一化植被指数对滑坡发生有明显促进作用，而水系对滑坡发生有明显抑制作用，结合该滑坡距水系距离，两者相对应。

图  10  灾害样本局部解释因子分析图

Figure  10.  SHAP-based local interpretation for disaster samples

下载: 全尺寸图片 幻灯片

潘家湾滑坡地处秭归县水田坝乡严坪村八组，滑坡面积9.14×10⁴ m²，体积27.43×10⁴ m³，为中型滑坡；变形迹象主要为前缘发育拉张裂缝，缝宽0.1～0.3 m，延伸长20～30 m，裂缝走向330°；潜在威胁15户40人，威胁经济损失1642万元，险情等级中型。现状较稳定，近期变形不明显，趋势较稳定。模型预测值为0.718，属于高易发区域。利用SHAP算法对该滑坡进行区域分析，可得结果如图10（b），潘家湾滑坡分布高程445～628 m，地层岩性为混合沉积岩，归一化植被指数值为0.69，距水系距离为500～1000 m。图10（b）表明：距水系距离、归一化植被指数值、高程和距道路距离对滑坡的发生有较为明显的促进作用，而地层岩性和坡度对滑坡的发生有抑制作用。

大块田滑坡位于秭归县水田坝乡，该滑坡面积25.65×10⁴ m²，体积2.565×10⁶ m³，为大型滑坡，潜在威胁17户51人，威胁经济损失2330万元，险情等级中级。模型预测值为0.671，属于高易发性区域，而该滑坡现状较稳定，变形不明显，模型预测与实际情况较为符合。利用SHAP算法对该滑坡进行区域分析，可得结果如图10（c），大块田滑坡分布高程440～780 m，地层岩性为硅质碎屑沉积岩，归一化植被指数值为0.52，距水系距离小于1000 m。基于图10（c）可以发现，距水系距离、距道路距离、归一化植被指数、高程和地层岩性对滑坡的发生有促进作用，而剖面曲率因子对滑坡的发生起抑制作用，其中岩性为碎屑沉积岩，在外力作用下容易发生破坏。

水田湾滑坡地处秭归县泄滩乡牛口村二组，该滑坡面积6.9×10⁴ m²，体积1.725×10⁶ m³，为大型滑坡；毁坏道路造成直接经济损失5万；潜在威胁对象为居民点、公路和柑橘林，共21户88人，威胁经济损失1815万元，险情等级中型。主要诱发因素为降雨和库水，变形迹象为前缘拉张裂缝，现状较稳定，近期变形不明显，趋势较稳定。模型预测值为0.822，属于极高易发区。利用SHAP算法对该滑坡进行区域分析，可得结果如图10（d），水田湾滑坡分布高程为155～261 m，该滑坡涉水，地层岩性为松散沉积物，归一化植被指数值为0.71。基于图10（d）可知，水田湾滑坡的主控因素为水系与地层岩性，该滑坡涉水，与SHAP分析结果相吻合。

第18号非灾害样本位于秭归县漆树坪村，其高程为1000～1200 m，地层岩性为混合沉积岩，距水系距离为500～1000 m，归一化植被指数值为0.79。根据图11（a）可以发现，除水系作用和高程作用外，其他因子对滑坡的发生均起到抑制作用，如距道路距离较远，不易受来往车辆震动和道路施工影响；归一化植被指数值较大，植被覆盖良好，有效抑制滑坡发生。

图  11  非灾害样本局部解释因子分析图

Figure  11.  SHAP-based local interpretation for non-disaster samples

下载: 全尺寸图片 幻灯片

第101号非灾害样本位于秭归县白云山村，其高程大于1200 m，地层岩性为碳酸盐沉积岩，距水系距离大于2000 m，距道路距离为1000～1200 m，归一化植被指数值为0.76。根据图11（b）可以发现，归一化植被指数、平面曲率及高程三个因子对该滑坡有促进因素，其他因子对滑坡的发生均起到抑制作用，如距道路距离、距水系距离均较远，不易受来往车辆震动和道路施工影响。

第200号非灾害样本位于秭归县赵家山村，其高程大于1200 m，地层岩性为混合沉积岩，距水系距离为1500～2000 m，距道路距离大于2400 m，归一化植被指数值为0.46。根据图11（c）可以发现，归一化植被指数因子对该滑坡有较明显促进因素，距道路距离、距水系距离以及地层岩性对滑坡的发生起到较明显抑制作用。

第400号非灾害样本位于秭归县六家包村，其高程为600～800 m，地层岩性为碳酸盐沉积岩，距水系距离为500～1000 m，距道路距离大于2400 m，归一化植被指数值为0.77，坡度处于20°～30°。根据图11（d）可以发现，坡度和距水系距离对该滑坡有较明显促进因素，距道路距离、归一化植被指数值以及地层岩性对滑坡的发生起到较明显抑制作用。这与第400号非灾害样本位置相符合，再根据各因子影响值，可以初步判定该样本点较为稳定。

4.   讨论

本文采用了3种不同的优化算法结合BP神经网络模型预测秭归县易发性。在对模型预测的易发性分区进行统计发现3种优化模型滑坡数量及占比均随着易发性的提高而增加，这说明3种优化模型均有效预测了研究区滑坡发生的概率，而在3种优化模型内部对比时发现，GTO-BP模型的预测效果最佳，对已知灾害样本的统计效果最好。同样，K折验证结果趋势及模型精度评价结果趋势与分区统计结果趋势相吻合，模型精度评价结果中个优化模型在提升模型性能方面有较为优异的效果，提升率均在10%以上，其中GTO-BP模型的各项评价指标提升率都最高。综合来看，GTO-BP模型在优化BP模型方面的效果有更为突出的表现。

通过4个灾害样本和4个非灾害样本的SHAP分析对比，证明SHAP算法的局部解释与滑坡所分布的特性相吻合，基于SHAP分析的局部解释可以用于解释单个滑坡的成因。本文研究结果显示：高程、地层岩性、归一化植被指数值、距道路距离及距水系距离对模型的影响较为明显，这与经验认知和前人研究结论相符合，在对滑坡进行可解释分析时，可重点关注这些评价指标。SHAP分析方法可以提高模型的透明度，展现评价因子重要性程度，因此在滑坡易发性分析时可以考虑结合SHAP分析从而提高研究结果的可利用程度。

在研究过程中发现，选择的灾害样本数量对模型的效果有着很大的影响，当前随着各位专家学者对模型的不断优化，模型的性能越来越强，如运行时间大幅度下降、运行效果提高等，模型对于样本数据的要求也越发严格^[41]，如果样本数据过少，易出现欠拟合的状态，反而优化后的模型表现欠佳，表现在AUC值下降等，当优化后的AUC值为0.4～0.7时，可以考虑两方面因素：一是优化算法与所选择模型的匹配度不高或有冲突；二是所选择的灾害样本数量过少，出现欠拟合的情况。

随着人工智能的兴起和在各领域的深度应用，将模型和研究相结合已经成为了必由之路^[42]，为了使模型在训练数据上的效果尽可能好，我们就需要使用优化算法帮助我们平衡模型的过拟合与欠拟合。因此优化算法和机器学习模型的结合成为目前的一个主流趋势^[24]，不仅在地质灾害易发性研究中，在其他研究领域各位专家学者也在积极将新算法和新模型相结合。没有一种优化算法可以适用于所有研究区和所有模型，因此模型和算法之间的匹配是后续研究的一个重点方向，可以考虑不同算法的组合，结合各个算法的优势点，从而获得更好的预测结果。同时，在进行算法与模型的耦合过程中可以加入注意力机制从而增强模型重点关注有效信息的能力。

5.   结论

（1）本文将贝叶斯优化算法、麻雀搜索算法（SSA）、大猩猩优化算法（GTO）与BP神经网络模型相结合，提出了贝叶斯-BP、SSA-BP与GTO-BP优化模型，其中GTO-BP模型有着最高的精度和模型性能，AUC值为0.8720，高于贝叶斯-BP的0.8653和SSA-BP模型的0.8452，有效提高了模型的全局搜索能力和寻优能力，更适合秭归县的滑坡易发性评价。

（2）根据GTO-BP模型，秭归县地区的极高及高易发区面积占比为30.71%，中易发区面积占比为17.21%，极低及低易发区面积占比为52.08%；极高及高易发区滑坡数占比为79.89%，中易发区滑坡数占比为8.64%，极低及低易发区滑坡数占比为11.47%。极高及高易发区主要集中于长江两岸、坡体较陡、植被覆盖稀疏及道路两侧。

（3）在对滑坡进行分析时借助SHAP算法进行可解释的滑坡成因探究，可解释的优化模型可以实现较复杂的模型训练和参数优化，将滑坡易发性预测和成因分析相结合，发现高程、地层岩性、归一化植被指数、距道路距离及距水系距离是秭归县滑坡地质灾害的主控因素，与极高及高易发区分布区域特征相吻合，在进行秭归县滑坡预测时可重点结合这几类影响因子。