0.   引言

泥石流是山区一种常见且极其严重的地质灾害，它是由水和大量松散物质（如沉积物、碎屑、泥沙等）组成的混合物，具有极强的破坏能力，在世界各地的山区造成巨大的人员伤亡和经济损失^[1]。汶川大地震触发了大量崩塌和滑坡，为泥石流发育提供了丰富的物源，地震灾区泥石流活动进入了活跃期^{[2 − 5]}，研究泥石流易发性及影响因子的贡献率，旨在为未来该区域泥石流灾害的防灾减灾提供参考依据。

以往对于区域泥石流易发性评价通常依据主观经验和统计方法，并结合遥感和地理信息系统展开，主要的方法包括：层次分析法，模糊综合评价法，信息量模型，贡献权重叠加模型等^{[6 − 9]}。随着人工智能的迅速发展，机器学习凭借其强大的非线性建模能力被广泛地运用于灾害领域，如滑坡、山洪、森林火灾以及少量的泥石流易发性研究等^{[10 − 13]}。与前述传统易发性评价方法相比，机器学习方法在数据处理阶段无需归一化批量处理数据、无需对评价因子权重进行主观赋值、能够客观计算评价因子贡献率及泥石流易发性概率。同时，由于机器学习模型众多且不同模型具有不同适用条件，对比分析不同机器学习模型的评价结果，并在此基础上开展泥石流易发性评价和因子贡献率及其原因分析，对于优选评价模型和提高评价结果准确性和可靠性具有重要意义。

本文以岷江上游为研究区，以小流域为评价单元，分别采用5种机器学习模型对汶川大地震前、后岷江上游泥石流易发性进行评价，优选出最佳模型，并在此基础上结合各因子的贡献率及其变化对地震前、后研究区泥石流发育规律进行分析与研究。

1.   研究区概况

岷江上游位于四川盆地向青藏高原东缘的过渡地带，在40～50 km水平范围内海拔从700 m陡变至5000 m，河流深切，地表起伏差异大，属于典型的高山峡谷区^[14]，见图1（a）。岷江上游干流全长约340 km，流域面积21580 km²，根据干流河道纵剖面特征可将干流由上游向下游分为北段、中段和南段，支流按汇入干流由南向北依次为：寿溪流域、渔子溪流域、杂谷脑河流域、黑水河流域、小姓沟流域，见图1（b）。基于Alos12.5 m DEM数据，根据ArcGIS 软件中的水文分析模块通过设定合理阈值提取小流域，结果表明5000为最佳的提取阈值。DEM自动提取的小流域在丘陵洼地等与实际不符合的区域需结合遥感数据进行手动修正，经手动校正后共提取出1848条小流域，作为泥石流易发性评价单元，见图1（c）。

图  1  研究区概况与评价单元

Figure  1.  Overview and evaluation unit in the study area

下载: 全尺寸图片 幻灯片

岷江上游地质条件复杂，新构造运动强烈，区内断层发育、地震活动频繁，地质灾害多发。“5•12”汶川大地震前岷江上游共有127条泥石流沟，受汶川大地震的影响，大量滑坡和崩塌堆积于坡脚，为泥石流发育提供了丰富的物源，震后累计新增320条泥石流沟，主要沿岷江干流及黑水河、杂谷脑河两岸分布。

2.   研究方法

2.1   易发性评价流程

如图2所示，易发性评价流程是根据不同易发性评价因子和泥石流在1848条流域的分布情况生成数据集，按7∶3的比例将其分为训练集和测试集^{[13 − 14]}，分别采用5种机器学习模型对其进行易发性评价，根据预测精度（predictive accuracy，ACC）和受试者工作特征曲线（receiver operating characteristic curve，ROC）下的面积（area under curve，AUC）选出最佳模型，并根据最佳模型确定震前震后各评价单元泥石流的易发性及不同评价因子的贡献率。

图  2  泥石流易发性评价流程

Figure  2.  Evaluation process of debris flow susceptibility

下载: 全尺寸图片 幻灯片

2.2   机器学习模型概述

本文选取具有代表性的5种机器学习模型，根据其内置数学原理和迭代方法的不同，分为浅层和集成两类，其架构如图3所示。通过ACC和ROC来对比各个模型的性能，选择出最优模型运用于泥石流易发性评价。

图  3  机器学习模型架构示意简图

Figure  3.  Schematic diagram of machine learning model architecture

下载: 全尺寸图片 幻灯片

（1）朴素贝叶斯（naive Bayes，NB）

朴素贝叶斯算法是一种经典的概率统计机器学习模型，它因应用简单和高学习效率的优点被广泛应用^{[15 − 16]}，模型的基础架构为贝叶斯定理。在本研究中，贝叶斯定理和全概率公式的结合被用于构建泥石流易发性评估模型。在模型的训练过程中，先验概率$P(y = j)$预先被计算，其次求出小流域x相对于小流域泥石流发生与否的状况类别$j$的类条件概率$ P(x\left| {y = j} \right.) $，具体过程如式（1）所示：

式中：f——小流域的特征序列，每个小流域x共包括 F维特征。

随后，将类条件概率和先验概率相乘即可得到小流域x为泥石流发生与否的状况类别j的后验概率$ P(y = j\left| x \right.) $。上述即是模型在测试过程中的原理，小流域最大的后验概率所评估的泥石流是否发生即为模型的最终输出。

（2）逻辑回归（logistic regression，LR）

逻辑回归是一种广泛运用于敏感性评估的多元回归算法^{[17 − 18]}，与属于生成性模型的朴素贝叶斯不同，它属于判别式模型；其中，模型中的因变量以数据编码的形式呈现，通常是“0”和“1”两个编码。因此，该模型适用于本研究以判定泥石流是否发生，将“0”定义为泥石流不发生，“1”则为泥石流发生。模型的构造过程可大致分为以下三步：（1）建立预测函数；（2）计算代价函数；（3）多次迭代用于模型调参，最终可在y=0和y=1之间拟合一条最佳回归线。

（3）决策树（decision tree，DT）

决策树算法通常被用于统计模型中^{[19 − 20]}；其基础架构由一个根节点、多个内部结点和叶节点组成，其中叶子节点代表分类的最终结果，标志着分支的结束；内部节点作为根节点和叶节点之间的连接纽带，代表了数据属性；根结点作为决策树模型的开端，应用“if-then”的逻辑并不断向下分支，其中分支所用到的样本特征即为影响泥石流发生的参数。在决策树的分支过程中，Gini系数作为一个标准进行特征优化，通过基尼系数，决策树模型不断深化，最终得到精确度最高的分类模型。如下式所示：

式中：$ {P_j} $——选中的小流域属于j类别的概率 。

（4）随机森林（random forest，RF）

随机森林是基于同构集成策略的集成学习模型^{[21 − 22]}，由多个决策树通过bagging算法组合而成。在本研究中，该模型的核心思想为小流域数据集的来回抽样，预先生成多个决策树模型，每个决策树模型从整体数据集中随机选取k个小流域进行训练，通过多数投票机制确立最后的分类结果（泥石流发生与否）。此外，n_estimator，即决策树的个数对于随机森林模型的准确率影响较大，应在模型应用前进行参数调优。

（5）极端梯度提升（extreme gradient boosting，XGBoost）

与上述的随机森林模型不同，XGBoost模型以gradient boosting算法^[23]作为核心思想进行基础决策树模型的集成，它最早由Chen等^[24]在2016年提出。该算法的构造过程是多棵决策树的先后拟合，每棵新决策树取以上一棵决策树的最优预测结果为构造基准，由此不断分裂出新的决策树，并将其与先前树的预测结果拟合，将多棵决策树的预测结果累加以得到最终的预测结果。为了防止过拟合现象的产生，该算法通过损失函数中的正则项控制分裂过程中的树的复杂程度，算法的目标函数如下式所示：

式中：${O_w}$——目标函数；

L——损失函数；

$ \theta $——正则化项；

${t_w}$——XGBoost构造过程中的单一决策树模型；

T——决策树的总数量；每一棵树都利用了数量为I的样本进行训练，

${y_i}$、$F({x_i})$——代表小流域${x_i}$的真实映射值和预测值。

2.3   评价因子选取与数据获取

根据小流域形成泥石流的孕灾条件从地形地貌、降雨、地质、植被覆盖和人类活动5个方面选取了震前11个因子，震后12个评价因子。

（1）地形地貌：面积坡度积分、侵蚀传递系数^[25]、面积高程积分、纵沟比降、形状系数和流域面积6个指标因子（根据12.5 m 的DEM栅格数据提取）；

（2）降雨：汶川地震前、后年均降雨量（国家气象科学数据中心获取岷江上游汶川地震前及2017年年均降雨量，为矢量数据）；

（3）地质：距断层距离、岩性、汶川地震地震烈度（由1∶ 200 000地质图提取岷江上游地区流域距断层距离和地层岩性信息；从岩性的抗侵蚀能力角度出发，根据《岩土工程勘察规范》（GB 50221—2001）^[26]，按照工程岩性分组要求将岷江上游地层岩性分为5个类别进行赋值；由国家地震科学数据中心网（http://data.earthquake.cn/index.html）获取汶川地震地震烈度的相关数据，为矢量数据）。

（4）植被覆盖：植被覆盖率（通过地理遥感生态网（http://www.gisrs.cn/）下载的遥感数据，根据公式计算出植被覆盖率，用研究区的矢量边界裁剪生成的栅格数据，得到研究区的植被覆盖率，为栅格数据）。

（5）人类活动：土地利用率（本文基于刘纪远^[27]等学者提出的土地利用类型分级原则，将研究区震前、震后各类土地利用类型进行1～4级的等级划分，为矢量数据）。

2.4   评价模型的验证

受试者工作特征曲线（receiver operating characteristic curve，ROC）是度量二分类模型评价效果的一种有效工具^[28]，其与预测精度（predictive accuracy，ACC）一起被广泛应用于验证机器学习模型的性能评价。ACC值和ROC曲线是依赖于混淆矩阵的统计度量^[29]，如表1所示，混淆矩阵揭示了模型结果与实际结果之间的差异，其中1表示泥石流发生，0表示未发生泥石流。

表  1  混淆矩阵

Table  1.  Confusion matrix

	预测结果
		1	0
真实结果	1	真阳性 （TP）	假阴性 （FN）
	0	假阳例 （FP）	真阴例 （TN）

下载: 导出CSV 

| 显示表格

ACC可通过式（4）估算：

ROC曲线说明了当阈值变化时真阳性率（true positive rate，TPR）和假阳性率（false positive rate，FPR）的变化^[30]，如式（5）（6）。通过设置不同的阈值，并将（$ FPR $，$ TPR $）绘制于二维坐标系中（其中，$ TPR $为纵轴，$ FPR $为横轴），便可得到ROC曲线，AUC值则为ROC与横轴的包络面积。理想的ROC曲线应接近左上角，通常AUC值越高，表明模型预测效果越好。根据先前的研究^{[29, 31]}，基于AUC值的模型性能可分为几个级别：0.5～0.6较差，0.6～0.7中等，0.7～0.8可接受，0.8～0.9优秀，0.9～1接近完美。

3.   泥石流易发性分析

3.1   模型优选

以小流域为单元，统计其对应的评价因子（震前11种，震后12种）和泥石流数据，分别生成震前和震后数据集（各1848条），运用5种机器学习模型分别对震前和震后泥石流易发性进行评价，并计算其对应的ACC和AUC值如表2和图4所示。

表  2  测试集中各模型ACC及AUC值

Table  2.  ACC and AUC values of the model on the test data set

		NB	RF	XGB	DT	LR
震前	ACC	0.92	0.93	0.93	0.91	0.93
	AUC	0.77	0.84	0.84	0.76	0.78
震后	ACC	0.83	0.86	0.85	0.81	0.85
	AUC	0.84	0.91	0.90	0.83	0.90

下载: 导出CSV 

| 显示表格

图  4  基于测试集的各模型ROC 曲线及AUC值

Figure  4.  ROC curves and AUC values of each model based on the test set

下载: 全尺寸图片 幻灯片

地震前后无论是浅层还是集成机器学习模型的ACC和AUC值均高于0.75，所选用的模型均表现较好。震前和震后集成机器学习模型中随机森林模型（ACC=0.93，AUC=0.84；ACC=0.86，AUC=0.91）略优于XGB模型，均高于浅层机器学习模型。集成机器学习模型相比于浅层机器学习模型，因其架构由多个分类器组合而成使得其训练过程更加聚合，方差更小，泛化能力更为优秀。综合上述5种机器学习模型在泥石流易发性评价中的表现，随机森林模型最佳。因此，后文均采用随机森林模型的评价结果对泥石流易发性和评价因子贡献率进行分析。

3.2   地震前后易发性评价结果

根据随机森林模型所计算的震前和震后泥石流易发性指数，采用自然断点法将其由低到高划分为极低、低、中、高和极高五个易发性等级，如图5所示，并对不同等级泥石流数量和发生率（发生率是指在某一易发性等级发生泥石流的小流域数量和该等级小流域总数量的比值）进行统计，如图6所示。

图  5  岷江上游地区泥石流易发性等级图

Figure  5.  The susceptibility grade of debris flow in the upper reaches of the Minjiang River

下载: 全尺寸图片 幻灯片

图  6  不同易发性等级内的泥石流数量与发生率

Figure  6.  The number and occurrence rate of debris flows within different susceptibility levels

下载: 全尺寸图片 幻灯片

震前灾害点主要分布在岷江上游干流和支流的中、下游，震后灾害点空间分布与震前大致相同，但灾害点数量明显增多，且干流下游增幅最为明显，同时，极高和高易发等级的空间分布与灾害点的空间分布特征一致（图5）。除极低易发等级以外，震后各易发性等级泥石流数量均增加，且增量随易发性等级的提高而增多；震前和震后泥石流发生率均随易发性等级的提高而增大，且等级越高增量越大，同时，震后泥石流发生率在各等级内均高于震前（图6）。

3.3   评价因子贡献率分析

基于随机森林模型，根据Python软件中feature_importances函数生成不同泥石流易发性影响因子的贡献率，计算式如式（7）：

式中：IP_i——第i个因子的贡献率；

i——泥石流易发性影响因子的序列号；

N——影响因子的总数。

结果如图7所示。可以看出：泥石流的发生是各个因子共同作用的结果，所选的12个因子对泥石流的发生都有着积极的贡献。其中地形因子6个（震前0.536，震后0.559），气象因子1个（震前0.102，震后0.108），地质因子5个（震前0.161，震后0.188），土地覆盖因子1个（震前0.062，震后0.087），人类活动因子1个（震前0.113，震后0.082）。即地形、地质和气象是影响泥石流发生的重要因子群，且除人类活动因子外，震后其他因子群的贡献率都有不同幅度的增加。对震前震后的单因子分析，除侵蚀传递系数和植被覆盖率因子以外，汶川地震后各因子贡献率均较震前有所减少，这主要是由于震前和震后所有因素的贡献率之和为1，震后地震烈度因子（贡献率为0.096）的加入导致其它因子的贡献率降低，但侵蚀传递系数贡献率增加显著。

图  7  评价因子贡献率

Figure  7.  Contribution rate of evaluation factors

下载: 全尺寸图片 幻灯片

为进一步探究泥石流的易发性与侵蚀传递系数的关系，分别对研究区各流域的泥石流数量、发生率和平均侵蚀传递系数进行统计，结果如图8所示。不同流域震后泥石流数量及发生率均高于震前，且按照干流上游向下游及支流流域汇入干流的顺序（由上游向下游），泥石流发生率和侵蚀传递系数均值均逐渐增大，表明泥石流与侵蚀传递系数的空间分布关系密切。

图  8  不同流域的泥石流数量、发生率和侵蚀传递系数均值

Figure  8.  Mean debris flow count, occurrence rate, and erosion transfer coefficient of debris flows in different watersheds

下载: 全尺寸图片 幻灯片

侵蚀传递系数代表了因构造运动导致侵蚀基准面下降后河道的溯源侵蚀过程，溯源侵蚀会从干流下游向上游、干流向支流、支流向小流域传递，从而影响河道和沟道的纵比降及高差，为泥石流发生提供动力条件。汶川地震导致大量崩塌和滑坡发生，形成了大量的松散堆积物，为泥石流的发生提供了丰富的物源，而汶川地震的震中位于研究区的河口附近，表征地震对小流域物源影响的地震烈度因子总体上从干流下游向上游逐渐降低。侵蚀传递系数和地震烈度的这种空间上的叠加作用加剧了震后干流和支流泥石流由下游向上游发育程度逐渐降低的空间分布规律。

4.   结 论

（1）地震前、后无论是浅层还是集成机器学习模型均表现较好（ACC均高于0.80，AUC均高于0.75）。集成机器学习模型均高于浅层机器学习模型，其中随机森林模型在地震前、后的表现最优，更适用于泥石流易发性评价。

（2）基于随机森林模型评价结果，震前和震后泥石流灾害点与极高和高易发等级的空间分布一致，主要分布在岷江上游干流和支流的中、下游，且在干流下游增幅最为明显。随着易发性等级的提高，震前、震后各等级内泥石流发生率均呈逐渐增大的趋势，且等级越高增量越大，震后泥石流发生率在各等级均高于震前。

（3）地震前、后侵蚀传递系数的贡献率均显著高于其它因子，与汶川大地震地震烈度空间分布特征叠加，为震后泥石流发生提供了良好的动力条件和丰富物源，加剧了震后干流和支流泥石流由下游向上游发育程度逐渐降低的空间分布规律。