一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法

专利名称：一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法

专利类型：发明专利

专利申请号：CN202410813851.1

专利申请（专利权）人：中国科学院、水利部成都山地灾害与环境研究所
权利人地址：四川省成都市天府新区群贤南街189号

专利发明（设计）人：蒋祯妮,王姣,崔鹏,张国涛,刘杰,杨治纬

专利摘要：本发明涉及一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，属于图像数据处理领域，包括：数据预处理步骤；模型参数建立及处理步骤：建立包括沟道数量、堆积扇面积大小、最大冲击力度和灾害发生频率的模型参数，并进行沟道数量识别、堆积扇面积大小计算、最大冲击力度计算和灾害发生频率获取；模型构建及评估结果指标化步骤：将四个参数的结果统一标准化处理后作为主成分分析法的输入变量，求取对模型整体贡献量大的主成分数量，获取各主成分在模型中的权重，构建模型结果。本发明可自动化、重复进行大区域大范围的灾害分析，缩减了人工勘察和现场监测等低效方法的分析时间和成本，结果生成稳定、标准化、可迭代优化而适用于不同环境和场景。

主权利要求：
1.一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述评估方法包括：数据预处理步骤；
模型参数建立及处理步骤：建立包括沟道数量、堆积扇面积大小、最大冲击力度和灾害发生频率的模型参数，基于单个流域内的沟道纵剖面特征定位沟口及堆积扇终点从而获取堆积扇横截中线及通过捕捉其谷信号推测沟道数量，通过堆积扇沟口和终点位置获得堆积扇半径及利用已识别堆积扇边界沟道计算堆积扇中线宽度及扩散角从而计算堆积扇面积，通过野外考察所得面积和高的关系计算巨砾高度和体积并结合所得堆积扇平均坡度计算灾害最大冲击力度，根据不同时期的沟道识别数量获取时期内灾害发生频数及推算灾害发生频率；
模型构建及评估结果指标化步骤：将四个参数的结果统一标准化处理后作为主成分分析法的输入变量，求取对模型整体贡献量大的主成分数量，获取各主成分在模型中的权重，构建模型，运行模型，获得泥石流活动评估结果。
2.根据权利要求1所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述基于单个流域内的沟道纵剖面特征定位沟口及堆积扇终点从而获取堆积扇横截中线及通过捕捉其谷信号推测沟道数量具体包括以下内容：获取沟道纵剖面数据：将各流域DEM经纬坐标数据转化为平面图像像素坐标格式的数组，读取该流域内的河流矢量数据，通过河网的分级筛选主河道数据，根据河段对应的上下节点标识属性和数据拓扑连接规律从上游至下游交叉匹配数据形成连续河流数据，并将河流数据坐标转化为与DEM相同的像素坐标，获取河流数据点对应DEM数据中的高程值及根据点经纬度坐标计算各点间实际地理距离，从而得到沟道纵剖面数据；
判断和优化沟口及堆积扇终点位置：根据各流域内沟道纵剖面数据计算邻点间坡度、坡度变化值、沟道平均坡度变化和坡度变化标准差，以预测最接近现实沟口的结果，并根据CI阈值方法确定最终所得流域的沟口位置和堆积扇终点；
提取堆积扇中线剖面：根据所得沟口和堆积扇终点坐标连接成线，计算其平均坡度和中点坐标，以中点为起点作正方向量垂线，即为流域中该堆积扇的中线剖面；
识别堆积扇横截中线上的沟道数量：识别堆积扇中线剖面中的局部最大值即主要峰信号和谷信号，保存预估谷信号和峰信号数据高程及剖面位置信息，利用高差变化筛除异常值，并设定沟道点高程最高不超过流域堆积扇以上地形高程，得到调整沟道数量，并结合实际解译的沟道数量数据，得到调整沟道数量和实际沟道数量的关系。
3.根据权利要求2所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述根据CI阈值方法确定最终所得流域的沟口位置和堆积扇终点包括：使用CI阈值寻找下游至上游地形突变的第一个正值特征点，通过方程
确定沟口坐标，其中， 、 和s分别代表坡度变化、平均
坡度、标准差，t表示临界t‑value，n表示数据大小，Point(x,y)表示符合不等式要求的最低点的点坐标，即沟口坐标；
根据最优的CI阈值得到最终所得流域的沟口位置，堆积扇终点则为流域内河流逆流方向上坡度变化急速转正的点，但如果生成流域边界时流域底部边界刚好与堆积扇边界重合，则提取流域底部边界与河流的交点作为堆积扇终点。
4.根据权利要求2所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述通过堆积扇沟口和终点位置获得堆积扇半径及利用已识别堆积扇边界沟道计算堆积扇中线宽度及扩散角从而计算堆积扇面积具体包括以下内容：计算堆积扇半径：根据沟口和堆积扇终点位置及高程，结合三角形勾股定理计算堆积扇半径；
计算堆积扇扩散角预测值：根据优化沟道的双向极点，计算堆积扇中线长，并通过公式计算堆积扇扩散角预测值 ，其中，arctan为反正切函数，Wmid‑fan为堆积扇中线长，R为堆积扇半径；
校正堆积扇扩散角：通过人工解译样本得到的堆积扇扩散角与计算得到的堆积扇扩散角预测值之间的关系，对堆积扇扩散角预测值进行校正得到校正堆积扇扩散角α；
推算堆积扇面积：根据堆积扇半径R和校正堆积扇扩散角α计算得到堆积扇面积为。
5.根据权利要求1所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述通过野外考察所得面积和高的关系计算巨砾高度和体积并结合所得堆积扇平均坡度计算灾害最大冲击力度具体包括以下内容：分割堆积扇上的巨砾：在确保最佳分割性能的前提下，采样聚焦于小区域及扇体位置居中的高分辨率影像或者无人机影像进行颗粒检测分割；
筛选合格图像数据：将输出结果依像素数量高低进行排序，计算特征物面积与图像总面积比筛除面积大于等于第一预设值且预测交并比和稳定系数大于等于第二预设值的特征物；
计算所用图像单个像素实际尺寸：根据相似三角形定理，依据图上距离和拍照设备参数计算图上单个像素代表的实际尺寸；
计算巨砾实际尺寸：根据图上像素的实际面积以及特征物的像素数量，得到特征物即巨砾的实际顶面面积；
建立野外巨砾面积与高的关系：根据野外考察所得巨砾长、宽和高的数据，得到巨砾长宽所得面积与高的关系为 ，其中， 表示巨砾的高， 表示巨砾的长‑宽所得面积；
推算巨砾体积：根据面积与高的乘积关系得到巨砾的体积为 ；
计算最大冲击力度：根据最大冲击力度计算公式得到不同流域的泥石流最大冲击力度。
6.根据权利要求2所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述根据不同时期的沟道识别数量获取时期内灾害发生频数及推算灾害发生频率具体包括以下内容：识别对照DEM中出现的沟道数量；
计算沟道不同时期沟道数量差：将对照DEM中的各流域的沟道数量与基准DEM中各流域的沟道数量相减得到不同时期沟道数量差；
计算不同DEM的时间差；
将沟道数量差除以时间差得到灾害发生频率。
7.根据权利要求1所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述模型构建及评估结果指标化步骤具体包括以下内容：将各参数变量进行批量归一化，在确保所有变量享有平等的可贡献程度的前提下，保留数据的分布特点；
将高维数据投影到低维空间中，实现PCA降维运算，并对标准化的数据进行拟合和转换，得到主成分得分；
计算每个主成分对数据总方差的独立贡献即解释方差比，累积解释方差比以选出构建模型所需的主成分；
生成从1到主成分个数的主成分序列，利用主成分特征向量表示原始变量即四个模型参数变量在主成分空间中的贡献，获取PCA分析的载荷矩阵；
根据协方差函数计算各变量间的协方差，并从协方差矩阵中提取主成分间的协方差值；
计算提取的主成分协方差值的总和，按照各协方差值占总和的比例确定各成分的归一化权重；
模型运算，获得泥石流灾害活动评估结果；
按需扩缩指标数值范围。
8.根据权利要求1‑7中任意一项所述的一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，其特征在于：所述数据预处理步骤具体包括以下内容：生成DEM数据：用ArcGIS的镶嵌工具分别镶嵌基准DEM数据和对照DEM数据，并向下重采样使镶嵌后的对照DEM与基准DEM分辨率一致，从而减少后续灾害频率计算中因数据分辨力差异所带来误差的影响；
划分流域并生成河网矢量数据：在ArcGIS中输入基准DEM数据，运行水文分析工具进行填挖、流向、集水区、汇流分析、河网链接、斯特拉勒河网分级、河网矢量化、流域划分、流域矢量化一系列操作，再叠加地形和影像数据，检查输出的河网和流域矢量数据的质量，并人工校正，生成出流域边界分明、河网分布准确的流域和河网数据；
制作各流域边界、DEM和河流数据：将所有流域根据ID分割和根据ID赋名，将以ID为名的流域数据制作成掩膜，输入基准DEM数据分别生成各流域DEM数据，并裁剪各流域内的河流数据；
野外勘测数据采样：使用卷尺测量堆积扇上巨砾的长宽高信息；
生成人工解译灾害数据编目：筛选出低云量的遥感影像数据，在已有流域信息的基础上，结合野外验证和遥感影像上的形态、颜色、纹理信息，推断不同含水率及颗粒组成下的泥石流类别，并赋值各流域灾害类别代号；
人工解译流域沟道数量：根据遥感影像数据，目视解译随机样本流域内的沟道数目；
人工解译流域内堆积扇扩散角：根据遥感影像数据，目视解译随机样本流域内的堆积扇扩散角；
利用照相机或者无人机拍摄包括剖面、地貌、灾害影响范围和周边环境的照片。 说明书 : 一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法技术领域[0001] 本发明涉及图像数据处理领域，尤其涉及一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法。背景技术[0002] 泥石流是一种常见的地质灾害，广泛发生于全球的山地、丘陵地区；作为一种固液混合流，它的破坏力极强，不仅会毁坏房屋和基础设施、产生巨大经济损失，还会对人类生命和财产安全构成严重威胁。泥石流的发生常受暴雨、高温等极端气候事件影响，具有流速快、冲击力大、突发性和随机性的特点，导致预测难、防治难和恢复重建困难。近年来，随着气候变化的加剧和人类活动的干扰，泥石流的频率增多、规模增大，尤其在旱区表现尤为突出，且群发性明显，对社会稳定、山区可持续发展构成极大挑战。由于泥石流常发育于偏远山区、环境险峻，其爆发不确定性大，导致数据采集困难、安全风险较高，泥石流规模和频率难以获取。但掌握这些与泥石流活动性相关的信息是泥石流防治和管理的基础。[0003] 目前，评估泥石流规模和频率的传统方法包括野外监测和堆积剖面测量，通过实地架设流量、流速监测设备和视频监控获取泥石流相关参数，或通过野外剖面测量和拍照获取泥石流规模信息。现行基于监测数据或堆积剖面测量泥石流规模、频率的方法，由于1）监测网络设备、安装、维护费用高，且堆积扇剖面测量需实地调查，费时、耗力；2）堆积扇活动性推测对数据量需求大、对数据分辨率要求高；3）均以单沟尺度应用为主；4）灾害发育时空和严重程度难以预测；5）受剖面采样位置、取样深度、监测记录时长影响。因此存在以下主要缺陷：1）设备、人力、时间成本高；2）应用范围受监测网络部署和现场勘察范围限制；3）难以或较少在大范围、垮流域情况下同步应用；4）难以满足短时、标准化、多样本、高效评估泥石流活动的需求、存在设备损毁或安全调查隐患；5）泥石流活动频率和规模分析结果准14确性和完整性难以保证。此外，测年法中树轮及 C方法取样及分析过程均较复杂、分析周期14长、旱区木本植物稀缺且土壤有机质含量少，C分析结果常为百年或千年尺度。导致这些方14法短时间内难以获得结果、耗财耗力、在旱区环境下适用性受限，且 C测年法难以推断近现代灾害活动频率，结果时间分辨率较低。数值模拟方式需要精确复刻孕灾环境，筛选致灾主控因素，对不同类型灾害的致灾机制和形成机理有全面了解，以还原泥石流形成过程并演算致灾规模，因此其结果受模型制约并与实际情况存在一定偏差。依赖遥感影像目视解译的方式，虽然适用范围大、费用较低，但由于分辨率参差、不同人解译标准和思维存在差异、大范围解译常需要多人合作，因此结果精度难以保证、受解译人员专业性影响较大、人工解译时间成本高。[0004] 然而，对比传统的监测法和剖面测量法，测年法和数值模拟方式只提供了新型研究方法，在推算泥石流成灾范围和活动期次上并未实现效率上的提升，适用环境和范围也未明显扩大；而遥感解译法尽管扩大了适用范围，但效率仍然较低，仍难以快速高效评估频繁发生于偏远山区的泥石流灾害。发明内容[0005] 本发明的目的在于克服现有技术的缺点，提供了一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，解决了现有技术存在的不足。[0006] 本发明的目的通过以下技术方案来实现：一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，所述评估方法包括：[0007] 数据预处理步骤；[0008] 模型参数建立及处理步骤：建立包括沟道数量、堆积扇面积大小、最大冲击力度和灾害发生频率的模型参数，基于单个流域内的沟道纵剖面特征定位沟口及堆积扇终点从而获取堆积扇横截中线及通过捕捉其谷信号推测沟道数量，通过堆积扇沟口和终点位置获得堆积扇半径及利用已识别堆积扇边界沟道计算堆积扇中线宽度及扩散角从而计算堆积扇面积，通过野外考察所得面积和高的关系计算巨砾高度和体积并结合所得堆积扇平均坡度计算灾害最大冲击力度，根据不同时期的沟道识别数量获取时期内灾害发生频数及推算灾害发生频率；[0009] 模型构建及评估结果指标化步骤：将四个参数的结果统一标准化处理后作为主成分分析法的输入变量，求取对模型整体贡献量大的主成分数量，获取各主成分在模型中的权重，构建模型结果，进行泥石流活动评估。[0010] 所述基于单个流域内的沟道纵剖面特征定位沟口及堆积扇终点从而获取堆积扇横截中线及通过捕捉其谷信号推测沟道数量具体包括以下内容：[0011] 获取沟道纵剖面数据：将各流域DEM经纬坐标数据转化为平面图像像素坐标格式的数组，读取该流域内的河流矢量数据，通过河网的分级筛选主河道数据，根据河段对应的上下节点标识属性和数据拓扑连接规律从上游至下游交叉匹配数据形成连续河流数据，并将河流数据坐标转化为与DEM相同的像素坐标，获取河流数据点对应DEM数据中的高程值及根据点经纬度坐标计算各点间实际地理距离，从而得到沟道纵剖面数据；[0012] 判断和优化沟口及堆积扇终点位置：根据各流域内沟道纵剖面数据计算邻点间坡度、坡度变化值、沟道平均坡度变化和坡度变化标准差，以预测最接近现实沟口的结果，并根据CI阈值方法确定最终所得流域的沟口位置和堆积扇终点；[0013] 提取堆积扇中线剖面：根据所得沟口和堆积扇终点坐标连接成线，计算其平均坡度和中点坐标，以中点为起点作正方向量垂线，即为流域中该堆积扇的中线剖面；[0014] 识别堆积扇横截中线上的沟道数量：识别堆积扇中线剖面中的局部最大值即主要峰信号和谷信号，保存预估谷信号和峰信号数据高程及剖面位置信息，利用高差变化筛除异常值，并设定沟道点高程最高不超过流域堆积扇以上地形高程，得到调整沟道数量，并结合实际解译的沟道数量数据，得到调整沟道数量和实际沟道数量的关系。[0015] 所述根据CI阈值方法确定最终所得流域的沟口位置和堆积扇终点包括：[0016] 使用CI阈值寻找下游至上游地形突变的第一个正值特征点，通过方程确定沟口坐标，其中， 、 和s分别代表坡度变化、平均坡度、标准差，t表示临界t‑value，n表示数据大小，Point(x,y)表示符合不等式要求的最低点的点坐标，即沟口坐标；[0017] 根据最优的CI阈值得到最终所得流域的沟口位置，堆积扇终点则为流域内河流逆流方向上坡度变化急速转正的点，但如果生成流域边界时流域底部边界刚好与堆积扇边界重合，则提取流域底部边界与河流的交点作为堆积扇终点。[0018] 所述通过堆积扇沟口和终点位置获得堆积扇半径及利用已识别堆积扇边界沟道计算堆积扇中线宽度及扩散角从而计算堆积扇面积具体包括以下内容：[0019] 计算堆积扇半径：根据沟口和堆积扇终点位置及高程，结合三角形勾股定理计算堆积扇半径；[0020] 计算堆积扇扩散角预测值：根据优化沟道的双向极点，计算堆积扇中线长，并通过公式 计算堆积扇扩散角预测值 ，其中，arctan为反正切函数，Wmid‑fan为堆积扇中线长，R为堆积扇半径；[0021] 校正堆积扇扩散角：通过人工解译样本得到的堆积扇扩散角与计算得到的堆积扇扩散角预测值之间的关系，对堆积扇扩散角预测值进行校正得到校正堆积扇扩散角α；[0022] 推算堆积扇面积：根据堆积扇半径R和校正堆积扇扩散角α计算得到堆积扇面积为。[0023] 所述通过野外考察所得面积和高的关系计算巨砾高度和体积并结合所得堆积扇平均坡度计算灾害最大冲击力度具体包括以下内容：[0024] 分割堆积扇上的巨砾：在确保最佳分割性能的前提下，采样聚焦于小区域及扇体位置居中的高分辨率影像或者无人机影像进行颗粒检测分割；[0025] 筛选合格图像数据：将输出结果依像素数量高低进行排序，计算特征物面积与图像总面积比筛除面积大于等于第一预设值且预测交并比和稳定系数大于等于第二预设值的特征物；[0026] 计算所用图像单个像素实际尺寸：根据相似三角形定理，依据图上距离和拍照设备参数计算图上单个像素代表的实际尺寸；[0027] 计算巨砾实际尺寸：根据图上像素的实际面积以及特征物的像素数量，得到特征物即巨砾的实际顶面面积；[0028] 建立野外巨砾面积与高的关系：根据野外考察所得巨砾长、宽和高的数据，得到巨砾长宽所得面积与高的关系为 ，其中， 表示巨砾的高，表示巨砾的长‑宽所得面积；[0029] 推算巨砾体积：根据面积与高的乘积关系得到巨砾的体积为 ；[0030] 计算最大冲击力度：根据最大冲击力度计算公式得到不同流域的泥石流最大冲击力度。[0031] 所述根据不同时期的沟道识别数量获取时期内灾害发生频数及推算灾害发生频率具体包括以下内容：[0032] 识别对照DEM中出现的沟道数量；[0033] 计算沟道不同时期沟道数量差：将对照DEM中的各流域的沟道数量与基准DEM中各流域的沟道数量相减得到不同时期沟道数量差；[0034] 计算不同DEM的时间差；[0035] 将沟道数量差除以时间差得到灾害发生频率。[0036] 所述模型构建及指标化步骤具体包括以下内容：[0037] 将各参数变量进行批量归一化，在确保所有变量享有平等的可贡献程度的前提下，保留数据的分布特点；[0038] 将高维数据投影到低维空间中，实现PCA降维运算，并对标准化的数据进行拟合和转换，得到主成分得分；[0039] 计算每个主成分对数据总方差的独立贡献即解释方差比，累积解释方差比以选出构建模型所需的主成分；[0040] 生成从1到主成分个数的主成分序列，利用主成分特征向量表示原始变量即四个模型参数变量在主成分空间中的贡献，获取PCA分析的载荷矩阵；[0041] 根据协方差函数计算各变量间的协方差，并从协方差矩阵中提取主成分间的协方差值；[0042] 计算提取的主成分协方差值的总和，按照各协方差值占总和的比例确定各成分的归一化权重；[0043] 模型运算，获得泥石流灾害活动评估结果；[0044] 按需扩缩指标数值范围。[0045] 所述数据预处理步骤具体包括以下内容：[0046] 生成DEM数据：用ArcGIS的镶嵌工具分别镶嵌基准DEM数据和对照DEM数据，并向下重采样使镶嵌后的对照DEM与基准DEM分辨率一致，从而减少后续灾害频率计算中因数据分辨力差异所带来误差的影响；[0047] 划分流域并生成河网矢量数据：在ArcGIS中输入基准DEM数据，运行水文分析工具进行填挖、流向、集水区、汇流分析、河网链接、斯特拉勒河网分级、河网矢量化、流域划分、流域矢量化一系列操作，再叠加地形和影像数据，检查输出的河网和流域矢量数据的质量，并人工校正，生成出流域边界分明、河网分布准确的流域和河网数据；[0048] 制作各流域边界、DEM和河流数据：将所有流域根据ID分割和根据ID赋名，将以ID为名的流域数据制作成掩膜，输入基准DEM数据分别生成各流域DEM数据，并裁剪各流域内的河流数据；[0049] 野外勘测数据采样：使用卷尺测量堆积扇上巨砾的长宽高信息；[0050] 生成人工解译灾害数据编目：筛选出低云量的遥感影像数据，在已有流域信息的基础上，结合野外验证和遥感影像上的形态、颜色、纹理信息，推断不同含水率及颗粒组成下的泥石流类别，并赋值各流域灾害类别代号；[0051] 人工解译流域沟道数量：根据遥感影像数据，目视解译随机样本流域内的沟道数目；[0052] 人工解译流域内堆积扇扩散角：根据遥感影像数据，目视解译随机样本流域内的堆积扇扩散角；[0053] 利用照相机或者无人机拍摄包括剖面、地貌、灾害影响范围和周边环境的照片。[0054] 本发明具有以下优点：一种基于堆积扇沉积特征识别的旱区泥石流活动评估方法，可自动化、重复进行大区域大范围的灾害分析，缩减了人工勘察和现场监测等低效方法的分析时间和成本，结果生成稳定、标准化、可迭代优化而适用于不同环境和场景。附图说明[0055] 图1本发明总体流程示意图；[0056] 图2数据组成及预处理示意图；[0057] 图3模型结构图；[0058] 图4沟道数量识别流程图；[0059] 图5堆积扇大小测算流程图；[0060] 图6最大冲击力度计算流程图；[0061] 图7灾害发生频率推算流程图；[0062] 图8指标模型构建流程示意图。具体实施方式[0063] 为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。因此，以下结合附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的保护范围，而是仅仅表示本申请的选定实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。下面结合附图对本发明做进一步的描述。[0064] 本发明具体涉及一种基于堆积扇沉积特征自动识别的旱区泥石流活动评估方法，其主要目的是基于堆积扇地形特征和沟道特征，遵循堆积扇分析和巨砾测量的传统原理，从依赖于人力目视识别的方式转化为高效的机器自动识别。其主要依托于python代码化人判断堆积扇特征的思维模式、识别过滤地形横断面中的谷信号（Troughsignal），实现自动识别堆积扇上沟道数量和堆积扇大小，通过差异化对比多时序数据推测灾害发生频率，结合分割任意物体模型（SegmentAnythingModel，简称SAM）开源分割技术识别堆积巨砾大小并还原灾害最大冲击力度，以指标形式综合展现成灾规模大小和发育频率。[0065] 如图1所示，其具体包括以下内容：[0066] 1、数据准备及预处理；[0067] 本发明使用的数据包括高程数据，高分辨率无人机或卫星影像，野外勘测数据以及人工遥感解译数据。这些数据在具体的参数建立、计算和优化中兼具不同的责任。由于建立的模型可复用，后续运用仅需当地高程数据、高分辨率无人机或卫星影像即可将本方法搬运到其他地区，或替换当地的地质和野外数据，使之因地制宜。[0068] 在本发明中使用的基准数字高程模型（DigitalElevationModel，简称DEM）数据为8组来源于（AlaskaSatelliteFacilityDistributedActiveArchiveCenter，简称ASFDAAC）的12.5米分辨率ALOSPALSARDEM数据集（2009年），以及用于对比分析获取灾害频率信息的美国国家冰雪数据中心（NationalSnowandIceDataCenter）提供的107幅光学影像生成的亚洲高山8米DEM数据集（2014‑2016年）作为对照数据。[0069] 如图2所示，数据准备及预处理步骤具体包括以下内容：[0070] （1）生成DEM数据：用ArcGIS的镶嵌工具分别镶嵌基准DEM数据和对照DEM数据（1组或多组，取决于数据量）。并向下重采样使镶嵌后的对照DEM与基准DEM分辨率一致，从而减少后续灾害频率计算中因数据分辨率差异所带来误差的影响。[0071] （2）划分流域、生成河网矢量数据：在ArcGIS中输入基准DEM数据，运行水文分析工具进行填挖、流向、集水区、汇流分析、河网链接、斯特拉勒（Strahler）河网分级、河网矢量化、流域划分、流域矢量化一系列操作。再叠加地形和影像数据，检查输出的河网和流域矢量数据的质量，加之适量的人工校正，生成出流域边界分明、河网分布准确的流域和河网数据。[0072] （3）制作各流域边界、DEM、河流数据：使用Python代码将所有流域根据ID分割，据ID赋名。将以ID为名的流域数据制作成掩模，输入基准DEM数据分别生成各流域DEM数据，利用GeoPandas库的clip模块剪裁各流域内的河流数据。[0073] （4）野外勘测数据采样：在实际野外调查中，使用卷尺测量堆积扇上巨砾的长、宽、高信息，并汇总整理。[0074] （5）生成人工解译灾害数据编目：筛选出低云量遥感影像数据，在已有流域信息的基础上，结合野外验证和遥感影像上的形态、颜色、纹理信息，推断不同含水率及颗粒组成下的泥石流类别，并赋值各流域灾害类别代号。其中各代号的意义为，0：无泥石流灾害（Nonetype）；1：高含砂水流（HCF）；2：水石流（DFLOOD）；3：黏性泥石流（VDF）。[0075] （6）人工解译流域沟道数量：根据遥感影像数据，目视解译随机样本流域内的沟道数目，并记录存储。[0076] （7）人工解译流域内堆积扇扩散角：根据遥感影像数据，目视解译随机样本流域内的堆积扇扩散角，并记录存储。[0077] （8）野外照片：野外考察过程中利用照相机或无人机拍摄包括剖面、地貌、灾害影响范围、周边环境等照片。[0078] 2、模型参数建立及处理；[0079] 如图3所示，本发明的模型参数由沟道数量、堆积扇大小、最大冲击力度和灾害发生频率组成。各参数的创新性体现在以下方面，分别是：基于单个流域内的沟道纵剖面特征定位沟口（堆积扇顶点）及堆积扇终点从而获取堆积扇横截中线及通过捕捉其谷信号推测沟道数量；通过堆积扇沟口和终点位置获得堆积扇半径及利用已识别堆积扇边界沟道计算堆积扇中线宽度及扩散角从而推算堆积扇面积；首次批量将图像识别方式应用于巨砾面积获取并通过野外考察所得面积与高关系推算巨砾高度和体积，结合所得堆积扇平均坡度计算灾害最大冲击力度；据不同时期的沟道识别数量获取时期内灾害发生频数及推算灾害发生频率。具体包括以下内容：[0080] A、沟道数量识别；[0081] 如图4所示，通过DEM数据处理、河流拓扑关系挖掘及坡度变化分析，自动识别沟口和堆积扇终点位置，并提取堆积扇中线剖面上的沟道数量。该方法的进步在于结合地形特征和统计阈值优化沟口识别，利用坡度变化和峰信号检测提取沟道数量，创新性地将地形信息自动化处理与人工校正结合，实现了高效、准确的地形特征分析和沟道识别。[0082] 具体步骤如下：[0083] （1）获取沟道纵剖面数据：将各流域DEM经纬坐标数据转化为平面图像像素坐标格式的数组。读取该流域内的河流矢量数据，根据ArcGIS生成河流数据自带的“grid_code（河网的分级）”筛选主河道数据，再根据“from_node（上节点）”“to_node（下节点）”属性根据数据拓扑连接规律从上游至下游交叉匹配数据形成连续河流数据。将上述河流数据坐标转化同DEM相同的像素坐标，获取河流数据点对应DEM数据中的高程值及根据点经纬坐标计算各点间实际地理距离，从而得到沟道纵剖面数据。[0084] （2）判断和优化沟口和堆积扇终点位置：根据各流域内沟道纵剖面数据，计算邻点间坡度、坡度变化值，以及沟道平均坡度变化、坡度变化标准差，用以预测最接近现实沟口的结果。本发明基于河流曲段裂点和地形变化的特点，使用CI阈值寻找由下游至上游地形突变的第一个正值特征点，沟口确定方程如下：[0085][0086] 其中， 、 和s分别代表坡度变化、平均坡度、标准差；t是临界t‑value，n是数据大小， 代表符合不等式要求的最低点的点坐标，即沟口坐标。[0087] 为优化数据结果，经多次实验显示，总体最优基准阈值为CI=98%，对比人工解译获得的沟口样本，此时预测的沟口与实际沟口产生的误差最小。为进一步得到更准确的阈值，根据不同流域内的灾害类型，发现该标准存在如下差异，即CI阈值在Nonetype类中取0.98，HCF类取0.997，DFLOOD类取0.965，VDF类取0.9925时，所得结果最优。根据最终所得流域的沟口位置保存沟口数据为shapefile文件。堆积扇终点则为流域内河流逆流方向上坡度变化急速转正的点，但若生成流域边界时流域底部边界恰好与堆积扇边界重合则可提取流域底部边界与河流的交点作为堆积扇终点以替代此步骤。[0088] （3）提取堆积扇中线剖面：根据所得沟口和堆积扇终点坐标，连接成线，计算其平均坡度和中点坐标，以中点为起点作正反向量垂线，即为流域中该堆积扇的中线剖面，提取剖面线上所有点数据分流域进行存储。[0089] （4）识别堆积扇横截中线上的沟道数量CN：使用python中的SciPy库的scipy.signal模块中的峰值检测功能识别堆积扇中线剖面中的局部最大值即主要峰信号和其对立信号即谷信号，保存预估谷数据和峰数据高程及剖面位置信息。利用高差变化筛除异常值，即移除信号中( )的数据，并设定沟道点高程最高不得超过流域堆积扇以上地形高程，得到调整沟道数量。基于所选实验区域单沟数据少，因此该方法适用于多沟改道型沟道的识别。结合实际解译的沟道数量数据，可获得调整沟道数量（AdjustedCN）和实际沟道数量（ActualCN）之间的数学关系如下，公式取得皮尔逊相关系数Pearson'sR为0.94，结果准确率较高。[0090][0091] 其中，Int表示向下取整函数，根据上述运算，优化了谷信号数据获取，并保存沟道数量数据为shapefile文件和优化信号数据为.csv格式文件。[0092] B、堆积扇面积大小计算；[0093] 如图5所示，通过计算堆积扇半径、扩散角及其校正值，测算堆积扇面积。其特点在于利用扇体计算公式和地形特征提取，提供了比人工测绘更快捷、稳定且重复性更高的堆积扇面积测算方法。[0094] 具体步骤如下：[0095] （1）计算堆积扇半径：根据已求得沟口和堆积扇终点位置及高程，根据三角形勾股定理求解堆积扇半径。[0096] （2）计算堆积扇扩散角预测值：根据优化沟道的双向极点，计算堆积扇中线长。而扇扩散角可通过以下方程计算。[0097][0098] 其中， 是单位为弧度的堆积扇扩散角预测值，arctan为反正切函数，是堆积扇中线长，R是堆积扇半径。[0099] （3）校正堆积扇扩散角。为保证上述计算的扩散角的准确性，检验人工解译样本得到的堆积扇扩散角数据与计算所得对应流域的扩散角数据之间是否存在数学关系。根据样本验证结果其之间存在以下关系。[0100][0101] 其中，α是单位为弧度的实际/校正堆积扇扩散角。所得皮尔逊相关系数Pearson'sR为0.89，因此该数学关系较准确。[0102] （4）推算堆积扇面积Afan。根据已知堆积扇半径（R）和校正堆积扇扩散角（α），可以计算得堆积扇面积（Afan）。此方式对比传统直接在遥感影像上人工描线测绘而言，快捷且输出稳定，不因不同操作人或不同测量标准和反复作业的不可重复性得到不同的结果。[0103][0104] C、最大冲击力度计算；[0105] 如图6所示，通过影像对干旱区堆积扇上巨砾进行颗粒分割，基于分隔特征物像素总数计算巨砾的顶面尺寸，通过野外调查所得巨砾尺寸关系公式推算巨砾高度和体积，最终结合地形和岩石参数估算泥石流最大冲击力度即侵蚀能力。其新颖之处在于利用高分辨率遥感数据批量识别、测量巨砾，通过挖掘巨砾尺寸相关变量间的关系填补二维影像中高度信息的空白，创新性地将图像分析与物理参数结合，实现对泥石流冲击力度的快速预测。[0106] 具体步骤如下：[0107] （1）分割堆积扇上的巨砾：在确保最佳分割性能的前提下，采用聚焦于小区域及扇体位置居中的高分辨率影像或无人机影像调用SAM进行颗粒检测分割，由于干旱区植被稀少的特点，图像中无除堆积体外的明显其他物象。[0108] （2）筛选合格图像数据：将输出结果依像素数量高低排序，计算特征物面积与图像总面积比筛除面积比大于等于20%、预测交并比和稳定系数均大于等于0.9的特征物，依流域保存合格数据及最大值到Excel文件。[0109] （3）计算所用图像单个像素实际尺寸。根据相似三角形定理，依据图上距离和拍照设备参数可计算图上单个像素代表的实际尺寸，计算公式如下：[0110][0111][0112][0113][0114] 其中，Paw和Pah是像素的实际宽和高，即遥感影像的空间分辨率；Ppw和Pph是图像上像素宽和高；d代表拍摄距离；而w和h分别为镜头焦平面的宽和高；f是焦距；W和H则是图像分辨率宽和高，单位为像素；ϵflight和ϵDEM是拍摄设备所处海拔和图片中心高程； 是单个像素的实际面积。[0115] （4）计算巨砾实际尺寸：根据图上像素的实际面积以及特征物的像素数量，可得到特征物即巨砾的实际顶面面积（长宽所得）。[0116][0117] 上式中N是一个特征物在图像上的像素数量，而 是地面物体的实际顶面面积。[0118] （5）建立野外巨砾面积与高关系：根据野外考察所得巨砾长、宽、高数据，发现巨砾长宽所得面积与高存在关系如下：[0119][0120] 其中， 是巨砾的高， 是巨砾的长‑宽所得面积，所得皮尔逊相关系数Pearson'sR为0.92，表明该方程关系较准确。该方法在二维影像中巨砾高度未知的情况下，可据公式较准确的估计巨砾高度，填补各流域内巨砾高度信息的空白。[0121] （6）推算巨砾体积：根据以上步骤可推出巨砾的高，用面积和高的乘积可算得巨砾的体积。[0122][0123] 其中， 是巨砾体积。该方式结合了可测巨砾面积和估算巨砾高度数据，用已知推测未知。[0124] （7）代入参数计算最大冲击力度τ*Max：使用最大冲击力度公式，进行计算。[0125] 代入公式需要的各参数值，读取根据上述所得巨砾体积和已知巨砾面积推测巨砾高度、以及岩石密度 、岩石内摩擦角 、堆积扇平均坡度θ、重力加速度 ，可推算出不同流域的泥石流最大冲击力度（τ*Max:Pa），自动连接所有数据及流域ID保存到Excel文件中。[0126] D、灾害发生频率获取；[0127] 如图7所示，灾害频率的计算是基于两次或多次不同时期沟道数量识别的结果差和时间差。相较于过去使用监测或灾史记录灾害发生次数而言，该方法仅需不同时期同一区域相同分辨率的DEM，具有一定的普适性。灾害发生频率计算结果的准确性依赖于沟道数量识别结果的准确性。[0128] 具体步骤如下：[0129] （1）识别对照DEM中出现的沟道数量：将对照DEM数据重复数据准备及预处理步骤中（3）步骤及模型参数建立及处理步骤中A部分的（3）和（4）步骤，得到对照时期的沟道数量数据。[0130] （2）计算沟道不同时期沟道数量差。将对照DEM中的各流域的沟道数量与基准DEM中各流域的沟道数量进行减法运算，公式如下所示：[0131][0132] 其中， 是沟道数量差， 是对照DEM对应的沟道数量， 是基准DEM对应的沟道数量。[0133] （3）计算不同DEM的时间差：时间差即基准与对照DEM年代属性的差，公式如下所示：[0134][0135] 其中，ΔT是以年为单位的时间差，TA是对照DEM对应的年代，TB是基准DEM对应的年代。[0136] （4）推算灾害发生频率f：灾害发生频率是沟道数量差和时间差的比，公式如下所示：[0137][0138] 若所得f即灾害发生频率（单位：次/年）的值偏低，可相应的乘以10或50或100，将其换算为次/十年或五十年或百年。[0139] 3、模型构建及指标化；[0140] 如图8所示，将以上所得四个参数的结果统一标准化处理后作为主成分分析法（Principalcomponentanalysis，简称PCA）的输入变量，求取对模型整体贡献量较大的主成分及数量，根据数据本身的特点获取各主成分在模型中的权重，构建内容全面的、合理取值范围的模型结果，实现指标的标准化并适用于不同区域或条件下的结果的可比性。相比据经验主观地率定参数权重而言，该方法的优势在于依据数据本身的差异性差异化分析各因素对欲构建模型的影响，科学客观且较全面地评价区域泥石流规模和活动性特征。[0141] 具体步骤如下：[0142] （1）标准化各参数变量：将各参数变量所得结果合并到同一个Excel页中，批量归一化，在确保所有变量享有平等的可贡献程度的前提下，保留数据的分布特点。[0143] （2）运行PCA检测：调用Scikit‑learn库sklearn.decomposition模块里的PCA检测功能，将高维数据投影到低维空间中，实现PCA降维运算，并对标准化的数据进行拟合和转换，返回主成分得分。[0144] （3）获取每个主成分的解释方差比率并筛选模型构建所需主成分：计算每个主成分对数据总方差的独立贡献即解释方差比（Explainedvarianceratio，简称EVR）。为保证选定的主成分保留足够的数据信息，累积解释方差比（Cumulativeexplainedvarianceratio，简称CEVR）即数个从大贡献率到小贡献率排列的EVR之和需超过总方差的80%‑90%阈值，以此选出构建模型所需主成分。[0145] （4）生成主成分序列并计算载荷矩阵：生成从1到主成分个数的主成分序列，利用主成分特征向量表示原始变量即上述四个模型参数变量在主成分空间中的贡献，获取PCA分析的载荷矩阵。该载荷矩阵表示各参数变量在主成分上的系数，代表着主成分的构成形式，公式如下所示：[0146][0147] 其中， 是第j个主成分， 是第k个变量的值， 是第j个主成分的第k个变量的解释方差，M是变量总数。[0148] （5）计算各参数变量间的协方差矩阵：根据协方差函数计算各变量间的协方差，并从协方差矩阵中提取主成分间的协方差值。[0149] （6）基于归一化协方差值分配权重并计算各成分权重。计算上述提取主成分协方差值的总和，按照各协方差值占总和的比例确定各成分的归一化权重。验证所得权重总和是否为1，以确保权重分配准确性。最终形成指征泥石流规模和活动性的指标（Indexscale‑activity）模型如下：[0150][0151] 其中， 是第j个主成分 的权重，m是数据集中样本总量。[0152] 在此发明所用数据集中，指标模型结果如下。[0153][0154][0155][0156][0157] （7）运行模型进行泥石流活动评估：输入泥石流相关数据，即获得CN、 、 和f所需的所有数据，运行Indexscale‑activity模型，获得泥石流活动Indexscale‑activity结果，其中，值越大，泥石流暴发频率越高，规模较大；[0158] （8）按需扩缩指标数值范围：根据数据结果和数值特征，为方便分析数据展示结果，可适当重分配数值范围为0‑10或0‑100，或其他。为利于样本数据结果分析，将本发明中所得指标结果按0‑10重赋值。[0159] 以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和完善，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

专利地区：四川

专利申请日期：2024-06-24

专利公开日期：2024-09-03

专利公告号：CN118379649B