广义极值分布函数参数估计与设计洪水计算方法分析

袁 超，王 盼，李幼凤

(中国电建集团西北勘测设计研究院有限公司， 陕西 西安 710065)

水利水电工程设计洪水计算常用统计分析方法，主要涉及频率分布模型的选择和参数估计方法的采用两方面的工作[1]。目前，人们还无法从水文机理上证明水文事件的概率(频率)分布函数。实际中，计算者选用现有的随机变量分布函数，根据收集的水文数据进行审查，通过选用合理的参数估计方法进行分布函数参数计算，经分布函数拟合度检验后，依据水文数据的经验频率与选用分布函数理论频率之间的拟合效果来评估水文序列频率分布函数的拟合效果[2]。

洪水频率分布模型大致可分为正态分布类、Γ分布类、极值分布类、指数分布类等[1]，其中广义极值分布是国际上水文频率计算采用较多的分布线型[1]。1998年，金光炎[3]教授对广义极值分布的统计性质、离均系数计算进行了系统分析研究。2006年，刘聪等[4]基于广义极值分布研究了最大风速频率分布模型。2011年，谢欣荣[5]采用均值、变差系数和偏态系数等三个统计参数计算广义极值分布曲线，并给出了离均系数表。金连根[6]，李扬[7]，许弘喆等[8]先后利用极值分布在设计风速、降水、潮水位推算中的应用。

洪水概率分布参数估计常见方法有矩法、极大似然法、概率权重矩法、线性矩法、权函数法及优化适线法等，除矩法外，上述各法均与所选定的分布模型有关[1]。Hosking(1990)在常规矩的基础上提出了线性距(L-moments)的概念，定义次序统计量线性组合的期望值为线性距。线性矩是概率加权矩的线性组合，是对概率权重矩的一种改进，与常规矩法相比，估计量具有良好的不偏性，变差系数和偏态系数的误差明显要小些[5]。2003年，陈元芳等[9]提出了两种可考虑历史洪水的样本线性距估计公式。2004年，段忠东等[10]比较了极值概率分布参数估计方法。2005年，金光炎[11]给出了矩、概率权重矩与线性矩之间的关系，分析了计算参数的近似公式及其精度。2007年，陈民等[12]分析了极值分布统计参数与线性矩的关系。2008年，陈元芳等[13]给出了具有历史洪水系列的线性矩公式，并对参数k值作了改进。国内学者对广义极值分布参数估计方法进行了比较[14-18]。

基于上述研究成果，本文对广义极值分布函数分布参数、统计特征参数、线性矩、离均系数等相互关系与计算方法进行了全面总结，并与国内常用皮尔逊Ⅲ型分布在统计参数、设计成果等方面差异性进行了对比分析。

1 广义极值分布

1.1 函数形式

1928年，Fisher和Tippett证明极值具有渐近分布函数(即其极限概率分布)，并概括了与原始分布对应的通常有三种类型的极限概率分布(渐近的极值分布)模型。其中，第Ⅰ型为指数原始分布，又称Gumbel分布(耿贝尔分布)；第Ⅱ型为柯西型原始分布，亦称Fréchet分布(弗雷歇分布)；第Ⅲ型为有界型分布，即Weibull分布(威布尔分布)[5]。

F(x)=exp[-(1-ky)1/k]k≠0

(1)

F(x)=exp[-exp(-y)]k=0

(2)

式中：y=(x-ξ)/α，k、α和ξ分别为形状参数、尺度参数和位置参数。k=0时为第Ⅰ型分布，k<0时为第Ⅱ型分布，k>0时为第Ⅲ型分布。其概率密度函数为：

f(x)=(1-ky)1/k-1exp[-(1-ky)1/k]/αk≠0

(3)

f(x)=exp[-y-exp(-y)]/αk=0

(4)

1.2 分布参数与统计特征参数的关系

(1) 当k≠0且k>-1/3时，分布参数与统计特征参数关系式为[5，12]：

α=|k|EXCv/[Γ(1+2k)-Γ2(1+k)]1/2

(5)

ξ=EX{1-sign(k)*[1-Γ(1+k)]Cv}/[Γ(1+2k)-Γ2(1+k)]1/2

(6)

EX=ξ+α[1-Γ(1+k)]/k

(7)

Cv=α[Γ(1+2k)-Γ2(1+k)]1/2/(|k|EX)

(8)

Cs=[Γ(1+3k)-3Γ(1+2k)Γ(1+k)+

2Γ3(1+k)]/[Γ(1+2k)-Γ2(1+k)]3/2

(9)

式中：EX为均值、Cv为离差系数、Cs为偏态系数，sign为符号函数，Γ(*)为gamma函数。

(2) 当k=0时，分布参数与统计特征参数关系式为[5]：

(10)

(11)

EX=ξ+αγ

(12)

(13)

Cs=1.139 547 1

(14)

式中：γ为欧拉常数，约等于0.577 215 7。

2 线性矩与参数估计

2.1 分布函数的线性矩

假定序列X服从特定分布函数，从序列中挑选出一组容量为n的样本，其次序统计量为x1∶n≤x2∶n≤…xn∶n，则变量所服从分布函数的前四阶线性矩可以表示为[18]：

(15)

式中：E(·)表示为期望，λ1为分布的线性位置(L-location)或为分布函数的均值；λ2为线性尺度(L-scale)，λ3和λ4为衡量分布函数的偏度和峰度的三阶、四阶线性矩。

线性矩法中，定义线性离差系数τ2(即L-Cv)、线性偏态系数τ3(即L-Cs)、线性峰度系数τ4(即L-Ce)分别表示为[11，18]：

(16)

2.2 样本线性矩的计算公式

当样本容量有限，变量X的n观测值为x1≤x2…≤xn，则λ1，λ2，λ3对应的样本矩l1，l2，l3及线性离差系数τ2、线性偏态系数τ3、线性峰度系数τ4计算公式如下[12-13]：

(17)

当样本为连续序列时，b0、b1、b2的计算公式如下[13]：

(18)

(19)

(20)

当样本为具有历史洪水的不连续序列时，b0、b1、b2的计算公式如下[13]：

(21)

(22)

(23)

式中：N为历史洪水考证期；a为历史洪水个数；l为实测洪水系列中特大洪水个数；n为实测洪水系列长度，n0=n-l+a为洪水系列总个数。在此指出，文献[7]中b1、b2中第一项求和公式m初值有误，应分别为2、3。

2.3 线性矩与分布参数的关系

(1) 当k≠0且k>-1/3时，线性矩、分布参数及统计特征参数关系式为[12]：

α=λ2k/[(1-2-k)Γ(1+k)]

(24)

ξ=λ1-α[1-Γ(1+k)]/k

(25)

EX=λ1

(26)

Cv=α[Γ(1+2k)-Γ2(1+k)]1/2/(|k|λ1)

(27)

(2) 当k=0时，分布参数与统计特征参数关系式为[12]：

α=ln2/λ2

(28)

(29)

EX=λ1

(30)

(31)

(3)k值改进公式。由于k=f(Cs)无显式表达式，常用试算或近似公式计算k值。当-0.5≤τ3≤0.5时，k值近似公式为[12-13]：

k≈7.8590c+2.9554c2

(32)

式中：c=2/(3+τ3)-ln2/ln3，上式可将计算误差控制在9×10-4以下[13]。

(1) 当k<0且k>-1/3时，k值改进公式如下[12-13]：

(33)

式中：A0=0.281 948 10，A1=-1.771 128 69，A2=0.694 409 03，A3=-0.219 400 01。

(2) 当k<0时，k值改进公式如下[12-13]：

(34)

式中：B0=0.283 820 48，B1=-1.796 269 58，B2=0.810 656 21，B3=-0.545 344 41，B4=0.738 599 57，B5=1.124 119 25，B6=1.969 619 41。

改进后的k值计算精度有较明显的提高[13]，因此，本文采用改进后的k值计算公式。

3 设计值推求

频率P的设计值Xp一般式如下：

Xp=EX(1+ΦpCv)

(35)

(1) 当k≠0且k>-1/3时，离均系数Φp的计算公式为[5]：

Φp=sign(k)×{Γ(1+k)-[-ln(1-P)]k}/[Γ(1+2k)-Γ2(1+k)]1/2

(36)

(2) 当k=0时，离均系数Φp的计算公式为[5]：

(37)

至此，利用式(18)—式(19)或式(20)—式(22)与式(17)可求得样本矩l1，l2，l3及线性偏态系数τ3。利用式(33)或式(34)求得k值，代入式(36)或式(37)求得离均系数Φp，代入式(9)求得偏态系数Cs；利用式(24)—式(27)或式(28)—式(27)，求得均值EX、离差系数Cv。

4 应用实例

某水文站有1940年—2019年年最大洪峰流量系列，调查历史洪水有1864年、1867年、1931年，1931年洪水量级与1980年实测洪水相当，不作特大值处理，加入连续系列进行频率计算。N=156，a=2；l=0，n=81，n0=n-l+a=83。由式(21)—式(23)求得参数b0、b1、b2，按线性矩法估计统计特征参数，结果见表1。为验证GEV分布的适用性，采用文献[11]中线性矩法估计P-Ш型频率曲线统计特征参数，结果见表1。GEV分布与P-Ш分布频率曲线图见图1。

表1 GEV分布与P-Ш分布统计特征参数估计成果表

图1 GEV分布与P-Ш分布频率曲线图(线性矩法)

从参数估计看(见表1)，线性矩法与矩法估计均值相同，线性矩法估计Cv与Cs值偏大。对于GEV分布与P-Ш分布，线性矩法估计均值相同与Cv值相近，Cs值相差较大。从计算结果看(见图1)，线性矩法求得的GEV分布与P-Ш分布频率曲线交叉，GEV分布频率曲线上部上翘，中部下凹，尾部可能会出现负值的不合理情况，即当重现期大于200年时，GEV分布计算设计成果偏大，且随着重现期的增大，两种线型偏离程度越大；当重现期在10年～200年时，GEV分布计算成果偏小，但两者相差不大，因此，采用GEV分布推求设计洪水具有一定的合理性，但推求校核洪水的合理性仍有待进一步研究。

5 结 论

本文对广义极值分布函数分布参数、统计特征参数、线性矩、离均系数等相互关系与计算方法进行了全面总结，并与国内常用皮尔逊Ⅲ型分布在统计参数、设计成果等方面的差异性进行了对比分析，得出以下几点结论：

(1) 线性矩法与矩法估计均值相同，线性矩法估计Cv与Cs值偏大。对于GEV分布与P-Ш分布，线性矩法估计均值相同与Cv值相近，Cs值相差较大。

(2) 与P-Ш分布相比，GEV分布曲线上部上翘，中部下凹，尾部可能会出现负值的不合理情况，即随着重现期的增大，两种线型偏离程度越大；在曲线中部，GEV分布计算成果偏小，但两者相差不大。

(3) 广义极值分布推求设计洪水具有一定的适用性，但推求校核洪水的可靠性仍有待进一步研究。