具体实施方式

下面结合附图对本发明实施例进行详细描述。应当明确，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。本领域技术人员应当知晓，下述具体实施例或具体实施方式，是本发明为进一步解释具体的发明内容而列举的一系列优化的设置方式，而该些设置方式之间均是可以相互结合或者相互关联使用的，同时，下述的具体实施例或实施方式仅作为最优化的设置方式，而不作为限定本发明的保护范围的理解。

在一个具体的实施例中，结合具体的实例，对本发明的方案的推导及运用方式进行阐述。

首先，以某市的实际地质数据为例，来分析河道距离、沟口距离与灾害点之间的关联关系。

1、灾害点与河道距离、沟口距离分析

以某市为例，对河流做1000m、2000m、3000m、4000m的多级缓冲区，统计不同缓冲区内的灾害点个数，统计结果如图1所示。在距离河道1000m以内的历史泥石流灾害点共104个，占比81.25％，在距离河道1000-2000m以内的历史泥石流灾害点共17个，占比13.28％，随着距离的增大，泥石流灾害点逐渐减少。表明大部分泥石流灾害都发生在河道口或者河道周围。

泥石流从沟口冲出，由于地形突然开阔，沟道的相对断面面积增大，基于瞬时流量相对恒定不变的前提，泥石流运动的直接表现为速度降低，并发生堆积，直到达到最大堆积长度和宽度，一般认为，在泥石流达到堆积扇的边缘时，其运动速度接近于零。泥石流在发生堆积时，其运动方向可以分解为纵向和横向两个方向，同理其流速也可以分解为纵向(垂直方向)和横向(水平方向)两个分速度。纵向分速度以沟口最大，堆积扇最前端最小，对应着泥石流最大堆积长度；沟口距离因子可以一定程度上展现泥石流的动量和破坏力大小。以流域出口为沟口，计算栅格像元距沟口的直线距离作为沟口距离。为了探究泥石流灾害与泥石流沟沟口距离之间的关系，以500m为缓冲距离，对沟口距离进行反距离归一化计算：

其中，R为缓冲距离阈值，d为沟口距离。

在一个优选的实施方式中，可以将该阈值R设置为例如500，则此时，L值计算变为：

将历史128个泥石流灾害点叠加至图层进行统计。距离沟口距离越近泥石流灾害发生频率越高，在0.8-0.6范围内灾害点共18个，占比14.29％，在0.8-1范围内灾害点共79个，占比63.7％。随着沟口距离量化值越大泥石流灾害发生频率也逐渐增大，如图2所示。表明泥石流灾害大部分都发生在沟谷沟口及附近。

通过将泥石流灾害点河道距离以及沟口距离进行分析，发现泥石流灾害大多发生在流域沟口位置以及靠近河道部分，而径流汇流模拟是基于沟谷、河道上的模拟。因此。以径流汇流模型为基础的泥石流灾害预警，能提高灾害预警的空间针对性。

2、泥石流灾害模型的构建

泥石流灾害是地质、地形、地貌、土壤、植被、降水与冰川融水等因素共同作用的结果。相等规模的降水或融雪，因为地形、地质土壤等条件不同而产生不同程度的泥石流灾害。由于地质、地形、地貌等因素在一定的时间内保持相对的稳定，因此，本方案中将其作为区域的背景因子，降水和融雪作为泥石流灾害的激发因子，构建泥石流灾害预警模型。

首先，依据DEM进行小流域划分。将背景因子纳入到不同的流域中，分析流域的背景因子对泥石流灾害的影响，构建基于流域单元的泥石流灾害的危险性概率模型。危险性概率结果是对泥石流活动程度的一种客观展现，危险性程度越高，表明相同的降水、融雪导致泥石流灾害的可能性越大。然后，引入降雨融雪—径流汇流模型模拟径流水深。将模拟水深作为泥石流灾害预警因子，建立基于径流模拟水深、危险性概率背景值因子、沟口距离因子的泥石流灾害预警模型。泥石流灾害预警模型构建思路如图3所示。

3、泥石流灾害危险性概率计算

3.1、流域划分

目前，泥石流灾害的危险性区划大多建立在单个栅格网格单元的基础上。单个栅格单元难以量化流域整体特征，如流域面积、纵比降等，而这些因子是构建泥石流模型的重要因子。因此，本发明利用区域DEM，基于水文D8法将DEM划分为大量的小流域，作为泥石流灾害的预警基本单元。

水文D8法是栅格水文分析的基本方法。栅格水流仅能从中心像元流动至其中一个邻近像元。由于每一个栅格像元有8个邻近像元，故称为水文D8法。水文D8法实际上是一种单流向水文分析方法，ArcGIS的水文分析模块即采用水文D8法。

3.2、泥石流灾害诱发因子的选取与逻辑回归分析

诱发泥石流灾害的因素较多，形成机理较为复杂。依据泥石流发生机理，在一个优选的实施方式中，本方案选择裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量因子作为泥石流灾害诱发因子，基于信息量结合逻辑回归的方法确定泥石流灾害危险性概率。当然，对于泥石流灾害诱发因子，也可以选用其他因素，或者使用本方案的上述多个因素的任意组合方式，具体可以依据灾害区域差异、灾害判断侧重点差异等进行设置，此处不再一一赘述。

各单因子信息量I(K,x_n)的计算方法如下：

式中：S为某市地区评价单元总数；N为某市地区发生泥石流灾害所分布的总单元数；S_i为某市内含有评价因子x_n的单元数；N_i为分布在评价因子x_n内特定类别的泥石流灾害单元数；K表示是否为泥石流灾害点。

上述单因子信息量I通过以往历史数据计算得到。计算出的各个单因子信息量，作为逻辑回归模型中的影响地址灾害发生的因子的历史数据，从而通过逻辑回归方式，结合上述历史数据计算获得的单因子信息量I，计算逻辑回归模型中的各影响因子对应的逻辑回归系数B₁…B_n。

逻辑回归模型建模方法如下：

上式中，P为泥石流灾害危险性概率，取值在[0，1]之间，(1-P)为泥石流灾害不发生的概率，X₁…X_n是影响地质灾害发生的因子，B₁…B_n是各影响因子对应的逻辑回归系数，A为回归方程常数项(常量)。

以某市为例，将剔除误差数据后的701个隐患点数据记录作为建模数据，131个有效灾害点数据作为检验数据。使用ArcGIS软件栅格数据提取工具提取地质隐患点和未发生地质灾害点致灾因子信息量值，以泥石流灾害是否发生为因变量，裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量信息量值为自变量在SPSS软件中建立逻辑回归模型，逻辑回归计算得到表1。

表1逻辑回归分析结果

各因子均通过显著水平为0.05的Wald检验，且因子间相关系数很小，相互独立；多次随机生成未发生泥石流灾害点并循环建模，模型结果稳定；最终确定泥石流灾害危险性概率如下：

上式中，H为泥石流灾害危险性概率值，X₁～X₁₁分别为裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量的信息量值。

将泥石流灾害概率分为五级，利用131个某市历史泥石流灾害点来对泥石流灾害概率进行精度验证，各灾害点所处泥石流灾害危险性概率所占比例如表2所示。

表2泥石流灾概率级别占比

在五级内的历史灾害点共110个，占总历史灾害点83.79％，四级内的灾害点共2个，占0.02％，在二级及一级内的灾害点共19个，占比11.48％。从表2中可见，验证数据与某市泥石流灾害危险性概率结果匹配程度较高。这说明本发明的泥石流灾害的危险性概率计算的方法合理性，为泥石流灾害预警模型的建立奠定了基础。

4、泥石流灾害预警模型

降雨是泥石流灾害发生的诱发条件，因此，本申请的方案中，作为一个优选的方式，将径流汇流、泥石流灾害危险性概率以及沟口距离三个因子作为泥石流预警主要因子对泥石流灾害进行预警。

4.1、泥石流灾害预警方法

泥石流灾害主要在河流沟谷附近。因此，以径流汇流模型模拟的水深替代传统的面雨量进行泥石流灾害预警，以提高预警的空间针对性。本发明提出一种通过建立汇流水深与泥石流灾害发生背景之间的关系来建立泥石流灾害预警模型的方法。设径流汇流水深分级为G，泥石流灾害危险性概率背景值为H，沟口距离因子L。利用二元逻辑回归的方法计算各个因素对泥石流灾害的计算权重，如下公式：

Y＝(aG+bH)×L

式中：Y为泥石流灾害预警结果，G径流汇流水深分级，H为泥石流灾害危险性概率背景值，a、b分别为径流汇流水深及泥石流灾害危险性概率背景值权重。

4.2、径流汇流模型

4.2.1、径流汇流模拟

在本发明中的一个优选的实施方式中，对纳维-斯托克斯方程进行栅格模式分解，按平流项、压力项和外力项分别进行时间片上每一个栅格的水深和水速计算。通过时间片的迭代，模拟流域的径流汇流过程。由于冰川消融产生的径流也是泥石流灾害发生的重要因素，将模型中的水量平衡公式设置为：

式中，分别为时间片上中心栅格c初始水深和输出水深； 分别为中心栅格时间片上的降水量、冰川消融量、下渗量、蒸发量；和为模型中平流项。

栅格降水量由气象站小时降水插值得到，冰川消融量由日消融量换算至小时消融量得到。与为平流项模型，其计算方法如下：

式中，为水流在邻近栅格b水向中心栅格c流动时所产生的速度分量，为水流在中心栅格c水向邻近栅格b流动时所产生的的速度分量，ΔL为两个栅格中心之间的距离。

基于与设周围邻域的栅格水流速度矢量为邻域栅格流体流动时到中心栅格的水量S_b→c为：

式中，(Δx，Δy)为单位像素栅格之间所产生的坐标偏差；v_max为水流流动时在栅格相等时间片上所产生的最大速度s。

计算径流中心栅格剩余水流量S_c为：

计算之后中心栅格水深为：

式中为中心栅格邻域在t_n时间片上的水深，为t_n时间片上邻域栅格b在流体运动时进入到中心栅格c中的水量比例，为t_n时间片上中心栅格剩余的水量比例。

水体在模拟交互运动的过程中满足动量守恒定律。设和表示中心栅格的初始水速和终止水速。则终止水速为：

水体在地表流动时，由于土壤及气候环境的影响，计算水体流动的下渗和蒸发量，在模型中认为同一位置上蒸发率与模拟时间段内保持不变。

式中，e为率定参数。

水体在地表运动时，下渗程度受多种因素影响，利用利普下渗曲线进行模拟计算：

f_t＝f(t)＝0.5×S·t^-0.5+A

式中S与A为与土壤性质有关的参数。

地表水在流动过程中受重力、压力的影响产生相应的加速度，用ΔH_b来代表单元之间水体流动的高程差，来计算相邻栅格对中心栅格所产生的速度矢量

式中j和i为水流流动模拟过程中邻域栅格与中心栅格的列序号和行序号差；α为与加速度、水的密度、摩擦力等相关的率定正数常量。

外力项和压力项统一利用水动力方法进行模拟。水体在地表流动的过程中时会受地表摩擦力影响而使水流速度发生纵向的变化，当实际水深低于给定的阈值时，设定一个衰减系数以减小模拟水速，衰减系数ε近似定义为：

式中d_min和d_max为水体流动所受摩擦力的下界深度和上界深度，单位为m，σ为流体的摩擦比例系数。

时间片上各栅格单元仅考虑与8邻域栅格水量交换，用于时间片内水流最大距离不超过栅格分辨率。为流域设定最大水速限制v_max，模拟水速大于v_max强制设置为v_max。在这种条件下，相邻两时间片的时间差Δt为：

Δt＝C/v_max

式中，C为栅格分辨率。

根据以上公式，每小时需要迭代的时间片数N为：

N＝3600/Δt

4.2.2、冰川融雪模型

冰川消融型泥石流也是泥石流灾害的主要类型之一，冰川消融所带来的水分也为周围流域沟谷形成泥石流灾害提供了一定的基础。本发明将冰川消融的部分加入径流模型模拟中，作为径流模拟初始水源部分。在现有技术中，目前大都采用能量平衡的方法对冰雪消融进行模拟，能量平衡的计算复杂并且需要大量资料，在资料稀缺的地方，能量平衡算法在当地难以推广，严重影响冰雪消融因素对泥石流影响的精确计算。在本发明的一个优选的实施方式中，采用如下方式作为冰川融雪模型并且加入至径流汇流模型中：

冰雪消融模型基本计算公式如下：

式中，M为日消融量(mm)，T为日平均气温(℃)，T₀为冰川消融温度阈值(℃)，MF为度日因子(mm·℃－1·d－1)。

将气温数据插值后代入上式计算得到冰川所在栅格像元的日消融量。栅格单元的降水量加上冰川消融量作为模型的水量输入。

4.2.3、径流汇流的水深分级

以某市为例，选择2015年8月19日00时起模拟18天，将气象站和气温数据代入径流汇流模型，模拟某地区的径流汇流过程，并以某江河段下游的实际观测水深与模拟水深进行对比。图4为模拟水深与实际水深对比曲线。从图4中可见，模拟水深曲线与实际水深曲线运动趋势一致。将实际水深数据与模拟水深数据进行pearson相关性分析，其相关性为0.838，为强正相关，表明本发明所采用的能径流汇流模型模拟出水深的变化趋势。

虽然模型模拟的水深变化趋势接近，但模拟水深的变化值与实际水深变化值不在相同量级上。这是因为，本发明的径流汇流模拟在干涸河道条件上进行的，未考虑河道初始水位；同时，模型输出的水深是以栅格像元为单位的，相当于将像元内的水平铺于矩形栅格单元上。因此，不能将模拟水深的绝对值直接代入模型中进行泥石流预警。在一个更为优选的实施方式中，本发明采用水深分级方法进行泥石流灾害预警模型的构建，表3为径流水深分级表。

表3径流水深分级

4.3、基于逻辑回归的泥石流灾害预警模型构建

本文基于区域历年泥石流灾害点、径流汇流模拟水深分级、泥石流灾害概率值来建立二元逻辑回归模型。以某市为例，选取2000年-2015年的历年泥石流灾害点数据参与逻辑回归运算，以泥石流灾害是否发生作为因变量，径流汇流模拟水深、泥石流灾害危险性概率背景值作为自变量。在SPSS软件中建立二元逻辑回归模型得到最终权重系数如下表：

表4逻辑回归方程表

其中，B为逻辑回归系数；S.E.为标准误差；Wald为卡方值；df为自由度；Sig.为显著性。

对泥石流灾害指标进行逻辑回归后，其所有因子显著性都小于0.05，基本可判定该逻辑回归模型具有显著性，逻辑回归模型通过检验。将B值进行量化，量化后作为预警因子的权重，根据下式,逻辑回归构建某市的泥石流灾害预警模型的表达式如下：

Y＝(0.52G+0.48H)×L

4.4、泥石流灾害预警精度评价

以某市为例，对本发明方法进行泥石流灾害预警精度评价。将2010-2019年的泥石流灾害点以及相同数量的随机点投入模型中，经计算获取Y值的范围，未发生灾害的随机点的T值分布与0到0.45之间，极少数大于0.45；而灾害点大部分分布于0.55到1之间。通过分析确定泥石流灾害预警Y值超过0.55发布黄色预警信号，超过0.65发布橙色预警信号，超过0.75发布红色预警信号。将泥石流灾害风险等级分为以下五个等级：

表5泥石流灾害预警等级划分

选取部分样本内历史灾害点以及2011年-2020年历年泥石流灾害点数据作为结果验证数据得到验证结果如下表6。

表6泥石流灾害预警验证结果

根据预报结果，选取案例历史灾害点数据共41个灾害点，预报等级在黄色以上灾害点为27个，占比65.85％。预报等级在橙色以上灾害点为26个，预报精度63.41％。预报等级在红色以上灾害点为26个，预报精度63.41％。2016年以前未预测到灾害点大部分原因是气象站较少，降雨数据缺失，导致部分地区径流汇流模拟结果欠佳即预警效果降低。自2016年以后由于开始逐渐新增气象站点，降水数据更为全面，预警模型精度也在一定程度上得到了提高。

图5为实际应用本发明提供的灾害预警方法后，所预测的该市2019年7月7日泥石流灾害预警图。与实际当天数据对比，当天发生泥石流灾害2处，不同颜色级别的预警共4个(其中包含2个不同预警级别的各2个)。这说明本文的泥石流预警模型具有较高的精度，能应用于实际的泥石流灾害预警预报工作中。

此外，在又一个具体的实施例中，本发明的方案还可以通过一种基于径流汇流模拟的泥石流灾害预警系统来实现，该系统包括：

泥石流灾害危险性计算模块，用于确定预警区域内，泥石流灾害危险性的概率背景值H；

水深计算模块，用于基于栅格化处理后的纳维-斯托克斯方程，进行径流汇流水深W的计算；

沟口距离因子计算模块，用于基于泥石流缓冲距离，获取泥石流灾害危险性相关的沟口距离因子L；

灾害预警结果计算模块，用于基于泥石流灾害危险性概率背景值H、径流汇流水深分级G、沟口距离因子L，计算泥石流灾害预警结果Y；

预警模块，用于基于泥石流灾害预警结果Y，进行泥石流灾害预警；

数据输入模块，用于获取预警区域泥石流灾害点数据。

优选的，所述泥石流灾害预警结果Y的计算方式为：

Y＝(aG+bH)×L

式中：Y为泥石流灾害预警结果，G为径流汇流水深分级，H为泥石流灾害危险性概率背景值，a、b分别为径流汇流水深及泥石流灾害危险性概率背景值的权重。

优选的，泥石流灾害危险性概率背景值H的获取方式为：

建立泥石流灾害危险性概率P与各单独诱发因子的逻辑回归关系：

其中，P为泥石流灾害危险性概率，X₁…X_n是影响地质灾害发生的单独诱发因子，B₁…B_n是各单独诱发因子对应的逻辑回归系数，A为回归方程常数项；

基于历史数据，计算各单独诱发因子对应的逻辑回归系数，此时的泥石流灾害危险性概率P作为泥石流灾害危险性概率背景值H。

更为优选的，该系统还可以实现如上所述的基于径流汇流模拟的泥石流灾害预警方法。

并且，本领域技术人员应当明了，本发明的方案还可以通过计算机可读介质，或者一种包含存储装置、处理器装置的设备，来实现以上系统，并且可以使得该系统在运行时，可以实现本发明说明书中所记载的基于径流汇流模拟的泥石流灾害预警方法，而上述的实现方式及其在此基础上的常规模块调整及更改，均应当视为落入本发明的保护范围之内。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random AccessMemory，RAM)等。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。