阅读说明：本技术 基于熵决策法的动态临界雨量计算方法 (Entropy decision method based dynamic critical rainfall calculation method ) 是由 林凯荣 梁汝豪 兰甜 陈海燕 于 2019-07-12 设计创作，主要内容包括：本发明涉及一种基于熵决策法的动态临界雨量计算方法；S1.收集整理流域多场洪水的实测累积雨量和对应洪峰流量数据,将其以前期土壤湿润程度和历时划分为多个样本系列；S2.对收集到的样本系列进行边缘分布拟合,用K-S检验选出拟合效果最好的分布；S3.利用Frank Copula函数构造联合分布,得到不同前期土壤湿润程度及不同历时条件下累积雨量和对应洪峰流量的联合概率分布；S4.将联合概率分布及相应参数分别代入贝叶斯效用函数和期望效用-熵风险函数公式,分别考虑主观因素、客观因素及主客观因素同时影响的情况,通过使效用函数和风险最小化得到不同前期土壤湿润程度及不同历时条件下的临界雨量值。本发明用风险决策方法解决了山洪灾害预警指标的确定问题。(The invention relates to a dynamic critical rainfall calculation method based on an entropy decision method; s1, collecting and settling actually-measured accumulated rainfall and corresponding flood peak flow data of multi-field flood of a drainage basin, and dividing the actually-measured accumulated rainfall and the corresponding flood peak flow data into a plurality of sample series according to the soil wetting degree and duration of a previous period; s2, performing edge distribution fitting on the collected sample series, and selecting the distribution with the best fitting effect by using K-S inspection; s3, constructing combined distribution by using a Frank Copula function to obtain combined probability distribution of accumulated rainfall and corresponding flood peak flow under different early-stage soil wetting degrees and different duration conditions; and S4, substituting the joint probability distribution and corresponding parameters into a Bayes utility function and an expected utility-entropy risk function formula respectively, considering the conditions that subjective factors, objective factors and subjective and objective factors are simultaneously influenced, and minimizing the utility function and the risk to obtain critical rain values under different soil wetting degrees in the early stage and different duration conditions. The invention solves the problem of determining the mountain torrent disaster early warning index by using a risk decision method.)

基于熵决策法的动态临界雨量计算方法

技术领域

本发明涉及山洪灾害预警预报领域，更具体地，涉及一种基于熵决策法的动态临界雨量计算方法。

背景技术

山洪灾害是指由短时间暴雨引起的突发性强、破坏力大的山丘区洪水。近年来随着我国经济、社会的快速发展和极端天气的加剧，山洪灾害已成为我国最严重的自然灾害之一。据统计，我国每年因山洪灾害造成的直接经济损失约占全国洪涝灾害损失的70％，伤亡人数约占80％。我国大部分地区属于季风气候，降雨时段集中，加之我国山区面积广大，水土流失严重，使得山洪灾害频发，造成严重的人员伤亡和财产损失，也阻碍了山区的经济社会发展。

由于山洪灾害具有突发性强、破坏力大的特点，加之山区监测站点较少、水文气象资料短缺等问题，使得山洪灾害难以进行有效地预警和防治。通过确定山洪灾害预警指标的方法可以有效解决这一难题。该方法是山洪灾害防治的重要环节，在山洪到来前只需比较预警指标与预报数据即可确定是否发生山洪灾害。目前常用的预警指标是临界雨量。临界雨量是指当降雨量达到某一量级时流域内发生山洪灾害，该降雨量即为临界雨量。当降雨量达到临界雨量时，地表形成径流，致使河流某一断面的流量超过临界流量，洪水淹没河流两岸的房屋农田，造成一定的社会经济损失。对于中小流域，山洪突发性强，准确的临界雨量可以使人们预知山洪的到来，从而提前发布洪水警报，开展疏散转移工作，最大程度地减小山洪带来的损失。若临界雨量偏小，会导致人员过早地转移，造成转移成本的增加以及社会秩序的混乱。若临界雨量偏大，则会出现山洪灾害已经发生而群众物资还未转移的情况，造成严重的损失。因此，通过科学方法计算出准确可行的动态临界雨量，能最大限度地减小山区的人员财产损失，为山洪灾害评价、预警、防治和决策部署工作提供技术支撑，对于山洪灾害的预警预报和防灾减灾工作意义重大。

发明内容

本发明为克服上述现有技术中临界雨量计算不准确的缺陷，提供一种基于熵决策法的动态临界雨量计算方法，通过该方法确定临界雨量具有较强的合理性和准确性，采用该方法能够求得风险最小时的临界雨量，使决策者在面对即将到来的山洪时，通过比较实测雨量与临界雨量，做出风险最小的决定。

为解决上述技术问题，本发明采用的技术方案是：一种基于熵决策法的动态临界雨量计算方法，包括以下步骤：

S1.收集整理流域多场洪水的实测累积雨量和对应洪峰流量数据，将其以前期土壤湿润程度和历时划分为多个样本系列。

本发明收集整理了流域多场实测洪水的降雨径流数据，考虑到临界雨量与前期土壤湿润程度的相关性，通常的方法是根据前期土壤湿润程度(AMCⅠ，AMCⅡ，AMCⅢ)对洪水径流数据进行划分，然后在每一种前期土壤湿润程度条件下按不同历时(3h，6h，12h，24h)进行细分。若收集到的降雨径流数据较少，此时根据前期土壤湿润程度划分数据会产生非常短的数据序列，难以选定准确的边缘分布。因此需将各场洪水的降雨数据输入水文模型中，模拟出不同前期土壤湿润程度条件下的流量系列，再将其按不同历时细分，得到多组累积雨量v与对应洪峰流量q的样本系列。

S2.对收集到的样本系列进行边缘分布拟合，用K-S检验选出拟合效果最好的分布。

进一步的，本发明选取指数分布、对数正态分布以及威布尔分布对累积雨量v和对应洪峰流量q样本系列进行拟合，再用K-S检验选出拟合效果最好的分布，具体计算步骤为：

S21.选取指数分布、对数正态分布以及威布尔分布对不同前期土壤湿润程度和不同历时条件下的累积雨量v和对应洪峰流量q样本系列进行边缘分布拟合，同时求出各个分布对应的参数。其中，指数分布、对数正态分布以及威布尔分布的函数表达式分别如下所示：

式中，x是随机变量，λ为指数分布的比例参数；

式中，x是随机变量，μ和σ分别为y＝ln(x)的均值和标准差；

式中，x为随机变量，α为尺度参数，β为形状参数，γ为位置参数；

三种边缘分布的拟合可分别用MATLAB软件中的函数expfit，lognfit和wblfit实现。

S22.用K-S检验法判断样本是否服从某一种分布，比较三种分布的P值，P值是K-S检验法中用来衡量显著性水平的统计量。在本发明中原假设为样本服从理论分布函数F(x)，若在显著性水平下(α＝0.05)P值小于0.05，则拒绝原假设，该样本不服从给定分布函数F(x)。反之，则接受原假设，且P值越大拟合效果越好。最终选取拟合效果最好的分布作为该样本系列的最优边缘分布。K-S检验可采用MATLAB软件中的函数kstest实现。

S3.利用Frank Copula函数构造联合分布，得到不同前期土壤湿润程度及不同历时条件下累积雨量和对应洪峰流量的联合概率分布。

进一步的，步骤S3具体包括以下步骤：

S31.选取Frank Copula函数构造联合分布，其定义如下：

式中，F(x,y)为随机变量(X,Y)在(x,y)处的联合分布函数，C为Copula函数，θ是Copula函数的参数，可通过Kendall相关系数τ计算得到，τ的计算公式如下所示：

随机变量v和q的联合累积分布函数F(x,y)可定义为：

F(x,y)＝C(F(x),G(y))＝C(u,v)

式中，F(x,y)为随机变量(X,Y)在(x,y)处的联合分布函数，C为Copula函数；其联合密度函数如下所示：

f(x,y)＝c[F(x),G(y)]f(x)g(y)

其中c是C的密度函数，f(x)和g(y)分别是随机变量X和Y的概率密度函数；

通过样本系列求得Kendall相关系数τ，通过上式可反推得到参数θ，再结合两个随机变量v和q的边缘分布即可求出联合分布。

在MATLAB软件中，Kendall相关系数τ可用函数corr求出，Copula函数的参数可用函数copulaparam计算得到，copula函数的联合概率分布则可通过函数copulacdf计算得到。

S4.将联合概率分布及相应参数分别代入贝叶斯效用函数和期望效用-熵风险函数公式，分别考虑主观因素、客观因素及主客观因素同时影响的情况，通过使效用函数和风险最小化得到不同前期土壤湿润程度及不同历时条件下的临界雨量值。

进一步的，步骤S4具体步骤包括：

S41.期望效用-熵决策模型能够将实施行动的客观风险与决策者的主观偏好结合起来。其函数表达式如下：

式中，mean_a∈A{|E[(u(X(a,θ))]|}≠0，Ha(θ)表示行动a对应状态θ的熵；X(a,θ)表示采取行动a时对应状态θ的结果，由四部分组成：X₁₁(准确发出警报，q≥Q*且v≥V_T)、X₁₂(漏报，q≥Q*且v＜V_T)、X₂₁(误报，q＜Q*且v≥V_T)和X₂₂(准确不发出警报，q＜Q*且v＜V_T)；λ∈[0,1]为“权衡系数”，反映了决策者行为的主观期望效用与客观不确定性之间的权衡。当λ＝0时，仅考虑决策者的主观偏好，期望效用将产生更大的影响；当λ＝1时，则不考虑决策者的期望效用，仅考虑熵，即客观不确定性的影响。但在实际决策中，往往需要同时考虑决策者的主观期望和客观风险，因此假设λ＝0.5作为两种因素同时考虑的情况。当采取行动a使函数取得最小值即风险最小时，行动a即为最优行动方案，此时的V_T即为最优临界雨量；

S42.利用贝叶斯效用函数计算仅考虑决策者的主观偏好时的临界雨量，其定义如下：

式中，q是流量值，m³/s；Q*是河流临界断面的临界流量值，m³/s；v是累积雨量值，mm；V_T为临界雨量值，mm；T是暴雨历时(3h，6h，12h，24h)；参数a＝10×10⁶，b＝200，c＝0.025，C₀＝10×10³，a’＝5×10⁶，b’＝800，c’＝0.03；

对于不同前期土壤湿润程度(AMCⅠ、AMCⅡ、AMCⅢ)及不同降雨历时(3、6、12、24小时)，临界雨量V_T可通过最小化预期效用损失函数来确定，即找到使得预期效用损失函数最小的V_T值V_T*，该值即为所求值；具体公式如下所示：

式中，f(q,v|T)是累积降雨量和对应洪峰流量的联合概率密度，U(q,v|V_T,T)是效用函数值；

S43.利用熵计算仅考虑客观不确定性时的临界雨量，若有两个连续随机变量X和Y，其联合熵和条件熵可分别如下式所示：

式中，f(x,y)为随机变量X和Y的联合概率密度；

式中，f(x|y)为条件概率密度，即Y取任意一个定值的情况下X取值的概率密度；根据条件概率密度的性质可得：

在临界雨量计算中条件熵可表达为以下形式，其中v为累积雨量，q为累积雨量对应的洪峰流量：

熵决策法是投资决策领域的经典方法，能够将实施行动的客观风险与决策者的主观偏好结合起来。在山洪到来时，对于决策者而言，是否发出洪水警报就是一种风险行动。在采取行动时，决策者的主观判断和客观事件的不确定性共同决定了该行动的风险大小，对于决策者来说，风险总是越小越好的，采用该方法能够求得风险最小时的临界雨量，使决策者在面对即将到来的山洪时，通过比较实测雨量与临界雨量，做出风险最小的决定。因此采用这一方法确定临界雨量具有较强的合理性和准确性，在未来的山洪灾害预警预报研究中也具有较大的前景性。本发明用风险决策方法解决了山洪灾害预警指标的确定问题，既考虑了决策者主观偏好的影响，也考虑了客观状态的不确定性，充分体现了风险行动的特性，且该方法计算复杂度较低，结构清晰明确，计算结果合理准确，可广泛用于山洪灾害预警预报工作。

与现有技术相比，有益效果是：本发明提供的一种基于熵决策法的动态临界雨量计算方法，引入风险决策领域中的熵决策法计算得到不同前期土壤湿润程度条件下的动态临界雨量，能最大限度地减小山洪暴发造成的人员财产损失，为山洪灾害评价、预警、防治和决策部署工作提供科学依据，对于山洪灾害的预警预报和防灾减灾工作意义重大。

附图说明

图1为本发明方法整体流程图。

图2为本发明计算累积雨量v及其对应洪峰流量q的三个时间段(点)示意图。

图3为本发明不同前期土壤湿润程度和不同历时累积降雨量v与对应洪峰流量q样本系列散点图。

图4为本发明累积雨量v与对应洪峰流量q样本系列的边缘分布拟合效果对比图。

图5为本发明前期土壤湿润条件下不同历时联合概率密度图。

图6为本发明前期土壤湿润条件下不同历时联合概率分布图。

图7为本发明利用熵决策方法得到的临界雨量计算结果图。

图8为本发明熵决策法临界雨量应用效果评估图。

具体实施方式

附图仅用于示例性说明，不能理解为对本发明的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本发明的限制。

实施例1：

如图1所示，一种基于熵决策法的动态临界雨量计算方法，包括以下步骤：

步骤1：本发明选择广东省佛冈县潖江上游流域作为实验区域。广东省地处亚热带季风气候区，雨热同期，且位于热带风暴的运动路径上，易受暴雨极端天气侵袭，再加上广东山地、丘陵面积广大，河流密集，易造成山洪灾害。山区密集的人口使得土地过度开发，乱砍滥伐致使水土流失加剧，也加重了广东省山洪灾害的影响。广东省清远市佛冈县位于广东省中部，是广东省三大暴雨中心区之一，地形地势自东北向西南倾斜，境内低山、丘陵、平原交错，河流众多，极易发生山洪灾害。因此，为进一步做好该区域山洪灾害的预警预报及防治工作，本发明提出一种动态临界雨量计算方法，为该区域山洪灾害的预警预报防灾减灾工作提供技术支撑。

步骤2：首先本发明收集了潖江上游流域近20年发生的23场洪水发生时的降雨径流数据，将23场洪水的降雨数据输入至率定好的TOPMODEL水文模型中，模拟出不同前期土壤湿润程度的流量系列，再将其按不同历时(3h，6h，12h，24h)细分，这样对于不同前期土壤湿润程度以及不同历时都有23组累积雨量v及对应洪峰流量q样本系列。本发明对累积雨量和对应洪峰流量的定义如图2所示。图中的黑色实线表示一次假设洪水过程，t0为暴雨开始时间，Tr为降雨历时(3h，6h，12h，24h)，Tc为流域平均汇流时间。因此用于拟合边缘分布的累积雨量v为t0到t0+Tr的累积降雨量，对应洪峰流量q为t0到t0+Tr+Tc的时间间隔内出现的最大流量值。最终得到的样本系列散点图如图3所示，其中，各个图中，不同形状的散点代表前期不同的土壤湿润程度，自上往下，依次代表的是AMC1(干旱)，AMC2(正常)和AMC3(湿润)；。

步骤3：将各组累积雨量v与对应洪峰流量q样本系列分别用指数分布、对数正态分布以及威布尔分布进行拟合，并对拟合效果进行K-S检验，检验结果如图4所示。图中纵坐标为K-S检验P值，P值越大说明边缘分布拟合效果越好。不同前期土壤湿润程度和不同历时累积雨量v与对应洪峰流量q的最优边缘分布及其参数详见表1。

表1累积雨量v与对应洪峰流量q样本系列的最优边缘分布参数值

步骤4：在MATLAB软件中利用Frank Copula函数构造联合分布。v-q样本系列的Kendall相关系数及Copula函数参数如表2所示，得到的联合概率密度和联合概率分布结果分别如图5和图6所示(以前期土壤湿润的情况为例)。表2累积雨量v与对应洪峰流量q的Kendall相关系数以及Copula函数参数

步骤5：通过对贝叶斯效用函数及v-q联合概率密度的乘积进行积分可计算得到效用函数的期望，期望取得最小值时的V_T即为最优临界雨量。贝叶斯效用函数仅考虑决策者的主观判断，计算结果如表3和图7所示。

步骤6：由前文所述的期望效用-熵风险函数的定义可知，当采取行动a使函数取得最小值即风险最小时，行动a即为最优行动方案，此时的V_T即为最优临界雨量。当λ＝1时，R(a)＝Ha(θ)，此时仅考虑客观状态的不确定性，结合熵的公式即可确定临界雨量。在实际决策中，往往需要同时考虑决策者的主观期望和客观风险，因此本发明还假设λ＝0.5作为两种因素同时考虑的情况，计算结果如表3和图7所示。

表3熵决策法临界雨量计算结果单位：mm

步骤6：本发明从实验区域历史径流系列中选取前期土壤湿润程度分别为干旱、正常、湿润的三场洪水进行应用效果分析。将熵决策法计算出的临界雨量结果与实际洪水的累积雨量相比较，评估临界雨量计算结果的准确性，结果如图8所示(图中带标记的实线表示利用熵决策法得到的临界雨量，黑色虚线表示验证洪水的实测累积雨量，红色竖直虚线表示洪水流量超过临界流量的时刻。横坐标表示距降雨起始时刻的历时)。根据评估结果可知，熵决策法计算出的临界雨量在总体上能够对山洪灾害进行准确的预警，可为该区域山洪灾害的预警预报防灾减灾工作提供技术支撑。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。