CN106557614A - 一种基于Halphen IB分布的洪水频率分析方法 - Google Patents

Abstract

本发明公开了一种基于Halphen IB分布的洪水频率分析方法和系统，属于水文分析计算领域。本发明方法率先将Halphen IB分布引入洪水频率分析中，构建洪水概率密度函数，并依据分布函数特点采用最大熵原理实现Halphen IB分布的参数估计，进而利用洪水概率密度函数推求T年一遇设计洪水值。本发明还实现了一种基于Halphen IB分布的洪水频率分析系统，本发明的洪水概率分析方法和系统非常适合水文频率分析，其拟合效果基本优于水文频率分析中其他分布。本发明为水文频率分析提供了一种更加有效的途径。

Description

一种基于Halphen IB分布的洪水频率分析方法

技术领域

本发明属于洪水预报领域，更具体地，涉及一种基于Halphen IB分布的洪水频率分析方法和系统。

背景技术

洪水一般用洪水发生的概率密度定义洪水大小，如20年一遇、50年一遇和百年一遇等。洪水频率分析的目的就是通过频率曲线的外延，推求T年一遇洪水的设计洪水值。如计算的设计洪水值较大，则规模过大，会增加投资，造成浪费；如计算的设计洪水值较低，则规模过小，又可能在不利水文条件下导致工程失事造成损失。

因此，开展高精度的洪水频率分析是水利水电工程设计与规划的首要问题，选择合适的频率分布线型与参数估计方法是其重要内容，目前多数研究仅选取单一的频率分布线型，如指数分布(Exponential，EXP)、Weibull分布、Gamma分布、Gumbel分布、广义极值分布(GEV)、皮尔逊III型分布(P-III)、对数皮尔逊III型分布(LP-III)及对数正态分布(LN)等分布进行洪水频率分析，其设计洪水结果具有较大的不确定性，设计洪水值的高估与低估会导致投资过多或安全风险增大等后果。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于Halphen IB分布的洪水频率分析方法和系统，其目的在于将Halphen IB分布引入洪水频率分析中，并依据分布函数特点采用最大熵原理实现分布函数的参数估计，在此基础上开展基于Halphen IB分布的设计洪水计算，由此解决现有分析技术的洪水频率分析精度不高的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种基于Halphen IB分布的洪水频率分析方法，该方法包括以下步骤：

(1)对特定水文站进行取样,采集年最大流量样本序列即x序列；

(2)采用Halphen IB分布构建洪水概率密度函数：

Halphen IB型分布的概率密度函数表达式为：

其中，变量x表示年最大流量，f(x)表示流量为x的概率密度；α和v是形状参数，m是尺度参数，且m＞0，α＞0；α、m和v共同决定Halphen IB分布的形状；

(3)基于最大熵原理，推导出Halphen IB分布函数的参数，包括以下子步骤：

(31)依据最大熵原理，洪峰序列x的概率密度函数f(x)可通过最大化熵值得到：

其中，g_i(x)为x的函数；C_i为g_i(x)的期望；

依据拉格朗日乘子法，f(x)可表示为：

f(x)＝exp(-λ₀-λ₁g₁(x)-λ₂g₂(x)…λ_mg_m(x)) (3)

其中：m为约束个数；λ_i,i＝0,1,2,...,m为拉格朗日乘子；

(32)基于最大熵原理，构造拟合洪峰序列x的Halphen IB型分布约束条件表达式为：

其中，E为期望；则概率密度函数可构造为：

(33)为了得到Halphen IB型分布中λ₀的表达式，将式(5)带入到(4a)中可得：

因此，可得λ₀的表达式为：

为了得到Halphen IB型分布参数的约束等式，对式(8)中的λ₁、α分别求导，可得：

由于上述两式无法求解三个参数，故对α求二阶导，有：

式(9a)、(9b)、(10)即为拉格朗日乘子与约束间的关系；

(34)为得到Halphen IB型分布拉格朗日乘子与参数间的关系，将式(7)代入式(5)，得：

为了得到与原函数相似表达式令得：

对比Halphen IB型分布的原函数即公式(1)：

可得Halphen IB型分布拉格朗日乘子与参数间的关系为：

(35)最后可以推求得到Halphen IB型分布参数与约束关系,由指数阶乘函数的性质：

以及(13)，整理式(9a)、式(9b)和式(10)可得：

进一步可得：

其中，归一化函数ef_v(·)是Halphen提出的指数阶乘函数；E为期望；var()为方差；将步骤(1)中采集河流历年最大流量序列数据即x序列代入以上公式，求出洪水概率密度函数的参数α、m和v；

(4)对洪水概率密度函数进行积分，得到分布函数F(x)：

其中，P表示流量小于X的洪水发生概率；利用分布函数推求T年一遇的洪水流量设计值X，其中

按照本发明的另一方面，提供了一种基于Halphen IB分布的洪水频率分析系统，该系统包括以下模块：

采样模块，用于对特定水文站进行取样,采集年最大流量样本序列即x序列；

建模模块，用于采用Halphen IB分布构建洪水概率密度模型：

Halphen IB型分布的概率密度函数表达式为

其中，变量x表示年最大流量，f(x)表示流量为x的概率密度；α和v是形状参数，m是尺度参数，且m＞0，α＞0；α、m和v共同决定Halphen IB分布的形状；

参数推导模块，用于基于最大熵原理，推导出洪水概率密度函数的参数满足下式：

其中，归一化函数ef_v(·)是Halphen提出的指数阶乘函数；E为期望；var()为方差；λ₁为拉格朗日乘子；将采样模块中采集河流历年最大流量序列数据即x序列代入以上公式，求出Halphen IB分布的参数α、m和v；

分析模块，用于对洪水概率密度函数进行积分，得到分布函数F(x)：

其中，P表示流量小于X的洪水发生概率；利用分布函数推求T年一遇的洪水流量设计值X，其中

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下技术特征及有益效果：

1、相对于常用的传统分布，Halphen IB分布函数簇本身包含了正态分布、Gamma分布和逆Gamma分布，具有足够的灵活性来模拟复杂多变的数据集，尾部特性十分良好，因此非常适于水文频率分析；

2、其拟合效果基本优于水文中的其它常用分布，最大熵原理能够有效地估计分布函数的参数，得到精度高的洪水频率分布函数。

附图说明

图1为本发明方法流程图；

图2为本发明方法拟合洪水概率密度曲线效果图；

图3为本发明方法拟合洪水概率密度分布曲线效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

(1)对特定水文站进行取样,采集年最大流量样本序列即x序列；

(2)采用Halphen IB分布构建洪水概率密度函数：

Halphen IB型分布的洪水概率密度函数表达式为

其中，变量x表示年最大流量，f(x)表示流量为x的概率密度；α和v是形状参数，m是尺度参数，且m＞0，α＞0；α、m和v共同决定Halphen IB分布的形状；

(3)基于最大熵原理，推导出洪水概率密度函数的参数满足下式：

其中，归一化函数ef_v(·)是Halphen提出的指数阶乘函数；E为期望；var()为方差；λ₁为拉格朗日乘子；将步骤(1)中采集河流历年最大流量序列数据即x序列代入以上公式，求出洪水概率密度函数的参数α、m和v；

(4)对洪水概率密度函数进行积分，得到分布函数F(x)：

其中，P表示流量小于X的洪水发生概率；利用分布函数推求T年一遇的洪水流量设计值X，其中

实施例：利用某流域典型水文站点A的年最大流量序列数据检验Halphen IB分布的拟合效果，将水文站点A的年最大流量序列数据代入参数推导公式中得到Halphen IB分布的参数α＝5.6685,m＝7.3715,v＝0.9448，将参数值代入洪水概率密度函数，绘制函数曲线，如图2所示为洪峰流量和其概率密度的拟合情况，由图2可以看出洪水概率密度函数拟合效果良好；对概率密度函数进行积分得到分布函数，绘制分布函数曲线，如图3所示，图中离散点和曲线分别为洪峰流量的经验值频率值和理论频率值，由图3可以看出分布函数拟合效果良好；

采用Halphen IB分布、正态分布、指数分布(EXP)、Gamma分布、Gumbel分布、广义正态分布(Generalized normal，GN)、P-III分布、广义Pareto分布(Generalized Pareto，GP)及Weibull分布等拟合典型水文站A的年最大洪峰流量序列，并基于K-S(Kolmogorov-Smirnov，K-S)检验法、均方根误差以及AIC准则(Akaike information criterion，AIC)对各分部的拟合结果进行比较分析，确定拟合最优的分布函数，进行洪水频率分析。取K-S检验显著性水平为α＝0.05，P值大于0.05时通过检验，RMSE和AIC值越小说明拟合效果越好。表1给出了各分布K-S检验统计量P值、RMSE和AIC值。结果表明，所有分布均通过了K-S检验。比较各分部的RMSE和AIC值可知，Halphen IB型分布的拟合效果优于常用分布。

表1

表2给出了上述分布函数拟合典型水文站点A的设计洪水值，结果表明，当重现期较大时，各分布计算结果差异明显，考虑到工程设计中，低估的设计值将大大增加大坝及下游的洪水风险，由表中数据比较分析可知Halphen IB型分布的频率分析结果优于常用分布。

表2

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于Halphen IB分布的洪水频率分析方法，其特征在于，该方法包括以下步骤：

(1)对特定水文站进行取样,采集年最大流量样本序列即x序列；

(2)采用Halphen IB分布构建洪水概率密度函数：

$f\left(x\right)=\frac{x^{-2v-1}\exp (-{\left(\frac{m}{x}\right)}^2+\mathrm{\α}(\frac{m}{x}\left)\right)}{m^{-2v}\underset{0}{\overset{\mathrm{\∞}}{\∫}}{x}^{2v-1}\exp (-x^2+\mathrm{\α}x)dt};$

其中，变量x表示年最大流量，f(x)表示流量为x的概率密度；α和v是形状参数，m是尺度参数，且m＞0，α＞0；

(3)洪水概率密度函数的参数α、m和v满足下式：

$-\frac{\underset{0}{\overset{\mathrm{\∞}}{\∫}}\ln x\exp (-{\mathrm{\λ}}_1\ln x-{\left(\frac{m}{x}\right)}^2+\mathrm{\α}(\frac{m}{x}\left)\right)dx}{\underset{0}{\overset{\mathrm{\∞}}{\∫}}\exp (-{\mathrm{\λ}}_1\ln x-{\left(\frac{m}{x}\right)}^2+\mathrm{\α}(\frac{m}{x}\left)\right)dx}=-E\left(\ln x\right)$

$\frac{{\mathrm{ef}}_{v+\frac{1}{2}}\left(\mathrm{\α}\right)}{{\mathrm{ef}}_v\left(\mathrm{\α}\right)}=E\left(\frac{m}{x}\right)$

$\frac{{\mathrm{ef}}_{v+1}\left(\mathrm{\α}\right){\mathrm{ef}}_v\left(\mathrm{\α}\right)-{\left({\mathrm{ef}}_{v+\frac{1}{2}}\right(\mathrm{\α}\left)\right)}^2}{{\left({\mathrm{ef}}_v\right(\mathrm{\α}\left)\right)}^2}=\mathrm{var}\left(\frac{m}{x}\right)$

其中，归一化函数ef_v(·)是Halphen提出的指数阶乘函数；E为期望；var()为方差；λ₁为拉格朗日乘子；将步骤(1)中采集河流历年最大流量序列数据即x序列代入以上公式，求出洪水概率密度函数的参数α、m和v；

(4)对洪水概率密度函数进行积分，得到分布函数F(x)：

$F\left(x\right)={\&Integral;}_0^{x}f\left(x\right)dx=P(X<x)$

其中，P表示流量小于X的洪水发生概率；利用分布函数推求T年一遇的洪水流量设计值X，其中

2.一种基于Halphen IB分布的洪水频率分析系统，其特征在于，该系统分为以下模块：

采样模块，用于对特定水文站进行取样,采集年最大流量样本序列即x序列；

建模模块，用于采用Halphen IB分布构建洪水概率密度模型：

$f\left(x\right)=\frac{x^{-2v-1}\exp (-{\left(\frac{m}{x}\right)}^2+\mathrm{\&alpha;}(\frac{m}{x}\left)\right)}{m^{-2v}\underset{0}{\overset{\mathrm{\&infin;}}{\&Integral;}}{x}^{2v-1}\exp (-x^2+\mathrm{\&alpha;}x)dt}$

其中，变量x表示年最大流量，f(x)表示流量为x的概率密度；α和v是形状参数，m是尺度参数，且m＞0，α＞0；

参数推导模块，用于基于最大熵原理，推导出Halphen IB分布函数的参数满足下式：

$\begin{array}{c}-\frac{\underset{0}{\overset{\mathrm{\&infin;}}{\&Integral;}}\ln x\exp (-{\mathrm{\&lambda;}}_1\ln x-{\left(\frac{m}{x}\right)}^2+\mathrm{\&alpha;}(\frac{m}{x}\left)\right)dx}{\underset{0}{\overset{\mathrm{\&infin;}}{\&Integral;}}\exp (-{\mathrm{\&lambda;}}_1\ln x-{\left(\frac{m}{x}\right)}^2+\mathrm{\&alpha;}(\frac{m}{x}\left)\right)dx}=-E\left(\ln x\right)\\ \frac{{\mathrm{ef}}_{v+\frac{1}{2}}\left(\mathrm{\&alpha;}\right)}{{\mathrm{ef}}_v\left(\mathrm{\&alpha;}\right)}=E\left(\frac{m}{x}\right)\\ \frac{{\mathrm{ef}}_{v+1}\left(\mathrm{\&alpha;}\right){\mathrm{ef}}_v\left(\mathrm{\&alpha;}\right)-{\left({\mathrm{ef}}_{v+\frac{1}{2}}\right(\mathrm{\&alpha;}\left)\right)}^2}{{\left({\mathrm{ef}}_v\right(\mathrm{\&alpha;}\left)\right)}^2}=\mathrm{var}\left(\frac{m}{x}\right)\end{array};$

其中，归一化函数ef_v(·)是Halphen提出的指数阶乘函数；E为期望；var()为方差；λ₁为拉格朗日乘子；将采样模块中采集河流历年最大流量序列数据即x序列代入以上公式，求出Halphen IB分布的参数α、m和v；

分析模块，用于对洪水概率密度函数进行积分，得到分布函数F(x)：

$F\left(x\right)={\&Integral;}_0^{x}f\left(x\right)dx=P(X<x)$

其中，P表示流量小于X的洪水发生概率；利用分布函数推求T年一遇的洪水流量设计值X，其中