技术领域
[0001] 本
发明属于洪水预报领域，更具体地，涉及一种基于广义第二类贝塔分布(Generalized beta distribution of the second kind，GB2)的洪水频率分析方法。
背景技术
[0002] 洪水一般用洪水发生的概率
密度定义洪水大小，如20年一遇、50年一遇和百年一遇等。洪水频率分析的目就是通过频率曲线的
外延，推求T年一遇洪水的设计洪水值。如计算的设计洪水值较大，则规模过大，会增加投资，造成浪费；如计算的设计洪水值较低，则规模过小，又可能在不利水文条件下导致工程失事造成损失。
[0003] 因此，开展高
精度的洪水频率分析是水利
水电工程设计与规划的首要问题，选择合适的频率分布线型与参数估计方法是其重要内容。
[0004] 目前多数研究仅选取单一的频率分布线型，如指数分布(Exponential，EXP)、Weibull分布、Gamma分布、Gumbel分布、广义极值分布(GEV)、皮尔逊III型分布(P-III)、对数皮尔逊III型分布(LP-III)及对数正态分布(LN)等分布进行洪水频率分析，其设计洪水结果具有较大的不确定性，设计洪水值的高估与低估会导致投资过多或安全
风险增大等后果。
发明内容
[0005] 针对
现有技术的以上
缺陷或改进需求，本发明提供了一种基于广义第二类贝塔分布(GB2)的洪水频率分析方法，其目的在于将广义第二类贝塔分布引入洪水频率分析中，并依据分布函数特点采用最大熵原理实现分布函数的参数估计，在此
基础上开展基于广义第二类贝塔分布的设计洪水计算，由此解决现有分析技术的洪水频率分析精度不高的技术问题。
[0006] 为实现上述目的，按照本发明的一个方面，提供了一种基于广义第二类贝塔分布的洪水频率分析方法，该方法包括以下步骤：
[0007] (1)采用年最大取样方法,采集年最大流量样本序列；
[0008] (2)应用广义第二类贝塔分布建立洪水概率密度模型：
[0009]
[0010] 其中，变量x表示年最大洪水流量，f(x)表示流量为x的概率密度；B()为贝塔函数；r1、r2、r3是形状参数，且r1>0,r2>0,r3＞0，r3决定全局形状，r1控制左尾形状，r2控制右尾形状，参数r1和r2共同决定分布的偏态特性；b为
位置参数，b＞0；
[0011] (3)推求基于最大熵原理的广义第二类贝塔分布的概率密度函数参数值，包括以下子步骤：
[0012] (31)引入最大熵原理，最大熵原理指出，某一随机变量x的概率密度函数f(x)可通过最大化熵值得到：
[0013]
[0014] 其中， gi(x)为x的函数；Ci为gi(x)的期望；
[0015] 依据拉格朗日乘子法，f(x)可表示为：
[0016] f(x)＝exp(-λ0-λ1g1(x)-λ2g2(x)…λmgm(x)) (3)
[0017] 其中：m为约束个数；λi,i＝0,1,2,...,m为拉格朗日乘子；
[0018] (32)基于最大熵原理，广义第二类贝塔分布约束条件表达式为：
[0019]
[0020]
[0021]
[0022] 其中，q为构造的约束条件参数；E为期望；
[0023] 则广义第二类贝塔分布的概率密度函数可构造为：
[0024]
[0025] 其中，λ2'为构造的拉格朗日乘子；
[0026] (33)推导拉格朗日乘子与广义第二类贝塔分布函数的约束条件之间关系，将式(5)带入到(4a)中可得：
[0027]
[0028] 令p＝b-q、 和t＝pxq，则上式可化简为：
[0029]
[0030] 令 则 因为y(0)＝0且y(∞)＝1，y∈[0,1]，则：
[0031]
[0032] 即有下式成立：
[0033]
[0034] 同时，计算λ0的另一种方法为：
[0035]
[0036] 令 随后对式(9)、(10)中的λ1和λ2分别求导，可得：
[0037]
[0038] 式中： 为digamma函数；
[0039]
[0040] 则拉格朗日乘子与广义第二类贝塔分布函数约束条件之间满足：
[0041]
[0042]
[0043] 对式(11)求二阶导可得下式：
[0044]
[0045]
[0046] (34)推导拉格朗日乘子与广义第二类贝塔分布函数参数之间关系，将式(8)带入到式(6)中，并与式(1)作比较，则下式成立：
[0047]
[0048] (35)推求广义第二类贝塔分布函数参数与约束条件之间关系，将式(16)代入式(13)和式(14)中，可得广义第二类贝塔分布函数参数与约束条件之间满足下式：
[0049]
[0050] 式中，var为方差；将步骤(1)中年最大流量样本序列,即x序列代入公式(17)，求出参数r1、r2、r3和b；
[0051] (4)将参数r1、r2、r3和b代入公式(1)中，即可得到洪水概率密度函数，可利用该函数得到洪水为某一量级x时的概率值f(x)；还可利用洪水概率密度函数，推求概率f(x)为T年一遇的设计洪水值x。
[0052] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下技术特征及有益效果：
[0053] 四参数的广义第二类贝塔分布包含了指数分布、威布尔(Weibull)分布及伽
马(Gamma)分布，具有足够的灵活性模拟复杂多变的数据集，因此非常适于水文频率分析。本发明通过引入广义第二类beta分布，并利用最大熵原理推求估计了分布函数参数，实现了精度较高的洪水频率分析。其拟合效果优于水文中的其它常用分布；最大熵原理能够有效地估计分布函数的参数，得到精度高的洪水频率分布函数。
附图说明
[0055] 图2为本发明方法拟合的洪水序列边缘分布图及概率密度曲线效果图。
具体实施方式
[0056] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及
实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0057] 如图1所示，本发明包括以下流程：
[0058] (1)采用年最大取样方法,采集年最大流量样本序列；
[0059] (2)应用广义第二类贝塔分布建立洪水概率密度模型：
[0060]
[0061] 其中，变量x表示年最大洪水值，f(x)表示流量为x的概率密度；B()为贝塔函数；r1、r2、r3是形状参数，且r1>0,r2>0,r3＞0，r3决定全局形状，r1控制左尾形状，r2控制右尾形状，参数r1和r2共同决定分布的偏态特性；b为位置参数，b＞0；
[0062] (3)基于最大熵原理，推导出广义第二类贝塔分布函数的参数计算公式：
[0063]
[0064]
[0065]
[0066]
[0067] 将步骤(1)中采集年最大流量样本序列数据作为序列x，代入以上公式，求出参数r1、r2、r3和b；
[0068] (4)将参数r1、r2、r3和b代入公式(1)中，得到洪水概率密度函数，可利用该函数输入需预测洪水流量x，计算得到洪水流量达到x的概率密度f(x)；还可利用洪水概率密度函数，输入洪水发生的概率密度f(x)，计算发生该概率下洪水流量x。
[0069] 实施例：现利用某流域典型水文
站点A、B和C的年最大日流量数据，检验广义第二类贝塔分布的模型效果。基于最大熵原理的方法，由A、B和C三个水文站点的年最大日流量序列求出模型的参数，表1给出了A、B和C三个水文站点洪水概率密度函数的参数，图2给出了基于GB2分布的3站年最大洪水序列边缘分布和概率密度曲线的拟合结果，结果显示GB2分布拟合效果良好，可用于设计洪水的计算。
[0070]
[0071] 表1
[0072] 采用GB2分布、正态分布、指数分布(EXP)、Gamma分布、Gumbel分布、广义正态分布(Generalized normal，GN)、P-III分布、广义Pareto分布(Generalized Pareto，GP)及Weibull分布等拟合典型水文站C的年最大日流量序列，并基于K-S(Kolmogorov-Smirnov，K-S)检验法、均方根误差以及AIC准则(Akaike information criterion，AIC)对各分部的拟合结果进行比较分析，确定拟合最优的分布函数，进行洪水频率分析。
[0073] 取K-S检验显著性水平为α＝0.05，P值大于0.05时通过检验，RMSE和AIC值越小说明拟合效果越好，表2给出了各分布K-S检验统计量P值、RMSE和AIC值。结果表明，所有分布均通过了K-S检验。比较各分部的RMSE和AIC值可知，GB2分布的拟合效果最优。
[0074]
[0075] 表2
[0076] 表3给出了基于上述分布函数计算的水文站点C的设计洪水值。结果表明，当重现期较大时，各分布计算结果差异明显，考虑到工程设计中，低估的设计值将大大增加大坝及下游的洪水风险，由表中数据比较分析可知GB2的频率分析结果更优。
[0077]
[0078]
[0079] 表3
[0080] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
