基于原型网络的云南怒江州泥石流灾害易发性评价与区划

2023-11-11韩俊王保云

中国地质灾害与防治学报 2023年5期

韩俊，王保云

（1.云南师范大学信息学院，云南昆明 650500；2.云南师范大学数学学院，云南昆明 650500；3.云南省高校复杂系统建模及应用重点实验室，云南昆明 650500）

0 引言

泥石流是一种自然地质灾害，发生时对当地居民的人身和财产安全造成重大的危害。我国泥石流研究防治工作始于20 世纪60 年代，早期的泥石流研究依靠实地考察，结合实地数据分析得出防治方法。唐邦兴等[1-2]通过对我国山地泥石流的考察研究，绘制了我国早期1∶600 万泥石流灾害分布与危险区划图。陈宁生等[3-5]、崔鹏等[6-8]对四川、云南、新疆、西藏等特定地区泥石流的考察，提出了有效的泥石流判别方法和防治策略。随着遥感技术的发展，从数字高程模型（DEM）和遥感数据中获取与泥石流相关的地形地貌类因子和物源类因子弥补了实地考察的不便，通过提取沟谷的地貌因子、物源和降雨等因子，使用动力学模型FLO-2D、Debris2D、MatDEM 等模拟特定沟谷的泥石流动力学过程，能够较准确的研究该沟谷或此类沟谷泥石流形成机理并提出预警和防治方法。

为了能够大面积地研究和评估，将研究区域通过汇水面积或其他阈值划分为小流域单元，统计各单元的泥石流相关因子，通过模型学习泥石流流域因子蕴含的特征来完成其他流域的分类预测。李益敏等[9]提取云南怒江州小流域单元的流域距断裂带距离、岩性、流域水系密度等11 个因子，建立模型预测得到泥石流易发性分区，并绘制了怒江小流域易发性分布图。孙滨等[10]选取高程、坡度和坡向等9 个影响因子，建立东川地区的小流域评价模型进行泥石流易发性评价。赵岩[11]详细构建了地貌、物源和激发三类泥石流因子的数据库，选取流域10 min 降雨量、植被覆盖、6 h 均降雨量等8 个因子，使用线性回归的方式来对泥石流发生频率进行预测。统计学的方法能够在区域尺度上进行预测，但划分流域的方法不一，并且专家学者收集的影响因子和建模方法不同，导致预警评价结果不一，很难量化其评价预测的准确性。为了规避这些“主观”选择因子的过程，将卷积神经网络自动提取图像特征的方法用于地质灾害影像的研究成为了新热点。刘坤香等[12]提取怒江州泥石流沟谷的DEM 图像，使用改进的残差网络学习图像特征，并用其他沟谷的分类预测。杨小兵等[13]基于深度学习方法对泥石流堆积扇的多光谱遥感影像进行训练和识别，以此找寻更多的泥石流隐患点。神经网络学习需要大量的样本数据，而区域尺度上的灾害数据属于小样本问题，小样本学习的提出用于解决样本量少的学习情况。茹颖[14]在遥感影像较少的数据集上，使用小样本学习对遥感图像分类得到了比传统卷积神经网络更好的效果。张萌月[15]在遥感影像识别目标较少的情况下，使用小样本学习方法相对传统的卷积神经网络有绝对优势。

本文以沟谷为评价单元，提取怒江州遥感影像建立泥石流沟谷数据集，使用基于原型网络（Prototypes Net）[16]的小样本学习方法对发生泥石流的沟谷影像特征进行学习，并将模型用于其余大量沟谷的泥石流易发性评价，借助ArcGIS 软件将每条沟谷易发性结果和易发性分区进行可视化。

1 研究区域与实验数据

云南省怒江傈僳族自治州是全国发生泥石流最严重的州（市）之一，怒江州地处云南省西北部横断山脉纵向岭谷区，辖泸水市、福贡县、贡山县和兰坪县，总面积14 584.95 km2。州内高黎贡山、怒山山脉、云岭山脉近南北向延伸，澜沧江（湄公河）、怒江（萨尔温江）、独龙江（伊洛瓦底江）等深切峡谷平行南下，河谷谷坡坡度一般35°～45°、部分达60°～70°。按云南省各年减灾年鉴及网络中的不完全统计，截至2022 年，有相关文字记载的泥石流灾害中，福贡县至少发生69 起泥石流灾害事件、贡山县至少67 起、兰坪县至少24 起、泸水市至少42 起。1950—2021 年，福贡县各年平均降水量达到3 460 mm、贡山县1 632 mm、兰坪县1 550 mm、泸水市2 596 mm[17]。较高的降水量、高山峡谷的特殊地形及其他因素导致了怒江州泥石流灾害的频发。怒江州的地理位置和地形分布情况如图1 所示。

图1 怒江州地理位置和地形分布情况Fig.1 Geographical location and topographical distribution of Nujiang Prefecture

泥石流遥感数据集的建立主要使用谷歌地图（Google Earth）遥感影像，地图的高程数据由SRTM 提供，空间分辨率30 m，遥感影像则由Landsat，Quick-Bird，Spot 等遥感卫星提供。Google Earth 拥有高精度的时间和空间分辨率，城市区域空间分辨率可达1m，农村地区可达2.5～5 m，边远山区可达15 m[18]。数据的提取和处理如图2 所示，利用数字高程模型（DEM）获取每条沟谷的掩膜信息，按照掩膜批量提取Google Earth 遥感图像，为了保留沟谷的大小特征，将提取的遥感图像按照DEM 比例进行重采样，并按照最大DEM沟谷图像的尺寸进行填充处理。

图2 Google Earth 遥感数据的提取与处理Fig.2 Extraction and processing of Google Earth remote sensing data

所使用的Google Earth 遥感和DEM 数据的信息如表1 所示。

表1 Google Earth 遥感和DEM 数据信息Table 1 Information on Google Earth remote sensing and DEM data

泥石流遥感数据集包含正样本、负样本和待评价沟谷，正样本的确定是从《云南减灾年鉴》[19]、相关文献和网络报道中获得发生泥石流山谷和村庄的精确地点，依据坐标在Google Earth 中与发生泥石流后最近的时间分辨率进行沟谷图像提取。负样本沟谷的确定依据为：在历史记录中不能检索到发生泥石流信息，影像中没有明显的冲积痕迹和堆积扇的形成且周围村庄密集。最终收集到泥石流沟谷50 条作为正样本，未发生泥石流沟谷42 条作为负样本，待评价沟谷样本600 条。

2 基于原型网络的沟谷泥石流灾害易发性评价方法

2.1 沟谷泥石流易发性评价流程

基于原型网络的沟谷泥石流易发性评价流程如图3 所示，流程包括数据收集、沟谷提取、流域划分、数据集划分、数据集预分类、模型训练和预测、易发性评价和易发性分区图，详细流程如下。

图3 基于原型网络的沟谷泥石流易发性评价流程Fig.3 Chartflow of susceptibility assessment of debris flow based on prototype networks

首先收集研究区域的DEM、Google Earth 遥感影像和历史泥石流数据，按上一节中的数据提取流程，以沟谷流域为评价单元，提取得到Google Earth 沟谷泥石流遥感影像数据集和沟谷流域分区。其次，考虑形态各异的沟谷潜在危险性存在特征差异，直接将沟谷分为正负样本两类训练，在忽略特征差异的同时网络模型也不稳定，按照影响泥石流的重要因子流域面积[20]，将正负样本按流域面积大小各预分为小、中、大3 类，0、1、2 属于正样本，3、4、5 属于负样本，如图3 中数据集预分类标签所示。以元学习的形式组织正负样本数据并输入到原型网络模型中进行训练，将待评价沟谷样本输入到已训练的模型中，与6 个沟谷分类进行相似度量，计算每个评价样本的泥石流易发性指数。最后将易发性指数分区，得到沟谷的泥石流易发性评价等级，在流域单元划分中可视化评价结果，得到易发性评价分区图。

从流程图中可以看出，与基于泥石流因子进行评价的流程相比，基于原型网络的沟谷泥石流易发性评价方法直接使用沟谷影像来学习泥石流沟谷的特征，减少了从小流域计算众多泥石流因子和泥石流因子筛选的繁琐过程，避免了选择不同因子带来的评价差异问题。原型网络在保证性能的前提下，能够快速高效的进行泥石流易发性的评价。

2.2 原型网络模型与泥石流易发性评价算法

原型网络是小样本学习中经典的元学习和度量学习方法[21]，网络模型如图4 所示。原型网络以元学习的方式组织泥石流数据集进行训练，通过特征提取器（feature extractor）提取支持集（support set）每一类沟谷的特征，计算每类沟谷特征均值作为该类的原型中心（Ck），将查询集（query set）中待测沟谷的特征（Q）与原型中心进行距离度量来完成待测沟谷的分类和预测。同时，按照度量距离计算待测沟谷从属类别的概率，并计算待测沟谷的泥石流易发性指数，从而得到待测样本的泥石流易发性评价等级。

图4 原型网络及分类模型Fig.4 Prototype network classification model

2.2.1 基于元学习方法的原型中心计算

元学习分为元训练和元测试。元训练将有限的训练样本分为不同的子任务进行训练，在小样本学习中称为C-wayK-shot 问题。具体为：在N类训练样本中，每类有少量样本，学习时每个子任务选取C类（C≤N），每类中选取K个样本，共C×K个样本组成支持集（support set），即子任务的训练集。每类选取K个样本后，从每类剩余样本中各选取Q个样本（K+Q小于等于类内样本总数），共C×Q个样本作为查询集（query set），即子任务的测试集。不断从总类N中抽取C-wayKshot 组成支持集进行训练，充分地利用了训练样本数据，使得模型在新任务时表现得更好。元测试则使用相同的方法来测试模型性能。

在训练和测试泥石流数据时，原型网络将VGG[22]、GoogleNet[23]、ResNet[24]等不带分类器的卷积网络作为特征提取器，计算支持集中每一类沟谷的原型中心作为该类沟谷的代表，原型中心计算如式（1）所示：

式中：Ck——每一类沟谷的原型中心；

Sk——训练批次中每一类沟谷的样本数；

(xi,yi)——某一样本和对应的标签；

f∅——特征提取器；

2.2.2 样本类别概率计算

如果将普通的神经网络分类器训练在小样本学习的任务中，几乎都是过拟合。基于非参数度量的方法（例如：最近邻、k最近邻、K-means、距离度量等）不需要优化参数，可以在元学习框架下构建端到端的小样本分类器。原型网络使用余弦或欧氏距离的距离度量方式作为分类器，通过距离度量沟谷样本的特征与每类原型中心的距离大小来完成分类，并通过式（2）计算沟谷样本属于某一类沟谷的概率：

式中：p∅——某一样本属于真实类别k的概率；

∅——提取器的参数。

d——距离度量函数；

Q——样本通过f∅提取的特征；

k′——同一批训练的每一类别。

在训练优化阶段，通过将式（3）中的目标函数J(∅)反向传播和随机梯度下降来对网络模型参数 ∅进行优化。

2.2.3 沟谷泥石流易发性等级评价

在进行沟谷的泥石流易发性评价时，通过计算沟谷的泥石流易发性指数来判定沟谷易发性等级。易发性指数通过所属类别概率来计算，文中原型网络使用欧氏距离来度量评价样本特征向量与类原型特征向量之间的相似程度，欧式距离越小，相似程度越大，使用式（2）计算所属类别的概率也就越大。

评价沟谷所属正样本中0、1、2 类其中一类沟谷的概率越高，说明样本泥石流易发性越高。使用式（4）将所属概率转化为易发性指数。具体方法为：将正样本每类所属概率之和与负样本每类所属概率之和作差，差值大于等于0，表示评价样本与正样本相似程度更大，取正样本中的最大值作为易发性指数，值越大易发性越大。反之，差值小于0，取负样本中最大值的相反数为易发性指数，值越小易发性越小。将沟谷的易发性指数按大小进行分区，即可得到沟谷的易发性等级。

式中：Ii——验证样本i的易发性指数；

m、k——正样本、负样本的类别数，文中m=k=3。

3 实验结果与分析

3.1 分类性能

本实验实现的硬件为CPU：Intel Xeon E5-2 650 v3、内存：128GB、GPU：NVIDIA GeForce RTX 3 090。软件为Ubuntu18.06、Pytorch1.9.0、python3.8、CUDA11.1。实验时设置训练输入为6-way 10-shot，即每轮选6 类，每类选10 张训练，测试时输入6-way 4shot 1query，即每轮测试选取6 类，每类4 张计算原型，1 张作为测试。选取9 个网络模型结构，包括VGG、GoogleNet、ShuffleNetV2[25]、MobileNetV2[26]、Conv4[27]、ResNet12[27]、ResNet18[24]、DenseNet[28]、Rir[29]作为原型网络的特征提取器，分类结果如表2 所示。

表2 多个特征提取器的分类性能Table 2 Classification performance of multiple feature extractors

表2 中，2 分类为只考虑正负样本的分类正确率，可以看出，提取沟谷特征时，Conv4 在原型网络中的6 分类和2 分类正确率优于其他模型。将正负样本沟谷输入此模型进行测试，得到6 分类结果的混淆矩阵和正负2 分类的混淆矩阵如表3 和表4 所示。

表3 6 分类测试混淆矩阵Table 3 Confusion matrix of six-classification

表4 正负2 分类测试混淆矩阵Table 4 Confusion matrix for positive and negative binaryclassification test

使用分类指标对分类结果进行衡量，计算6 分类的卡帕系数（Kappa）和2 分类的准确率（Accuracy）、精确率（Precision）、召回率（Recall）、F1 Score，得到表5 指标表。

表5 2 分类和6 分类指标表Table 5 Binary-classification and six-classification indicator table

表2—4 表明，Conv4 表现出了良好的分类性能，作为模型的特征提取器能够带来更好的泥石流易发评价结果。

3.2 评价结果及分析

保持同样的模型参数，将600 张怒江州待评价沟谷及正负样本输入网络，计算每一个沟谷样本属于6 类沟谷的概率，按式（4）计算每个样本的易发性指数Ii。使用自然断点法，按照易发性指数将易发区分为5 个等级：低易发区[-1,-0.760 7]、较低易发区(-0.760 7,0]、易发区(0,0.603 6]、较高易发区(0.603 6,0.828 8]、高易发区(0.828 8,1]。最后，在ArcGIS 中使用DEM 进行易发性分区可视化，得到如图5 所示分区图。

图5 怒江州泥石流易发性分区图Fig.5 Susceptibility zoning map of debris flow in Nujiang Prefecture

从图5 中可以看出，预测为高易发的沟谷主要分布在贡山县、福贡县、泸水市中，兰坪县较少。易发性沟谷分区中，低易发区249 条、较低易发区96 条、易发区100 条、较高易发区124 条、高易发区123 条。使用ArcGIS 计算沟谷的高程差、主沟长度、沟谷投影面积和坡降比。计算时，高程差为ArcGIS 中使用沟谷的最高海拔减去沟口海拔得出，主沟长度为ArcGIS 中使用流量计算出最长的河流长度，沟谷面积为使用CGCS-2000 投影坐标得到的栅格面积，坡降比即沟床比降，为高程差与主沟长的比值，反映沟谷的整体陡度。按照易发性分区将高程差、主沟长度、沟谷面积和坡降比分布可视化，如图6 所示，作为对比，正样本分布如图6 中虚线空心箱线图所示。

图6 怒江州泥石流易发性分区中沟谷高程差、沟长、面积、坡降比的分布情况Fig.6 Distribution of the valley elevation difference,valley length,area,and slope ratio in the debris flow susceptibility zone in Nujiang prefecture

从分布情况中可知，在预测的沟谷中，低易发区、较低易发区、易发区、较高易发区在高程差、沟长、面积的分布上呈现下降趋势，坡降比呈上升趋势。易发区和较高易发区分布在高程差小、沟长短、面积小和坡降大的沟谷中。高易发区沟谷分布则较广，面积、高程差、沟长和坡降比从小到大都有高易发沟谷，这和正样本沟谷的分布有着相似的特点。

李益敏[9]使用泥石流因子评价方法，以5 km2的汇水面积为阈值，将怒江州划分为1 414 个小流域单元，选取流域单元的距断裂带距离、岩性、Melton 指数[30]、流域延伸率、流域高差率、河流弯曲系数、流域水系密度、平均植被覆盖度、年均降水量、距道路距离、距居民点距离共 11 个评价因子，采用确定性系数模型 CF和多因子叠加权重确定法，对怒江州小流域单元进行泥石流易发性评价分区。与本文易发性分区对比，存在着相同点：预测为高易发及较高易发的沟谷（其文中为小流域）都集中在澜沧江以西，兰坪县较少，这与泥石流历史分布都较为符合。预测为高易发及较高易发沟谷（流域）数分别占全部沟谷（流域）的35.69%和34.72%，非常相近。

3.3 典型沟谷分析

在《云南省减灾年鉴》中检索到4 条发生了人员伤亡或多次发生泥石流的沟谷，地理位置如图7 中标签A、B、C、D 处所示，编号A 为贡山县普拉底乡东月谷，于2010 年8 月18 日凌晨1 时30 分暴发特大山洪泥石流灾害，造成11 212 人受灾，死亡96 人，失踪53 人，重伤9 人，轻伤30 人，转移安置灾民3 282 人，直接经济损失1.4 亿元。编号B 为贡山县普拉底乡咪谷河，于2010年7 月26 日凌晨发生泥石流灾害，造成3 人死亡、8 人失踪、11 人受轻伤，直接经济损失200 万元。编号C 为贡山县独龙江乡巴坡村沟谷，于2010 年4 月4 日12 时发生泥石流灾害，造成1 人失踪，公路与通讯中断，巴坡村2014 年2 月20 日、2017 年8 月12 日都有发生泥石流的记录。编号D 为兰坪县金顶镇七联村练登大沟，于2010 年9 月19 日发生泥石流，对下游基础设施造成破坏，直接经济损失360 万元，又于2019 年8 月7 日发生泥石流。4 条沟谷的遥感图像和沟口泥石流堆积情况如图7 中所示。

4 条泥石流沟谷的所属概率和易发性指数如表6所示。

表6 4 条泥石流沟谷原型网络计算的所属概率和易发性指数Table 6 Probability and susceptibility index of four debris flow valleys with prototype networks

在预测中东月谷、咪谷河和练登大沟与1 类泥石流沟谷相似度最高，巴坡村与0 类泥石流沟谷相似度最高，东月谷、咪谷河和巴坡村三条沟谷的易发性指数都属于高易发性区间，这说明三条沟谷具有非常高的泥石流易发性。练登大沟易发性指数在易发区间，具有较高的泥石流易发性。

基于因子统计分析是较为成熟的分析方法，文中收集了沟谷相关的高程差、主沟长度、沟谷面积、沟谷周长、沟谷平均坡度、坡降比、Melton 指数、切割密度和圆状率共9 个地形因子。其中，Melton 指数[30]由Melton M A提出，其值为流域高程差与面积平方根的比值，反映了集水区的动力学及其对泥石流的敏感性。切割密度为主沟长度与流域面积的比值，反映集水区内岩石的抗风化能力和地貌发育情况[31]。圆状率由沟谷周长与主沟长度比值所得，影响汇流时间。对9 个因子使用随机森林的分类方法进行训练，模型正负样本分类准确率为61.11%，对4 条沟谷进行预测，预测结果如表7所示。

表7 4 条泥石流沟谷因子分析方法计算的所属概率和易发性指数Table 7 Probability and susceptibility index of four debris flow valleys with factor analysis methods

可以看出，虽然使用因子分析和原型网络对4 条沟的预测类别不同，但除了巴坡村预测错误，其他3 条沟仍被预测为正样本。同样的，与李益敏[9]基于因子分析的易发性分区相比，4 条沟谷在其评测中的易发性结果为：东月各为低易发、咪谷河为高易发、巴坡村为低易发、练登大沟为极低易发。对于四条已经发生泥石流的沟谷来说，特别是东月各沟谷发生特大泥石流灾害的特点，4 条沟谷在本文原型网络评价中为高易发和较高易发，这与历史情况更为接近。

泥石流的发生是地貌条件、物质条件和激发条件共同作用的结果[11]，选取若干小流域易发性评价实验[32-34]中认为重要的地貌条件和物质条件因子，对上面4 条泥石流沟谷进行分析，分析如表8 所示[35-36]。

表8 4 条泥石流沟谷的地貌条件和物质条件因子分析Table 8 Factors analysis of geomorphological and material conditions of four debris flow valleys

在地形地貌和物质条件对泥石流形成的影响中，流域面积大汇水量多，主沟长、坡度大、高程差大、坡降比大、Melton 指数大都会为泥石流的冲积提供动能。岩性对泥石流发生的影响中，片岩、混合岩、砂岩都容易风化产生松散物质，花岗岩、板岩、千枚岩、石灰岩、其他碳酸盐岩虽然比片岩、混合岩更硬，但是对泥石流也有较高的敏感性[37-40]。灌林和草地的含沙量比林地多，水渗透能力也比林地强，能够提供更多的松散物质和更强的雨水侵蚀能力[41-42]。淋溶土和松软薄层土也提供了松散物质[43]。东月谷、咪谷河、练登大沟3 条沟谷中，疏林地、草地、松软薄土和淋溶土、砂岩和片岩等为泥石流提供了大量松散物质，练登大沟半个坡面为人工矿区，对植被破坏大，大量裸露的松散泥土，容易被雨水冲刷和搬运。其次，3 条沟谷都有较长的主沟长度、较大的汇水面积、高程差、平均坡度和Melton 指数，汇集大量雨水且促成强大的冲积能量，从图7 中的堆积扇可以看出，冲积泥沙量巨大，阻断河流形成堰塞湖、毁坏公路和农田、摧毁房屋，造成巨大的经济损失。以人员伤亡和经济损失最大的东月谷为例，暴发泥石流时冲出体积约为 60×104m3[44]，破坏力极强。巴坡村虽然沟长和面积较小，但同样有足够的松散物质，并且坡度极大，降雨量达到一定程度，极容易形成小型泥石流或者滑坡，每当松散物质得到补充，只要一定的雨水条件激发，就会周期性的发生泥石流。

4 讨论

文中考虑到形态各异的沟谷潜在的特征差异及模型的稳定性，按面积大小对泥石流沟谷进行预分类，再进行原型网络的训练，可能存在更为合理的预分类，让神经网络能更好地区分沟谷潜在的图像特征从而使得评价更为准确。其次，能够提取到沟谷图像的还有高分卫星遥感和数字高程模型，这2 种图像对泥石流特征是否有着更好的表征能力，或者使用它们的组合能够带来更好的评测效果，是今后研究工作的一个重点。最后，以沟谷为单元进行评价，从历史信息中能够确定发生泥石流的沟谷数量极其有限，特别是极小部分多次发生泥石流的高频沟谷，网络并没有学习到这种高频特征。

李益敏等[9]使用因子评价的方法中，评价为高易发性的流域在贡山县、福贡县和泸水市中，但主要集中在贡山县，而本文中贡山县稍多，同时在贡山县、福贡县和泸水市中的分布比其更加“均匀”，更加符合历史泥石流灾害的分布规律。但是，评价的准确度需要从后续泥石流灾害事件中得到验证。而且，小样本深度学习的方法依赖于数据和模型，在后期的工作中，继续收集更多的泥石流沟谷图像，并不断改进模型，提高模型的精度，以此来减少模型性能带来的误差，提高预测精度。