CN110009002A - 基于多维降雨特征空间最优决策的山洪快速预警预报方法 - Google Patents

Abstract

本发明公开了一种基于多维降雨特征空间最优决策的山洪快速预警预报方法，首先基于三层蒸散发原理计算逐次历史雨洪过程的初始土壤水饱和度S；其次统计逐次历史雨洪过程的降雨量P；接着，结合流域DEM数据，量化逐次历史雨洪过程的降雨空间分布特征K；基于DEM数据计算流域的地貌单位线；并以P，K和S为相互正交的坐标轴构建多维降雨特征空间坐标系，构建降雨数据点集，依据洪水过程是否超警，将对应的降雨数据点集划分为超警子集A与未超警子集B两类，结合SVM原理求解两类子集之间的最优分类决策方程。若需要进行预警则以该降雨为地貌单位线的输入，计算流域出口处的径流过程，进行实时洪水预报。

Description

基于多维降雨特征空间最优决策的山洪快速预警预报方法

技术领域

本发明涉及水文技术领域，尤其涉及一种基于多维降雨特征空间最优决策的山洪快速预警预报方法。

背景技术

我国汇流面积200至3000km2的中小河流近9000条。近年来，受气候变化影响，由局地强降水造成的中小河流突发性洪水频繁发生，已成为造成人员伤亡的主要灾种。据统计，我国中小河流洪涝灾害和山洪地质灾害损失约占全国洪涝灾害经济损失的70％～80％，死亡人数占2/3左右，对人民群众生命财产安全构成了严重威胁，因此对中小流域突发性洪水进行准确可靠地预警预报成为亟待解决的重要问题。

由于大部分中小河流源短流急，洪水具有历时短、上涨快、难预报的特点，基于水文模型构建的预报方案常常难以对快速涨落的山洪做出及时有效的反应，因此直接基于降雨数据进行山洪快速预警的方法为中小流域洪灾防治提供了一个新的思路。2006年11月我国国务院批复的《全国山洪灾害防治规划》中指出：当一个流域或区域某一时段内的降雨量达到或超过某一量级时，该流域或区域将发生山溪洪水、泥石流、滑坡等山洪灾害，把这时的降雨量，称为该流域或区域的临界雨量。由于可以依据实时监测的降雨来决策是否发布相应等级的山洪预警，因此该方法不需要过度依赖对于未来降雨状况的准确估计，同时也不需要进行预报方案的构建，这可以为应对灾害争取更多的时间，所以直接基于降雨数据来进行山洪快速预警的方法被广泛应用。目前基于降雨进行山洪预警的临界指标主要局限于降雨量，但对实际的山洪灾害预警而言，每场降雨的时程分配均不相同，这可能使实际的降雨过程与拟定山洪灾害预警指标所依据的雨型不一致，从而导致预警不准确。除了降雨时间分布之外，降雨空间分布的不均匀性也可能对流域中的坡面产流过程产生较大影响，由降雨空间分布的不均匀性带来的产流峰值的变化差距可以超过40％。即同一个流域，降雨过程不同，但成灾流量及其它条件相同，所需要的临界降雨量是不同的。基于遥感、地理信息以及数字流域等技术的发展，采用数值矩阵描述地表高程变化的栅格数字高程模型(DEM，Digital Elevation Model)逐步成熟，并得到了广泛的应用。尤其是在地形复杂的山区性中小流域中，DEM数据因其能够较为准确地考虑流域内地形变化而具有重要的应用价值。如何利用DEM数据来对临界降雨量计算方法进行完善，从而在一定程度上考虑地理特征空间分异性对于临界降雨量计算的影响，也是基于临界降雨量进行山洪预警预报走向分布化，精细化过程中的重点和难点之一。

为了进一步促进水文模型中产流计算的发展，需要更深入理解地理特征空间分异性对于产流计算的影响，研究蓄满与超渗空间组合模式的分布式产流计算方法。

要将降雨的时空分布特征作为预警指标，首先面临的挑战就是将降雨的时空分布特征进行量化。在目前的实际应用中，主要采用“分类”的方法来考虑不同降雨过程之间时程分布特征的区别，即考虑雨峰形状各异的不同雨型。在空间分布方面依据降雨中心的位置大致划分为发生在流域的上游，中游以及下游的降雨。在研究降雨的时空分布特征对洪水产生过程带来的影响时，分类的方法可以基本满足要求。但是将降雨的时空分布特征作为山洪预警的临界指标时，则需要将其量化为一个确切的数值，即“定量”。基于多维降雨特征进行山洪预警，关键是要确定发布预警时的降雨特征临界状态。在多维降雨特征构成的空间中，依据洪水过程是否需要预警对具有不同特征的降雨过程进行分类，找到最优分类决策面，即是否需要发布预警的临界面。而在临界面中，土壤饱和度-降雨量-降雨空间分布-降雨时程分布之间相互影响，在改变其中任意一个因子后，其余的因子的值也会随之发生变化，进而影响到多维空间中最优分类决策方程的求解。这对一维数值空间中进行临界判断的传统决策方法提出了挑战，因此如何定量的描述降雨特征之间的相关关系，并准确地求解最优分类决策方程是值得进一步探讨的问题。因此土壤饱和度，降雨量，降雨空间分布，降雨时程分布是影响临界雨量的主要因素，若是忽略这些因素对于临界雨量的影响，则会导致临界雨量计算的不准确，从而为山区性中小流域的山洪防治留下隐患。

因而对于多维降雨特征的忽视不利于基于临界雨量的中小流域的山洪防治的发展。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及到的缺陷，提供一种基于多维降雨特征空间最优决策的山洪快速预警预报方法。

本发明为解决上述技术问题采用以下技术方案：

基于多维降雨特征空间最优决策的山洪快速预警预报方法，包括以下步骤：

步骤1，基于三层蒸散发原理计算逐次历史雨洪过程的初始土壤水饱和度S；

步骤2，统计逐次历史雨洪过程的降雨量H；

步骤3，结合流域的数字高程模型数据、即DEM数据，量化逐次历史雨洪过程的降雨空间分布特征K；

步骤4，量化逐次历史雨洪过程中的降雨时程分布特征G；

步骤5，基于DEM数据计算流域的地貌单位线；

步骤6，以H，K，S和G为基础构建多维降雨特征空间坐标系，并构建降雨数据点集，依据洪水过程是否超警，将对应的降雨数据点集划分为超警子集A与未超警子集B两类，结合向量机原理即SVM原理求解两类子集之间的最优分类决策方程；

步骤7，计算实时降雨过程的降雨量、降雨空间分布特征、土壤水饱和度和降雨时程分布特征，并基于最优分类决策方程判断该实时降雨过程是否需要发布预警，若需要进行预警则以该降雨为地貌单位线的输入，计算流域出口处的径流过程，进行实时洪水预报。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤1的具体步骤如下：

步骤1.1，整理得到逐小时不间断历史降雨径流资料，并依据径流过程的涨落均匀摘录出若干次雨洪过程；

步骤1.2，以逐小时不间断降雨资料为输入，采用三层蒸散发原理计算得到逐小时不间断的土壤水饱和度变化过程，并从中提取出摘录出的若干次雨洪过程中开始时的土壤水饱和度，记为S。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤2中根据以下公式计算每次历史雨洪过程的降雨量H：

式中，T为雨洪过程的总时段数，t为时段的编码，从1到T；P_t为t时段内降雨。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤3的具体步骤如下：

步骤3.1，计算流域中每一个栅格单元中的汇流流速；

利用流域DEM数据将流域划分为若干个正交的栅格单元，计算出流域中每一个栅格单元的坡度、流向和汇流累计值，并提取流域中的河道栅格和坡地栅格；其中，根据以下公式计算栅格单元中的汇流流速v；

式中，D_p为坡面栅格中的坡度；θ为预设的第一经验系数阈值；A为河道栅格单元中的累计汇流面积，D_r为河道栅格单元中的坡度；a、ε、β分别为预设的第二、第三、第四经验系数阈值；

步骤3.2，根据每一个栅格单元中的汇流流速，得到各个栅格单元中的水流运动到流域出口所需要的时间，得到流域中的汇流时间栅格Raster_Htime；

步骤3.2.1，根据以下公式计算水流经过每一个栅格单元中所需要的时间Ctime：

式中，Csize为栅格单元的边长，即流域DEM数据的分辨率；dir为栅格单元的流向，即该栅格单元中的水流流动的方向，dir等于1、2、3、4、5、6、7、8时栅格单元的流向分别为东，东南，南，西南，西，西北，北，东北；

步骤3.2.2，根据以下公式计算编号为i的栅格单元中的水流运动到流域出口所需要的时间Htime_i：

式中，N为各个栅格单元中的水流运动到流域出口所经过的栅格单元的总数；i为各个栅格单元中的水流运动到流域出口所经过的栅格单元的编号，从1到N；Ctime_i为水流经过编号为i的栅格单元所需要的时间；

步骤3.2.3，将流域中所有栅格单元的Htime_i组合得到汇流时间栅格Raster_Htime；

步骤3.3，对逐次历史雨洪过程中站点观测的雨量数据进行空间插值，得到降雨空间分布栅格，并结合汇流时间栅格得到降雨空间分布特征K：

式中，F为流域中栅格单元的总数，c为流域中栅格单元的编码，从1到F；SP_c为编号为c的栅格单元中的降雨量，由逐次历史雨洪过程中站点观测的雨量数据进行空间插值得到，MaxP为流域中所有的栅格单元中最大的降雨量；Htime_c为编号为c的栅格单元中的汇流时间，MaxT为流域中所有的栅格单元中最大的汇流时间。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤4的具体步骤如下：

步骤4.1，将逐次历史雨洪过程中的降雨数据由小到大排序，得到NSeq_P；

NSeq_P＝{P₀,P₁…P_t…P_T}

式中，P_t表示t时段内的降雨量，t从0变化到T；且P₀<P₁<…<P_t<…P_T，即随着下标t的增大，P_t的值也在不断增大；

步骤4.2，基于NSeq_P计算得到累积时段雨量序列ASeq_P；

ASeq_P＝{AP₀,AP₁…AP_t…AP_T}

式中，AP_t为时段为t时的累积时段雨量，且AP₀<AP₁<…<AP_t<…AP_T，即随着下标t的增大，AP_t的值也在不断增大，最大的AP_T的值为H；

步骤4.3，计算Seq_P中的累积时段雨量百分比，得到累积时段雨量百分比序列ANSeq_P；

NSeq_P＝{APF₀,APF₁…APF_t…APF_T}

式中，t从0变化到T；APF₀<APF₁<…<APF_t<…APF_T，即随着下标t的增大，APF_t的值也在不断增大，最大的APF_T的值为1；APF_t的计算公式如下：

步骤4.4，基于ASeq_P和ANSeq_P计算G；

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤5的具体步骤如下：

步骤5.1，基于步骤3.2中计算得到的汇流时间栅格统计得到流域中的栅格数—时间关系；

在Raster_Htime的基础上，以Δtime为间隔进行分组统计，得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到流域中的栅格数—时间关系；

式中，Tcount_z表示汇流时间大于等于Z×Δtime、小于(Z+1)×Δtime的栅格单元的数量；Z为汇流时间分组的编号，从0到SZ，其中SZ＝MaxT/Δtime；Δtime为预设的时间间隔阈值；τ_c为编号为c的栅格单元的汇流时间；MaxT为流域中各个栅格单元的汇流时间中的最大值；F为流域中栅格单元的总数；c为流域中栅格单元的编号，从1到F；L为指示函数，当τ_c∈[Z×Δtime,(Z+1)×Δtime)时L为1，否则L为0；

步骤5.2，将步骤5.1中计算得到的流域中的栅格数—时间关系转换为面积—时间关系；

将汇流时间在[Z×Δtime,(Z+1)×Δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围Tarea_z，即表示落在Tarea_z范围内的降雨会在[Z×Δtime,(Z+1)×Δtime)时间范围运动至流域的出口处；

Tarea_z＝Tcount_z×Csize²

即得到流域中的面积—时间：

步骤5.3，将步骤5.2中得到的面积—时间关系转换成流量—时间关系；

式中，UnitR为单位降雨深，一般取1mm；Tflow_z为Tarea_z范围内的单位降雨于[Z×Δtime,(Z+1)×Δtime)时间范围在流域出口处形成的径流；

即得到流域内的流量—时间关系：

步骤5.4，在流量—时间关系的基础上结合线性水库调蓄演算公式得到流域的地貌单位线；

UH_tu＝Y×Tflow_tu+(1-Y)×UH_tu-1

式中，tu为单位线时间序列编号、由1到TU，TU为UH_tu的值小于0.01的时刻；UH_tu为tu时的单位线纵坐标值；UH_tu-1为tu-1时单位线纵坐标值；Tflow_tu的值由tu与0、SZ以及TU的相对大小决定，Y为调蓄系数，F为预设的线性水库的演算系数。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤6的具体步骤如下：

步骤6.1，通过实地洪痕调查或者查阅流域管理机构的防灾预警规划得到安全保护对象所在位置的警戒流量Warn_Q；

步骤6.2，以逐次历史雨洪过程中的降雨数据为输入，采用新安江模型计算出流域出口的流量变化过程，确定最高流量Max_Q；

步骤6.3，若洪水过程中的Max_Q≥Warn_Q，则判断该场洪水为超警洪水，将该场洪水过程对应的数据点划分为超警子集A中，标记为1；否则将其归于未超警子集B中，标记为-1；

步骤6.4，基于支持向量机原理求解两类子集之间的最优分类决策方程：

令分类决策方程为W^T×X+γ＝0，式中，W^T为法向量，决定了分类决策方程所在的平面即决策平面的方向，γ为位置向量，决定了决策平面的位置，X为特征向量，其形式为[ACP,K,S,G]；

基于A和B两类子集与决策平面之间的位置关系构建约束条件：

式中，是决策平面和数据点x_i之间的距离，x_i表示为[U_i,P_i,K_i,S_i]，y_i为数据点x_i的标记，在A中为1，在B中为-1；

求解方程：

得到最优的法向量W^T _per和位置向量γ_per，从而得到最优的分类决策方程Per_Function：

Per_Function＝W^T _per×X+γ_per。

作为本发明基于多维降雨特征空间最优决策的山洪快速预警预报方法进一步的优化方案，所述步骤7的具体步骤如下：

步骤7.1，以实时降雨过程为输入，分别计算得到实时降雨过程的的降雨量H_s、降雨空间分布特征K_s、土壤水饱和度S_s和降雨时程分布特征G_s；

步骤7.2，以实时降雨过程的H_s、K_s、S_s和G_s为基础构建实时降雨的特征向量X_s＝[H_s，K_s，S_s,G_s]，并将之带入Per_Function中，若计算得到的数值大于1，则判定该次实时降雨需要发布预警，否则不发布预警；

步骤7.3，若实时降雨需要发布预警，则以该次实时降雨数据为输入，计算出地貌单位线后进行卷积计算，得到流域出口的流量过程，并发布预警。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明提供的一种基于多维降雨特征空间最优决策的山洪快速预警预报，以影响临界雨量的物理因子为基础，量化了土壤饱和度，降雨量，降雨空间分布，降雨时程分布对于临界雨量的影响作用，进而提出一种基于多维降雨特征空间最优决策的山洪快速预警预报。这样既保证了计算结果的精度与可靠性，同时解决了地理特征空间分异明显的山区性中小流域中的山洪快速预警预报问题。且本方法主要应用流域数字高程模型，数据来源稳定可靠，方法中变量之间的函数关系明确，有利于流域中山洪预警预报的快速自动发布，通过数字流域技术以简化提取步骤，同时，保证了结果的客观合理性，可以进一步促进数字水文学以及山区性中小流域山洪防治研究的深入发展。

附图说明

图1是计算流程示意；

图2是大河坝流域高程分布栅格；

图3是大河坝流域坡度栅格；

图4是大河坝流域坡向栅格；

图5是大河坝流域汇流累计栅格；

图6是大河坝流域河道栅格和坡面栅格；

图7是大河坝流域流速分布；

图8是大河坝流域流向指示示意图；

图9是大河坝流域汇流时间分布；

图10是大河坝流域地貌单位线；

图11是分类决策示意图。

具体实施方式

下面以大河坝流域为例，结合附图对本发明的技术方案做进一步的详细说明：

如图1所示，本发明公开了一种基于多维降雨特征空间最优决策的山洪快速预警预报方法，包括以下步骤：

步骤1，基于三层蒸散发原理计算逐次历史雨洪过程的初始土壤水饱和度S；

步骤1.1，整理得到逐小时不间断历史降雨径流资料，并依据径流过程的涨落均匀摘录出若干次雨洪过程；

步骤1.2，以逐小时不间断降雨资料为输入，采用三层蒸散发原理计算得到逐小时不间断的土壤水饱和度变化过程，并从中提取出摘录出的若干次雨洪过程中开始时的土壤水饱和度，记为S。

步骤2，统计逐次历史雨洪过程的降雨量H，其中，根据以下公式计算每次历史雨洪过程的降雨量H：

式中，T为雨洪过程的总时段数，t为时段的编码，从1到T；P_t为t时段内降雨。

步骤3，如图2所示，结合流域的数字高程模型数据、即DEM数据，量化逐次历史雨洪过程的降雨空间分布特征K；

步骤3.1，计算流域中每一个栅格单元中的汇流流速；

如图3、图4、图5所示，利用流域DEM数据将流域划分为若干个正交的栅格单元，计算出流域中每一个栅格单元的坡度、流向和汇流累计值，并提取流域中的河道栅格和坡地栅格，如图6所示；其中，根据以下公式计算栅格单元中的汇流流速v；

式中，D_p为坡面栅格中的坡度；θ为预设的第一经验系数阈值；A为河道栅格单元中的累计汇流面积，D_r为河道栅格单元中的坡度；a、ε、β分别为预设的第二、第三、第四经验系数阈值；

步骤3.2，如图7所示，根据每一个栅格单元中的汇流流速，得到各个栅格单元中的水流运动到流域出口所需要的时间，得到流域中的汇流时间栅格Raster_Htime；

步骤3.2.1，根据以下公式计算水流经过每一个栅格单元中所需要的时间Ctime：

式中，Csize为栅格单元的边长，即流域DEM数据的分辨率；dir为栅格单元的流向，即该栅格单元中的水流流动的方向，dir等于1、2、3、4、5、6、7、8时栅格单元的流向分别为东，东南，南，西南，西，西北，北，东北，如图8所示；

步骤3.2.2，根据以下公式计算编号为i的栅格单元中的水流运动到流域出口所需要的时间Htime_i：

式中，N为各个栅格单元中的水流运动到流域出口所经过的栅格单元的总数；i为各个栅格单元中的水流运动到流域出口所经过的栅格单元的编号，从1到N；Ctime_i为水流经过编号为i的栅格单元所需要的时间；

步骤3.2.3，如图9所示，将流域中所有栅格单元的Htime_i组合得到汇流时间栅格Raster_Htime；

步骤3.3，对逐次历史雨洪过程中站点观测的雨量数据进行空间插值，得到降雨空间分布栅格，并结合汇流时间栅格得到降雨空间分布特征K：

式中，F为流域中栅格单元的总数，c为流域中栅格单元的编码，从1到F；SP_c为编号为c的栅格单元中的降雨量，由逐次历史雨洪过程中站点观测的雨量数据进行空间插值得到，MaxP为流域中所有的栅格单元中最大的降雨量；Htime_c为编号为c的栅格单元中的汇流时间，MaxT为流域中所有的栅格单元中最大的汇流时间。

步骤4，量化逐次历史雨洪过程中的降雨时程分布特征G；

步骤4.1，将逐次历史雨洪过程中的降雨数据由小到大排序，得到NSeq_P；

NSeq_P＝{P₀,P₁…P_t…P_T}

式中，P_t表示t时段内的降雨量，t从0变化到T；且P₀<P₁<…<P_t<…P_T，即随着下标t的增大，P_t的值也在不断增大；

步骤4.2，基于NSeq_P计算得到累积时段雨量序列ASeq_P；

ASeq_P＝{AP₀,AP₁…AP_t…AP_T}

式中，AP_t为时段为t时的累积时段雨量，且AP₀<AP₁<…<AP_t<…AP_T，即随着下标t的增大，AP_t的值也在不断增大，最大的AP_T的值为H；

步骤4.3，计算Seq_P中的累积时段雨量百分比，得到累积时段雨量百分比序列ANSeq_P；

NSeq_P＝{APF₀,APF₁…APF_t…APF_T}

式中，t从0变化到T；APF₀<APF₁<…<APF_t<…APF_T，即随着下标t的增大，APF_t的值也在不断增大，最大的APF_T的值为1；APF_t的计算公式如下：

步骤4.4，基于ASeq_P和ANSeq_P计算G；

步骤5，如图10所示，基于DEM数据计算流域的地貌单位线；

步骤5.1，基于步骤3.2中计算得到的汇流时间栅格统计得到流域中的栅格数—时间关系；

在Raster_Htime的基础上，以Δtime为间隔进行分组统计，得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到流域中的栅格数—时间关系；

式中，Tcount_z表示汇流时间大于等于Z×Δtime、小于(Z+1)×Δtime的栅格单元的数量；Z为汇流时间分组的编号，从0到SZ，其中SZ＝MaxT/Δtime；Δtime为预设的时间间隔阈值；τ_c为编号为c的栅格单元的汇流时间；MaxT为流域中各个栅格单元的汇流时间中的最大值；F为流域中栅格单元的总数；c为流域中栅格单元的编号，从1到F；L为指示函数，当τ_c∈[Z×Δtime,(Z+1)×Δtime)时L为1，否则L为0；

步骤5.2，将步骤5.1中计算得到的流域中的栅格数—时间关系转换为面积—时间关系；

将汇流时间在[Z×Δtime,(Z+1)×Δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围Tarea_z，即表示落在Tarea_z范围内的降雨会在[Z×Δtime,(Z+1)×Δtime)时间范围运动至流域的出口处；

Tarea_z＝Tcount_z×Csize²

即得到流域中的面积—时间：

步骤5.3，将步骤5.2中得到的面积—时间关系转换成流量—时间关系；

式中，UnitR为单位降雨深，一般取1mm；Tflow_z为Tarea_z范围内的单位降雨于[Z×Δtime,(Z+1)×Δtime)时间范围在流域出口处形成的径流；

即得到流域内的流量—时间关系：

步骤5.4，在流量—时间关系的基础上结合线性水库调蓄演算公式得到流域的地貌单位线；

UH_tu＝Y×Tflow_tu+(1-Y)×UH_tu-1

式中，tu为单位线时间序列编号、由1到TU，TU为UH_tu的值小于0.01的时刻；UH_tu为tu时的单位线纵坐标值；UH_tu-1为tu-1时单位线纵坐标值；Tflow_tu的值由tu与0、SZ以及TU的相对大小决定，Y为调蓄系数，F为预设的线性水库的演算系数。

步骤6，如图11所示，以H，K，S和G为基础构建多维降雨特征空间坐标系，并构建降雨数据点集，依据洪水过程是否超警，将对应的降雨数据点集划分为超警子集A与未超警子集B两类，结合向量机原理即SVM原理求解两类子集之间的最优分类决策方程；

步骤6.1，通过实地洪痕调查或者查阅流域管理机构的防灾预警规划得到安全保护对象所在位置的警戒流量Warn_Q；

步骤6.2，以逐次历史雨洪过程中的降雨数据为输入，采用新安江模型计算出流域出口的流量变化过程，确定最高流量Max_Q；

步骤6.3，若洪水过程中的Max_Q≥Warn_Q，则判断该场洪水为超警洪水，将该场洪水过程对应的数据点划分为超警子集A中，标记为1；否则将其归于未超警子集B中，标记为-1；

步骤6.4，基于支持向量机原理求解两类子集之间的最优分类决策方程：

令分类决策方程为W^T×X+γ＝0，式中，W^T为法向量，决定了分类决策方程所在的平面即决策平面的方向，γ为位置向量，决定了决策平面的位置，X为特征向量，其形式为[ACP,K,S,G]；

基于A和B两类子集与决策平面之间的位置关系构建约束条件：

式中，是决策平面和数据点x_i之间的距离，x_i表示为[U_i,P_i,K_i,S_i]，y_i为数据点x_i的标记，在A中为1，在B中为-1；

求解方程：

得到最优的法向量W^T _per和位置向量γ_per，从而得到最优的分类决策方程Per_Function：

Per_Function＝W^T _per×X+γ_per。

步骤7，计算实时降雨过程的降雨量、降雨空间分布特征、土壤水饱和度和降雨时程分布特征，并基于最优分类决策方程判断该实时降雨过程是否需要发布预警，若需要进行预警则以该降雨为地貌单位线的输入，计算流域出口处的径流过程，进行实时洪水预报：

步骤7.1，以实时降雨过程为输入，分别计算得到实时降雨过程的的降雨量H_s、降雨空间分布特征K_s、土壤水饱和度S_s和降雨时程分布特征G_s；

步骤7.2，以实时降雨过程的H_s、K_s、S_s和G_s为基础构建实时降雨的特征向量X_s＝[H_s，K_s，S_s,G_s]，并将之带入Per_Function中，若计算得到的数值大于1，则判定该次实时降雨需要发布预警，否则不发布预警；

步骤7.3，若实时降雨需要发布预警，则以该次实时降雨数据为输入，计算出地貌单位线后进行卷积计算，得到流域出口的流量过程，并发布预警。

本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，包括以下步骤：

步骤1，基于三层蒸散发原理计算逐次历史雨洪过程的初始土壤水饱和度S；

步骤2，统计逐次历史雨洪过程的降雨量H；

步骤3，结合流域的数字高程模型数据、即DEM数据，量化逐次历史雨洪过程的降雨空间分布特征K；

步骤4，量化逐次历史雨洪过程中的降雨时程分布特征G；

步骤5，基于DEM数据计算流域的地貌单位线；

步骤6，以H，K，S和G为基础构建多维降雨特征空间坐标系，并构建降雨数据点集，依据洪水过程是否超警，将对应的降雨数据点集划分为超警子集A与未超警子集B两类，结合向量机原理即SVM原理求解两类子集之间的最优分类决策方程；

步骤7，计算实时降雨过程的降雨量、降雨空间分布特征、土壤水饱和度和降雨时程分布特征，并基于最优分类决策方程判断该实时降雨过程是否需要发布预警，若需要进行预警则以该降雨为地貌单位线的输入，计算流域出口处的径流过程，进行实时洪水预报。

2.根据权利要求1所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤1的具体步骤如下：

步骤1.1，整理得到逐小时不间断历史降雨径流资料，并依据径流过程的涨落均匀摘录出若干次雨洪过程；

步骤1.2，以逐小时不间断降雨资料为输入，采用三层蒸散发原理计算得到逐小时不间断的土壤水饱和度变化过程，并从中提取出摘录出的若干次雨洪过程中开始时的土壤水饱和度，记为S。

3.根据权利要求2所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤2中根据以下公式计算每次历史雨洪过程的降雨量H：

式中，T为雨洪过程的总时段数，t为时段的编码，从1到T；P_t为t时段内降雨。

4.根据权利要求3所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤3的具体步骤如下：

步骤3.1，计算流域中每一个栅格单元中的汇流流速；

利用流域DEM数据将流域划分为若干个正交的栅格单元，计算出流域中每一个栅格单元的坡度、流向和汇流累计值，并提取流域中的河道栅格和坡地栅格；其中，根据以下公式计算栅格单元中的汇流流速v；

式中，D_p为坡面栅格中的坡度；θ为预设的第一经验系数阈值；A为河道栅格单元中的累计汇流面积，D_r为河道栅格单元中的坡度；a、ε、β分别为预设的第二、第三、第四经验系数阈值；

步骤3.2，根据每一个栅格单元中的汇流流速，得到各个栅格单元中的水流运动到流域出口所需要的时间，得到流域中的汇流时间栅格Raster_Htime；

步骤3.2.1，根据以下公式计算水流经过每一个栅格单元中所需要的时间Ctime：

式中，Csize为栅格单元的边长，即流域DEM数据的分辨率；dir为栅格单元的流向，即该栅格单元中的水流流动的方向，dir等于1、2、3、4、5、6、7、8时栅格单元的流向分别为东，东南，南，西南，西，西北，北，东北；

步骤3.2.2，根据以下公式计算编号为i的栅格单元中的水流运动到流域出口所需要的时间Htime_i：

式中，N为各个栅格单元中的水流运动到流域出口所经过的栅格单元的总数；i为各个栅格单元中的水流运动到流域出口所经过的栅格单元的编号，从1到N；Ctime_i为水流经过编号为i的栅格单元所需要的时间；

步骤3.2.3，将流域中所有栅格单元的Htime_i组合得到汇流时间栅格Raster_Htime；

步骤3.3，对逐次历史雨洪过程中站点观测的雨量数据进行空间插值，得到降雨空间分布栅格，并结合汇流时间栅格得到降雨空间分布特征K：

式中，F为流域中栅格单元的总数，c为流域中栅格单元的编码，从1到F；SP_c为编号为c的栅格单元中的降雨量，由逐次历史雨洪过程中站点观测的雨量数据进行空间插值得到，MaxP为流域中所有的栅格单元中最大的降雨量；Htime_c为编号为c的栅格单元中的汇流时间，MaxT为流域中所有的栅格单元中最大的汇流时间。

5.根据权利要求4所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤4的具体步骤如下：

步骤4.1，将逐次历史雨洪过程中的降雨数据由小到大排序，得到NSeq_P；

NSeq_P＝{P₀,P₁…P_t…P_T}

式中，P_t表示t时段内的降雨量，t从0变化到T；且P₀<P₁<…<P_t<…P_T，即随着下标t的增大，P_t的值也在不断增大；

步骤4.2，基于NSeq_P计算得到累积时段雨量序列ASeq_P；

ASeq_P＝{AP₀,AP₁…AP_t…AP_T}

式中，AP_t为时段为t时的累积时段雨量，且AP₀<AP₁<…<AP_t<…AP_T，即随着下标t的增大，AP_t的值也在不断增大，最大的AP_T的值为H；

步骤4.3，计算Seq_P中的累积时段雨量百分比，得到累积时段雨量百分比序列ANSeq_P；

NSeq_P＝{APF₀,APF₁…APF_t…APF_T}

式中，t从0变化到T；APF₀<APF₁<…<APF_t<…APF_T，即随着下标t的增大，APF_t的值也在不断增大，最大的APF_T的值为1；APF_t的计算公式如下：

步骤4.4，基于ASeq_P和ANSeq_P计算G；

6.根据权利要求5所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤5的具体步骤如下：

步骤5.1，基于步骤3.2中计算得到的汇流时间栅格统计得到流域中的栅格数—时间关系；

在Raster_Htime的基础上，以Δtime为间隔进行分组统计，得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到流域中的栅格数—时间关系；

式中，Tcount_z表示汇流时间大于等于Z×Δtime、小于(Z+1)×Δtime的栅格单元的数量；Z为汇流时间分组的编号，从0到SZ，其中SZ＝MaxT/Δtime；Δtime为预设的时间间隔阈值；τ_c为编号为c的栅格单元的汇流时间；MaxT为流域中各个栅格单元的汇流时间中的最大值；F为流域中栅格单元的总数；c为流域中栅格单元的编号，从1到F；L为指示函数，当τ_c∈[Z×Δtime,(Z+1)×Δtime)时L为1，否则L为0；

步骤5.2，将步骤5.1中计算得到的流域中的栅格数—时间关系转换为面积—时间关系；

将汇流时间在[Z×Δtime,(Z+1)×Δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围Tarea_z，即表示落在Tarea_z范围内的降雨会在[Z×Δtime,(Z+1)×Δtime)时间范围运动至流域的出口处；

Tarea_z＝Tcount_z×Csize²

即得到流域中的面积—时间：

步骤5.3，将步骤5.2中得到的面积—时间关系转换成流量—时间关系；

式中，UnitR为单位降雨深，一般取1mm；Tflow_z为Tarea_z范围内的单位降雨于[Z×Δtime,(Z+1)×Δtime)时间范围在流域出口处形成的径流；

即得到流域内的流量—时间关系：

步骤5.4，在流量—时间关系的基础上结合线性水库调蓄演算公式得到流域的地貌单位线；

UH_tu＝Y×Tflow_tu+(1-Y)×UH_tu-1

式中，tu为单位线时间序列编号，由1到TU，TU为UH_tu的值小于0.01的时刻；UH_tu为tu时的单位线纵坐标值；UH_tu-1为tu-1时单位线纵坐标值；Tflow_tu的值由tu与0、SZ以及TU的相对大小决定，Y为调蓄系数，F为预设的线性水库的演算系数。

7.根据权利要求6所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤6的具体步骤如下：

步骤6.1，通过实地洪痕调查或者查阅流域管理机构的防灾预警规划得到安全保护对象所在位置的警戒流量Warn_Q；

步骤6.2，以逐次历史雨洪过程中的降雨数据为输入，采用新安江模型计算出流域出口的流量变化过程，确定最高流量Max_Q；

步骤6.3，若洪水过程中的Max_Q≥Warn_Q，则判断该场洪水为超警洪水，将该场洪水过程对应的数据点划分为超警子集A中，标记为1；否则将其归于未超警子集B中，标记为-1；

步骤6.4，基于支持向量机原理求解两类子集之间的最优分类决策方程：

令分类决策方程为W^T×X+γ＝0，式中，W^T为法向量，决定了分类决策方程所在的平面即决策平面的方向，γ为位置向量，决定了决策平面的位置，X为特征向量，其形式为[ACP,K,S,G]；

基于A和B两类子集与决策平面之间的位置关系构建约束条件：

式中，是决策平面和数据点x_i之间的距离，x_i表示为[U_i,P_i,K_i,S_i]，yi为数据点x_i的标记，在A中为1，在B中为-1；

求解方程：

得到最优的法向量W^T _per和位置向量γ_per，从而得到最优的分类决策方程Per_Function：

Per_Function＝W^T _per×X+γ_per。

8.根据权利要求7所述的基于多维降雨特征空间最优决策的山洪快速预警预报方法，其特征在于，所述步骤7的具体步骤如下：

步骤7.1，以实时降雨过程为输入，分别计算得到实时降雨过程的的降雨量H_s、降雨空间分布特征K_s、土壤水饱和度S_s和降雨时程分布特征G_s；

步骤7.2，以实时降雨过程的H_s、K_s、S_s和G_s为基础构建实时降雨的特征向量X_s＝[H_s，K_s，S_s,G_s]，并将之带入Per_Function中，若计算得到的数值大于1，则判定该次实时降雨需要发布预警，否则不发布预警；

步骤7.3，若实时降雨需要发布预警，则以该次实时降雨数据为输入，计算出地貌单位线后进行卷积计算，得到流域出口的流量过程，并发布预警。