基于Stacking集成学习的泥石流易发性评价:以四川省雅江县为例

2021-04-07苏刚秦胜伍乔双双扈秀宇陈阳车文超

世界地质 2021年1期

苏刚，秦胜伍，乔双双，扈秀宇，陈阳，车文超

吉林大学建设工程学院，长春 130000

0 引言

泥石流是山区最容易发生的地质灾害之一，主要由于突然性的暴雨等极端天气引发。近年来，泥石流已被公认为是全世界最主要的自然风险，会造成许多人员伤亡和经济损失，例如2008年在四川汶川就曾发生过由于地震引发的严重泥石流灾害，造成了巨大的损失[1]。雅江县隶属于四川省甘孜州，为高山深切峡谷地貌，地形起伏较大，雅砻江将研究区深切为东西两半，加上夏季暴雨影响，存在大量泥石流地质灾害。为减少泥石流灾害带来的损失，有必要对泥石流进行易发性评价，以达到预防目的，为雅江县泥石流监测预警和风险管理提供直观有效参考。

20世纪70年代以来，学者们对泥石流进行了系统的研究。在泥石流易发性评价方面，区域泥石流易发性研究通常以经验和统计方法为基础，结合3S技术展开。其中，启发式方法(如层次分析法[2])的精度主要取决于专家对某个地区泥石流发生真正原因的了解，具有一定的主观性。而一些统计学方法如信息量模型[3]和频率比[4]则通过定量分析克服了一定的主观影响。近年来，随着人工智能技术的迅速发展，机器学习技术凭借其强大的非线性建模能力广泛应用于泥石流易发性评价，如支持向量机[5]、神经网络[6]和逻辑回归[7]等。虽然人工智能及机器学习技术的快速发展，为泥石流易发性评价提供了全新的研究思路。但是，多数研究仅采用了一种单独方式进行预测，由于泥石流易发性评价问题的复杂性，单一机器学习模型可能由于随机性而导致泛化能力降低。因此，人们寻求通过集成模型来进一步提高模型的预测精度。目前，基于bagging[8]和boosting算法[9]的滑坡、泥石流易发性评价应用较多，也取得了不错的成果，特别是随机森林算法[10]。但bagging和boosting集成方式一般考虑同质分类器，如随机森林采用决策树，很少有研究结合常用的不同原理的机器学习方法来构建高质量模型。

Stacking(堆叠)是Wolpert于1992年在“Stacked generalization”论文中提出的一种异质集成学习算法[11],旨在以一种取长补短的方式结合不同机器学习方法的优势来消除误差，从而进一步提升预测精度。Stacking集成方法已在许多领域取得了显著的成果，如碰撞损伤严重程度分析[12]、电力负荷预测[13]等。

因此，笔者以泥石流多发的雅江县为研究区，以3种机器学习方法为基模型(支持向量机、神经网络和随机森林)，通过Stacking集成学习框架建立了一种多模型融合的泥石流易发性评价模型，并同时采用支持向量机、神经网络和随机森林进行对比验证，以期能够客观验证Stacking融合模型的效果，又能对雅江地区泥石流灾害预警提供直观有效的参考。

1 研究区概况

研究区(雅江县)位于四川省西北部，甘孜州腹地，其面积约为7 681.5 km2，地理坐标为100°19′～101°26′E，29°03′～30°30′N。研究区属青藏高原亚湿润气候区，全年冬长无夏，春秋相连，无霜期短，年平均降水量783.2 mm，主要集中在5—10月。研究区河流呈树枝状分布，其中雅砻江由北向南纵贯全境，将全县深切为东西两半，河谷深切呈“V”字型，主要支流有鲜水河、庆大河、曲汝河、密西沟、祝桑河、霍曲河、马岩河7条。地势总体上北高南低，东北和西北部为山原地貌；中部为河谷地貌；西南部是极高山地貌，海拔在5 000 m以上。研究区地震加速度0.15 g～0.2 g(g=9.8 m/s2)，基本烈度为Ⅶ度，地层较为单一，主要出露三叠系和新生界第四系地层[14]。土壤类型主要有淋溶土、半淋溶土、高山土及水成土。

随着西部大开发战略的实施，雅江县的城市建设日益发展。但由于研究区处于高山深切峡谷地貌，土地资源极其有限，在城镇化建设和发展中，向坡要地、建筑加载和工程切坡等人类工程活动强烈，再加上山体表层多松散岩石与土壤，在强降雨的诱发下泥石流等地质灾害极容易发生，严重威胁着当地居民的生命和财产安全[15]，如2011年暴雨诱发的“7·13四川甘孜泥石流”，造成国道318线上22处路基被摧毁，2处路基塌方，15条乡村道路被冲毁180余公里，桥涵冲毁90余座。因而有必要对研究区进行泥石流易发性评价，为雅江县泥石流预防提供参考。

本研究泥石流灾害数据来源于2010年四川省雅江县地质灾害调查数据集[15]，研究区共发育158处泥石流灾害，主要分布于研究区沿雅砻江两岸及其支流的沟谷地区，为典型的高寒高海拔山区沟谷型泥石流，按照行政区划八角楼乡分布最多，为40条，红龙乡和瓦多乡最少，各1条，其他乡均在3～17条之间，研究区泥石流分布见图1，本文按照70%、30%的比例将灾害点随机分为两组，训练点110个，验证点48个，用于模型的建立和验证。

图1 研究区泥石流分布

2 研究方法

2.1 评价因子的选取

根据研究区的实际状况和对以往文献的查阅，在充分考虑研究区资料获取难易程度及研究尺度后本文选取了12个泥石流评价因子，即高程、坡度、坡向、地形起伏度、平面曲率、剖面曲率、年平均降雨、到河流的距离、到道路的距离、地形湿度指数(Topographic Wetness Index,简称TWI)、归一化差分植被指数(Normalized Difference Vegetation Index，简称NDVI)及土壤类型，研究区因子图(图2)，各数据来源见表1。

a.高程; b.坡度;c.坡向; d.地形起伏度; e.平面曲率; f.剖面曲率; g.年平均降雨; h.到河流距离; i.到道路的距离; j.TWI; k.NDVI; l.土壤类型。

表1 数据来源

高程是易发性研究中重要的评价因子，通过影响降雨分布、植被覆盖、土地利用及土壤水分间接影响泥石流的发生，特别是对于山高谷深的地方，降雨、植被等随着高程出现垂直变化的特点[16]。研究区高程范围为2 171～5 142 m，根据等间距将其分为4类，分别为2 171～3 000 m，3 000～3 600 m，3 600～4 200 m，>4 200 m。

大多数泥石流与坡度直接相关，影响泥石流的发生、发展及其运动规模。在山高沟深、地形较陡的地方，水流比较容易汇流到一起。因此在暴雨条件下，坡度较陡的山地更容易发生泥石流灾害。研究区坡度变化范围为0°～64°，利用等间距分为5类，分别为0°～10°，10°～20°，20°～30°，30°～40°，>40°。

坡向是泥石流的重要评价因子，描述边坡暴露于阳光、风和降雨的方向，从而影响成坡材料的性能，与阴坡相比，阳坡冰雪融化快、岩石风化较为强烈，所以阳坡泥石流发生概率大于阴坡。由于研究区起伏较大，几乎没有平坦区域，所以将研究区坡向分为8类，分别为0°～22.5°和337.5°～360°为正北，22.5°～67.5°为东北，67.5°～112.5°为正东，112.5°～157.5°为东南，157.5°～202.5°为正南，202.5°～247.5°为西南，247.5°～292.5°为西，292.5°～337.5°为西北。

地形起伏度是指地表某一局部区域内地形高度的差异，直接代表地形的陡峭程度，是地貌的重要判断指标。地形起伏度通过ArcGIS焦点统计功能计算，邻域设置为高100 m、宽100 m。研究区地形起伏度变化范围为0～170 m，其分为5类，分别为0～20 m，20～35 m，35～50 m，50～65 m，>65 m。

曲率是指地表任意一点的变形状态。平面曲率为曲率水平分量，负值表示下凹，正值表示上凸。剖面曲率是地面坡度的导数，描述了地表在垂直方向的弯曲状态，决定物质相对运动速度，正值表明加速，负值表明减速，反映地面土壤侵蚀的程度。平面曲率根据自然间断点法分为3类:-4.19～-0.1，-0.1～0.1，-0.1～0.5；剖面曲率同样根据自然间断点法分为3类:-7.87～-0.1，-0.1～0.1，0.1～7.46。

突发性水源是泥石流的重要诱因，为松散固体源提供了动力条件，降雨造成的地表水在冲洗疏松沉积物的过程中逐渐演变成泥石流。此研究应用ArcGIS克里金插值(高斯函数)获得了年平均降雨量分布图，并将降雨根据等间距分为5类:687～720 mm，720～740 mm，740～760 mm，760～780 mm，780～805 mm，从东南向西北逐渐减小。

河流为泥石流提供了重要的水文基础和流通通道。此研究利用欧氏距离建立了河流缓冲区以分析泥石流和河流之间的关系，并将到河流距离分为4类:0～200 m，200～2 000 m，2 000～5 000 m，>5 000 m。

到道路的距离在一定程度上反映了人类活动的影响，如建造道路、砍伐树木和人工削坡都可能造成坡面岩体裸露、风化严重，从而产生泥石流的物源及流通基础。将到道路距离分为4类:0～200 m，200～2 000 m，2 000～5 000 m，>5 000 m。

TWI是对径流长度和径流面积的定量描述，直接反映土壤含水量随空间变化的规律，也间接反映研究区的水土流失状况。研究区地形湿度指数根据自然间断点法分为4类:-0.62～6.45，6.45～9，9～14，14～35。

NDVI也是泥石流易发性评价常用的评价因子，是反映出植被空间分布密度的最佳因子。植被的根系可以对土壤进行加固，可以预防土壤的侵蚀，而水土流失正是泥石流的重要诱因之一。NDVI值在-1到1之间，正值表示植被覆盖，而且随着NDVI的增加植被覆盖越密；0表示为岩石或裸土等;负值表示地面是云、水和雪等。研究区NDVI根据自然间断点法分为4类:-0.58～0，0～0.2，0.2～0.39，0.39～0.71。

物源条件是泥石流发育的必要条件，岩土类型决定了泥石流灾害的规模和破坏形式，但研究区地层比较单一，以三叠系及新生界第四系地层为主，故决定用土壤类型代表研究区的物源条件[14]。研究区土壤类型有淋溶土、半淋溶土、水成土和高山土，其中水成土占地面积仅0.43%，且无泥石流，故建模时未考虑。

2.2 信息增益比

在泥石流易发性评价建模中，并不是所有的评价因子都能够较好地预测泥石流分布，甚至一些预测能力较差的评价因子会反向影响模型的性能，降低模型的预测能力。因此，有必要对评价因子的预测能力进行定量研究，以消除预测能力不强或很弱的评价因子。其中，信息增益比在选择因子中具有广泛的适用性，有助于评价和排序输入变量的重要性，以信息理论为基础，在信息增益的基础上增加了一个惩罚系数，跟踪信息熵的减少来量化评价因子的重要性，是数据挖掘中衡量因子预测能力的标准技术[16]。

已知数据训练样本T由n个样本组成，那么样本T在自变量F(泥石流评价因子)下的信息增益比IGR(T,F)可以通过以下公式进行计算:

(1)

IG(T,F)=Info(T) Info(T,F)

(2)

(3)

(4)

(5)

其中，n(Li,T)表示在T训练集中变量第i类(如泥石流L1，非泥石流L2)中样本的数量。Tj(j=1,2,…,m)是在T训练集下某一自变量F(泥石流影响因子)的第j类的数量。

Info(T)为各类别的总体信息熵，Info(T,F)是样本T在自变量F下的信息熵。IG(T,F)为各评价因子的信息增益，Spli(T,F)表示分裂信息度量，是将训练样本T划分为m个子集所产生的潜在信息值。

2.3 Stacking集成方法

Stacking是一种异质集成学习算法[11]，异质集成是指通过组合若干原理不同的分类器，来组建一个强分类器，以增强泛化能力。此研究Stacking算法采用一种两层结构，Stacking学习框架如图3所示，第一层的分类器称为基模型，第二层的分类器称为元模型，具体训练过程为：

(1)将原始数据分为训练集和测试集；

(2)利用训练集训练第一层的若干基模型，分别得到若干预测结果(使用五折交叉验证)；

(3)将第一层得到的若干预测结果作为第二层模型的输入特征，而训练集的标签仍然作为新数据集的标签；

(4)使用新数据集训练第二层的模型得到最终的预测结果。

图3 Stacking学习框架

元模型的训练集由基模型的预测值组成，如果直接用基模型的训练集来生成第二层训练集，就导致数据在两层模型中重复使用，产生过拟合。为了避免这种现象，有必要对数据使用进行规划，对每个基模型采用五折交叉验证[13]，具体过程为：

(1)将原始训练数据集随机分成数量相同的五块子数据集，每个数据块不能重叠；

(2)对于单个基模型，轮流使用4个数据块作为训练集，剩下一个作为测试集，每个基模型针对自己的测试集都可输出一个预测结果，并且这5个结果可以拼成新的数据集，新数据集与原始训练集大小相同；

(3)将生成的新数据集作为元模型的输入特征，从而实现了所有训练数据从原始特征到第二层特征的变换，而且每个基模型预测的数据均未参与到该模型的训练。

这种数据使用令所有数据在训练中仅使用一次，能有效防止训练数据的泄露，避免过拟合。

2.4 Stacking采用的机器学习模型

立足于基模型的预测能力及差异性，在Stacking融合模型第一层应选择预测能力强且有较大理论差异的模型作为基模型。第二层应选择泛化能力强的模型，将多个基模型对于训练集的偏置找出并进行纠正。因此，Stacking融合模型第一层基模型分别为支持向量机、神经网络和随机森林，第二层选择逻辑回归作为元模型，这几种算法不仅有充分的理论支持，在泥石流易发性评价中也得到过广泛使用。

支持向量机是一种监督学习方法，是由结构风险最小化概念发展而来的具有优秀泛化能力的机器学习方法，支持向量机设计了一种非线性变换，令训练数据通过核函数投影到高维特征空间，使得投影后的训练数据线性可分[5]。近年来，支持向量机逐渐在泥石流易发性评价中获得广泛使用。参数设置为，模式选择专家模式，使用RBF作为内核，停止标准为0.001，规则化参数为10，gamma为0.1。

神经网络是在模仿人脑神经冲动传导原理基础上建立的一种非线性数学模型，输入层、隐含层和输出层组成了神经网络的结构，每层由神经元组成，也可称之为节点[17]。神经网络应用于泥石流易发性评价具有不错的表现。参数设置为，结构为1-2-1层，输入层节点数与泥石流条件因子数相同，输出层节点数为2。隐含层节点数取输入节点数(12)和输出节点数(2)之积的平方根，取5。最后一层用S形激活函数输出泥石流发生的概率。

随机森林是集成学习中一种袋装思想的算法策略。bagging方法通常考虑同质弱学习器，独立、并行地学习这些弱学习器，并按照一定的确定性平均过程进行组合。随机森林以决策树为基本学习器，建立了袋装集成算法，同时引入了特征的随机抽样，进一步提升了模型的抗噪声能力，有效地防止了过拟合的发生，在处理高维、大数据集时，仍然可以高速处理和学习[10]。随机森林通过SPSS Modeler 18.0软件构建，参数设置为500个决策树，最大树深度为10，节点个数m为评价因子个数(12个)的平方根，故将m取值为3。

逻辑回归是二分类问题中应用最为广泛的机器学习模型，具有良好的稳定性[7]。原理是将样本的特征与样本发生的概率联系起来，计算结果是通过样本的特征来拟合计算出一个事件发生的概率。逻辑回归通过SPSS Modeler 18.0软件构建，参数设置为采用多项式逻辑回归，方法选择进入法，模型类型为主效应。

2.5 受试者工作特征曲线

受试者工作特征曲线又称ROC曲线(Receiver Operating Characteristic Curve),通过选择所有阈值来绘制整条曲线，该方法曾广泛用于验证二分类模型的性能[8,10]，ROC曲线下面积(Area Under Curve，简称AUC)代表模型的可预测性，其值范围是0.5～1，1代表是完美模型，0.5则表示是较差的模型[17]。利用训练集数据可以得到模型成功率ROC曲线，曲线下面积即模型的成功率；测试集数据得到预测率ROC曲线，用于验证模型的预测是否准确，曲线下面积就是预测率。

3 结果与讨论

3.1 信息增益比

评价因子信息增益比结果(保留四位小数)如图4所示。其中有4个因子有较强的预测能力，分别为到河流距离、到道路距离、高程、剖面曲率。这和研究区的实际情况是一致的。雅砻江上游水源丰富，而且大多数河流都位于山谷中，当暴雨发生时很容易形成泥石流。河流为泥石流提供了重要的水文基础和水流通道。同样，到道路的距离很重要。道路在某种程度上反映了人类工程活动和地貌环境的影响。高程和剖面曲率则代表了地形因素对泥石流的影响。而坡向(0.007 3)、平面曲率(0.004 6)这两个因子均<0.01，预测能力可以忽略，因此决定在之后的建模中去掉这两个评价因子。

图4 评价因子信息增益比结果

3.2 易发性分区结果

本研究使用自然间断点分级法对泥石流易发性图进行分类，共分为5个区域，分别为极低易发区、低易发区、中易发区、高易发区和极高易发区，最后得到泥石流易发性图如图5所示。

这4个模型制作的易发性图从总体上看趋势相同，主要区别在雅江县西侧及南侧区域。西部区域如德差乡、柯拉乡、西俄罗镇和红龙乡等地，在支持向量机和神经网络的易发性图中，绝大部分是中等以下区域，而随机森林则有局部地区为高等以上。南部地区主要是恶古乡地区，神经网络的易发性图极高等占比较大，而随机森林、支持向量机的易发性图显示出较低的易发性，而Stacking融合模型则综合了这3个基模型的特征，修正了各模型的误差。

Stacking融合模型制作的易发性图的高易发和极高易发区主要集中在由北向南纵贯全境的雅砻江全线两岸及其部分支流附近，中部主要是雅砻江沿岸，依次经过普巴绒乡、呷拉镇、雅江县城、米龙乡、麻郎措乡、恶古乡、八衣绒乡、波斯河乡、吉居乡、牙衣河；其中西部包括德差乡大部分，柯拉乡的霍曲河附近，西俄罗镇的吉珠沟附近，红龙乡的东北地区；东部包括八角楼乡的中部地区，祝桑乡的大部分地区。

高易发区地表沟谷发育，泥石流灾害点距河流距离分布大多<2 000 m，泥石流发育多沿沟谷流通，致使大型河流对研究区的影响较高。特别是雅砻江沿线河流深切，山体高大，谷坡陡峻，河床狭窄，可以为泥石流提供势能，同时提供泥石流运移的通道。高程大部分分布在2 171～3 600 m，少部分位于3 600～4 200 m，4 200 m以上几乎都是极低易发区。坡度分布为20°～40°。此段区域主要为低山至河流过渡区域，高差较大，坡度适中，能够为泥石流提供势能。土壤类型主要为淋溶土和半淋溶土。NDVI分布为0～0.39，地表主要是植被较少、风化严重的裸地。

3.3 泥石流易发性与实际对比分析

模型分区面积占比情况见表2，模型分区泥石流灾害点分布情况见表3。4个模型生成的易发性图的空间分布模式相似(图5)。大部分高易发性的网格单元集中在雅砻江沿线的河谷中。所有模型在极高和高易发性级别产生的网格单元都较少，而极低和低的易发性级别所占比例最大(表2)。

a.支持向量机; b.神经网络; c.随机森林; d.Stacking融合模型。

表2 模型分区面积占比情况

表3 模型分区泥石流灾害点分布情况

由表2～3可知，Stacking融合模型制作的易发性图大约有13.13%区域对泥石流易发性极低，其余40.16%，23.49%，18.11%，5.11%的区域分别表现为低、中、高和极高易发性。位于极低、低、中、高、极高易发性区的泥石流点分别有0，3，5，42个及108个，位于高易发性区以上的泥石流点占94.94%，位于极高易发性区以上的泥石流点占68.35%。3个基模型生成的泥石流易发性图位于高易发性区以上的泥石流点占比由大到小分别为随机森林(93.04%)、支持向量机(90.51%)和神经网络(89.87%)，位于极高易发性区以上的泥石流点占比由大到小分别是神经网络(64.56%)、随机森林(62.66%)和支持向量机(59.49%)。相比之下，Stacking融合模型生成的泥石流易发性图在这两部分的泥石流点占比均高于3种基模型，有相当高的准确率，这和ROC曲线的预测率有较好的一致性。

3.4 ROC曲线验证结果

使用训练集的3种模型的ROC曲线见图6，使用验证集的3种模型的ROC曲线见图7，模型成功率及预测率见表4。各模型预测率分别为支持向量机(0.897)，神经网络(0.906)，随机森林(0.916)，Stacking融合模型(0.955)。Stacking融合模型通过集合各个模型的优势，舍弃基模型的误差，预测率达到了0.955，相对3个基模型分别提升了支持向量机(0.058)，神经网络(0.049)，随机森林(0.039)。