技术领域
[0001] 本
发明属于
水文学中的极端水文事件风险分析领域，更具体地，涉及一种基于最大熵的极端水文事件风险分析方法和系统。
背景技术
[0002] 极端水文事件风险分析是依据历史极值水文信息建立概率分布统计模型，揭示极端水文事件的
时空分布规律、演变趋势以及预估极值水文事件的重现期，反演一定设计标准的设计暴雨、设计洪水值。极端水文事件风险分析的核心是确定水文极值事件的分布函数形式。目前极端水文事件风险分析和设计洪水计算存在的主要问题是：受分布函数形式、偏态特性、尾部特征的限制，不同分布函数推求的概率风险和反演的设计洪水值差异明显。为此，本发明引入广义分布函数簇，参与计算。如何选择最优的分布函数是一个难点和热点问题。现有模型选择方法主要是依据均方根误差RMSE、AIC准则和BIC准则。但是，AIC准则和BIC准则受限于所选模型概率
密度函数参数的数目，对于参数较多的模型，不能给出相容估计，无法总是正确的选择最优分布。
[0003] 由此可见，
现有技术存在受限于所选模型概率密度函数参数数目、对于参数较多的模型不能给出相容估计、无法总是正确的选择最优分布的技术问题。
发明内容
[0004] 针对现有技术的以上
缺陷或改进需求，本发明提供了一种基于最大熵的极端水文事件风险分析方法和系统，由此解决现有技术存在受限于所选模型概率密度函数参数数目、对于参数较多的模型不能给出相容估计、无法总是正确的选择最优分布的技术问题。
[0005] 为实现上述目的，按照本发明的一个方面，提供了一种基于最大熵的极端水文事件风险分析方法，包括如下步骤：
[0006] (1)对采集的流域多年的实测流量值，使用年最大取样，得到长序列水文数据；
[0007] (2)将多种广义极值分布函数作为待选模型，对长序列水文数据进行拟合后对待选模型进行参数估计，得到待选模型的参数估计值，根据待选模型的参数估计值确定待选模型的概率密度分布函数；
[0008] (3)计算待选模型的概率密度分布函数的熵值，将熵值最大的待选模型作为最优模型；
[0009] (4)利用最优模型进行极端水文事件风险分析，得到流域百年一遇、千年一遇的设计洪水值。
[0010] 进一步地，多种广义极值分布函数包括：广义Gamma分布、广义第二类Beta分布、Halphen A分布、Halphen B分布和Halphen IB分布。
[0011] 进一步地，步骤(2)中进行参数估计的具体实现方式为：
[0012] 对长序列水文数据进行拟合后构建最大熵原理拉格朗日乘子与多种广义极值分布函数参数的函数映射关系，根据函数映射关系建立待选模型的参数估计方程组，利用待选模型的参数估计方程组计算待选模型的参数估计值。
[0013] 进一步地，待选模型的概率密度分布函数的熵值为：
[0014]
[0015] 其中，H(x)为连续信息熵的熵函数即待选模型的概率密度分布函数的熵值，f(x)为待选模型x的概率密度分布函数，a和b分别为计算熵值时的下限和上限。
[0016] 进一步地，步骤(4)包括：
[0017] 利用最优模型的概率密度分布函数和拟合后的长序列水文数据绘制概率流量曲线，利用概率流量曲线进行极端水文事件风险分析，将概率流量曲线中概率为1％时对应的流量作为流域百年一遇的设计洪水值，将概率流量曲线中概率为0.1％时对应的流量作为流域千年一遇的设计洪水值。
[0018] 按照本发明的另一方面，提供了一种基于最大熵的极端水文事件风险分析系统，包括：
[0019] 第一模
块，用于对采集的流域多年的实测流量值，使用年最大取样，得到长序列水文数据；
[0020] 第二模块，用于将多种广义极值分布函数作为待选模型，对长序列水文数据进行拟合后对待选模型进行参数估计，得到待选模型的参数估计值，根据待选模型的参数估计值确定待选模型的概率密度分布函数；
[0021] 第三模块，用于计算待选模型的概率密度分布函数的熵值，将熵值最大的待选模型作为最优模型；
[0022] 第四模块，用于利用最优模型进行极端水文事件风险分析，得到流域百年一遇、千年一遇的设计洪水值。
[0023] 进一步地，多种广义极值分布函数包括：广义Gamma分布、广义第二类Beta分布、Halphen A分布、Halphen B分布和Halphen IB分布。
[0024] 进一步地，第二模块中进行参数估计的具体实现方式为：
[0025] 对长序列水文数据进行拟合后构建最大熵原理拉格朗日乘子与多种广义极值分布函数参数的函数映射关系，根据函数映射关系建立待选模型的参数估计方程组，利用待选模型的参数估计方程组计算待选模型的参数估计值。
[0026] 进一步地，待选模型的概率密度分布函数的熵值为：
[0027]
[0028] 其中，H(x)为连续信息熵的熵函数即待选模型的概率密度分布函数的熵值，f(x)为待选模型x的概率密度分布函数，a和b分别为计算熵值时的下限和上限。
[0029] 进一步地，第四模块包括：
[0030] 利用最优模型的概率密度分布函数和拟合后的长序列水文数据绘制概率流量曲线，利用概率流量曲线进行极端水文事件风险分析，将概率流量曲线中概率为1％时对应的流量作为流域百年一遇的设计洪水值，将概率流量曲线中概率为0.1％时对应的流量作为流域千年一遇的设计洪水值。
[0031] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
[0032] (1)本发明通过引入多种广义极值分布、利用最大熵原理进行参数估计、模型选择，从分布形式-参数估计方法-模型选择三方面突破了极端水文事件风险分析的理论
瓶颈。提高了风险分析及设计洪水计算的精确度和合理性，有效降低了设计值过高造成的资源浪费以及设计值过低导致的破坏性风险。使得本发明不受限于所选模型概率密度函数参数数目、对于参数较多的模型可以给出相容估计、总是正确的选择最优分布。
[0033] (2)本发明作为新的型模选择方法，与现有方法进行相互验证，验证所提模型的合理性。突破了现有方法极值分布概率密度函数形状单一、
频率曲线
外延后设计结果偏差较大以及无法充分反映水文极值序列偏态、重尾特性的理论瓶颈。避免了工程设计中高估T年一遇设计值造成过于安全带来的经济损失，以及低估设计值导致的堤防水库溃堤漫坝风险，发挥了水库汛期的运行效益，提高了设计洪水计算合理性与可行性。
附图说明
[0034] 图1是本发明
实施例提供的一种基于最大熵的极端水文事件风险分析方法的
流程图；
[0035] 图2是本发明实施例提供的五种分布的概率流量曲线图。
具体实施方式
[0036] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0037] 如图1所示，一种基于最大熵的极端水文事件风险分析方法，包括如下步骤：
[0038] (1)对采集的流域多年的实测流量值，使用年最大取样，得到长序列水文数据；
[0039] (2)将多种广义极值分布函数作为待选模型，对长序列水文数据进行拟合后对待选模型进行参数估计，得到待选模型的参数估计值，根据待选模型的参数估计值确定待选模型的概率密度分布函数；
[0040] (3)计算待选模型的概率密度分布函数的熵值，将熵值最大的待选模型作为最优模型；
[0041] (4)利用最优模型进行极端水文事件风险分析，得到流域百年一遇、千年一遇的设计洪水值。
[0042] 进一步地，多种广义极值分布函数包括：广义Gamma分布、广义第二类Beta分布、Halphen A分布、Halphen B分布和Halphen IB分布。5种广义极值分布函数的适应性强，适用于极端水文事件的风险分析。
[0043] 进一步地，步骤(2)中进行参数估计的具体实现方式为：
[0044] 对长序列水文数据进行拟合后构建最大熵原理拉格朗日乘子与多种广义极值分布函数参数的函数映射关系，根据函数映射关系建立待选模型的参数估计方程组，利用待选模型的参数估计方程组计算待选模型的参数估计值。
[0045] 进一步地，待选模型的概率密度分布函数的熵值为：
[0046]
[0047] 其中，H(x)为连续信息熵的熵函数即待选模型的概率密度分布函数的熵值，f(x)为待选模型x的概率密度分布函数，a和b分别为计算熵值时的下限和上限。
[0048] 进一步地，步骤(4)包括：
[0049] 利用最优模型的概率密度分布函数和拟合后的长序列水文数据绘制概率流量曲线，利用概率流量曲线进行极端水文事件风险分析，将概率流量曲线中概率为1％时对应的流量作为流域百年一遇的设计洪水值，将概率流量曲线中概率为0.1％时对应的流量作为流域千年一遇的设计洪水值。
[0050] 选择熵值最大的模型作为最优模型的依据为：
[0051] 根据最大熵原理，在没有约束条件情况下，数据集中的每个数据出现的概率都是相等的，此时函数的熵值为最大，随着约束条件的增加，熵值会不断减小。下面通过离散信息熵和连续信息熵进行验证：
[0052] 利用离散信息熵的表达式：
[0053]
[0054] 式中，H(p)为离散信息熵的熵函数；pi为每个事件发生的概率。熵的单位与公式中对数的底有关。
[0055] 在没有其他约束条件的情况下，由：
[0056]
[0057] 用拉格朗日乘子法求条件极值，构造H(p)的拉格朗日乘子式L(pi)，求解H(p)的最大值如下：
[0058]
[0059] 式中，L(pi)为拉格朗日乘子式； 是约束条件；待定系数λ是拉格朗日乘子。
[0060] 分别对pi进行求导，并让求导后的值为0，求解出pi和λ，得到如下等式：
[0061]
[0062] 解得：
[0063] pi＝eλ-1
[0064] 由：
[0065]
[0066] 可得：
[0067]
[0068] 说明，在没有约束条件情况下，数据集中的每个数据出现的概率都是相等的，此时函数的熵值为最大。
[0069] 利用连续信息熵的表达式：
[0070]
[0071] 式中，H(x)为连续型信息熵的熵函数；f(x)为变量x的概率密度函数。假设概率分布f(x)满足条件：
[0072]
[0073] 式中，Ar(x)是关于x的任意函数；Cr是Ar(x)的期望值。
[0074] 通过上述离散熵方程可知，对于连续熵函数，同样用拉格朗日乘子法求条件极值，构造f(x)的拉格朗日乘子式，可求解出f(x)的表达式如下：
[0075]
[0076] 则：
[0077]
[0078] 将f(x)表达式带入Hm(f)中可得到：
[0079]
[0080] 式中，λr为约束确定的拉格朗日乘数；r表示数量。
[0081] 假设q(x)是满足n个约束条件的另一个概率分布，其中：n≥m
[0082]
[0083] 令：
[0084]
[0085] 构造函数：
[0086] g(x)＝xlnx-x+1
[0087] 式中，x＞0。
[0088] 对g(x)求导可得：
[0089] g′(x)＝lnx则g(x)在(0，1)上单调递减，在(1，∞)上单调递增。
[0090] 则：
[0091] g(x)≥g(1)＝0可得：
[0092] xlnx≥x-1
[0093] 上式可改写为：
[0094]
[0095] 其中，h(x)＞0。
[0096] 令： 可得：
[0097]
[0098] 等式两边同乘上q(x)，可得：
[0099]
[0100] 即：
[0101]
[0102] 则：
[0103]
[0104] 将q(x)与f(x)表达式带入上式，得：
[0105]
[0106] 其中：
[0107]
[0108]
[0109] 由于n≥m，则：λ0、λ1、...λm存在且有意义，则上式可写出：
[0110]
[0111] 可得：
[0112] Hn(q)≤Hm(f)，n≥m
[0113] 以上说明，随着约束条件的增加，熵值会不断减小。
[0114] 上述验证过程中，选取5种分布模型中熵值最大的模型最为最优模型，是因为此模型是在包含已知信息的情况下，约束条件增加最少的模型，同样也是进行未知假设最少的模型。与普拉斯不充分推理原则相符，也与最大熵原理相符。即：在包含已知信息的情况下，应当不做任何未知假设，把未知事件当成等概率事件处理。因此，理应为最优拟合模型。
[0115] 为了更清晰的展现本发明的目的和结构及技术方案，以下以清江流域结合附图对本发明进一步详细说明，具体实施步骤包括：
[0116] 步骤1，获取清江流域43年来的年径流资料，运用年最大值法计算出每年代表年径流量，得到长序列水文数据，作为计算拟合数据。
[0117] 步骤2，根据步骤1，选择广义Gamma分布、广义第二类Beta分布(GB2)、Halphen A分布、Halphen B分布以及Halphen IB分布5种广义分布模型，进行拟合。
[0118] 步骤3，根据步骤2，利用最大熵原理，获得各模型的估计参数，进而推导5种模型的概率密度表达式。
[0119] 步骤4，根据步骤3，利用连续信息熵的表达式：
[0120]
[0121] 将上述5种分布模型的概率密度函数分别带入连续信息熵公式中，求得5种分布模型的熵值，根据最大熵原理，选取5种分布模型中，熵值最大的模型作为最优模型。所得结果如表1。
[0122] 表1 最大熵原理计算结果表
[0123] 分布函数 熵值H(x)GG 3.435
GB2 3.668
Hal-A 2.827
Hal-B 2.442
Hal-IB 2.212
[0124] 步骤5，根据步骤4，计算5种分布的AIC值、BIC值、RMSE值，选出最优模型，与最大熵原理所选模型进行比较分析，得出最优模型。所得结果如表2、图2。
[0125] 表2 AIC、BIC、RMSE值计算结果表
[0126] 分布函数 RMSE AIC BICGG 0.0259 -256.8 -248.4
GB2 0.0196 -263.1 -254.7
Hal-A 0.0269 -255.5 -247.1
Hal-B 0.0288 -252.9 -244.5
Hal-IB 0.0379 -211.9 -203.5
[0127] 步骤6，根据步骤5，画出5种分布的概率-流量曲线，并结合最优模型分布，推求5种分布预测的清江流域百年一遇、千年一遇设计洪水值。所得结果如表3。
[0128] 表3 清江流域百年一遇、千年一遇设计洪水值表
[0129] 分布函数 百年一遇 千年一遇GB2 12630.92 19606.77
GG 11635.92 14959.42
Hal-A 11448.49 14480.81
Hal-B 11174.42 13945.87
Hal-IB 11322.21 14532.21
[0130] 根据表1所示，熵值最大的模型为GB2分布，作为最大熵原理选出的最优分布，根据表2所示，根据AIC值、BIC值、RMSE值所选模型为GB2分布，与最大熵原理结论相符合。可得出清江流域年径流量的最优分布为GB2分布。根据图2所示，GB2分布的概率-流量曲线明显更
符合事件情况。根据表3所示，清江流域百年一遇、千年一遇设计洪水值分别为12630.92m3/s、19606.77m3/s。可以看出，GB2预测值与其他分布预测值有明显差距，可以更好的拟合真实情况。
[0131] 综上所述，本发明引入多种广义极值分布函数作为待选模型；运用最大熵原理进行参数估计，确定分布函数的具体概率密度函数表达式；根据离散熵、连续熵公式，推导最大熵选择分布函数的原理，计算分布函数的熵值，利用熵值的大小来选择最优分布函数形式；进而推求流域百年一遇、千年一遇的设计洪水值。本发明提高了风险分析及设计洪水计算的精确度和合理性，有效降低了设计值过高造成的资源浪费以及设计值过低导致的破坏性风险，突破了现有方法极值分布概率密度函数形状单一、频率曲线外延后设计结果偏差较大以及无法充分反映水文极值序列偏态、重尾特性的理论瓶颈。
[0132] 本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
