袁 超，王 盼，李幼鳳

(中國電建集團西北勘測設(shè)計研究院有限公司， 陜西 西安 710065)

水利水電工程設(shè)計洪水計算常用統(tǒng)計分析方法，主要涉及頻率分布模型的選擇和參數(shù)估計方法的采用兩方面的工作[1]。目前，人們還無法從水文機理上證明水文事件的概率(頻率)分布函數(shù)。實際中，計算者選用現(xiàn)有的隨機變量分布函數(shù)，根據(jù)收集的水文數(shù)據(jù)進行審查，通過選用合理的參數(shù)估計方法進行分布函數(shù)參數(shù)計算，經(jīng)分布函數(shù)擬合度檢驗后，依據(jù)水文數(shù)據(jù)的經(jīng)驗頻率與選用分布函數(shù)理論頻率之間的擬合效果來評估水文序列頻率分布函數(shù)的擬合效果[2]。

洪水頻率分布模型大致可分為正態(tài)分布類、Γ分布類、極值分布類、指數(shù)分布類等[1]，其中廣義極值分布是國際上水文頻率計算采用較多的分布線型[1]。1998年，金光炎[3]教授對廣義極值分布的統(tǒng)計性質(zhì)、離均系數(shù)計算進行了系統(tǒng)分析研究。2006年，劉聰?shù)萚4]基于廣義極值分布研究了最大風速頻率分布模型。2011年，謝欣榮[5]采用均值、變差系數(shù)和偏態(tài)系數(shù)等三個統(tǒng)計參數(shù)計算廣義極值分布曲線，并給出了離均系數(shù)表。金連根[6]，李揚[7]，許弘喆等[8]先后利用極值分布在設(shè)計風速、降水、潮水位推算中的應(yīng)用。

洪水概率分布參數(shù)估計常見方法有矩法、極大似然法、概率權(quán)重矩法、線性矩法、權(quán)函數(shù)法及優(yōu)化適線法等，除矩法外，上述各法均與所選定的分布模型有關(guān)[1]。Hosking(1990)在常規(guī)矩的基礎(chǔ)上提出了線性距(L-moments)的概念，定義次序統(tǒng)計量線性組合的期望值為線性距。線性矩是概率加權(quán)矩的線性組合，是對概率權(quán)重矩的一種改進，與常規(guī)矩法相比，估計量具有良好的不偏性，變差系數(shù)和偏態(tài)系數(shù)的誤差明顯要小些[5]。2003年，陳元芳等[9]提出了兩種可考慮歷史洪水的樣本線性距估計公式。2004年，段忠東等[10]比較了極值概率分布參數(shù)估計方法。2005年，金光炎[11]給出了矩、概率權(quán)重矩與線性矩之間的關(guān)系，分析了計算參數(shù)的近似公式及其精度。2007年，陳民等[12]分析了極值分布統(tǒng)計參數(shù)與線性矩的關(guān)系。2008年，陳元芳等[13]給出了具有歷史洪水系列的線性矩公式，并對參數(shù)k值作了改進。國內(nèi)學者對廣義極值分布參數(shù)估計方法進行了比較[14-18]。

基于上述研究成果，本文對廣義極值分布函數(shù)分布參數(shù)、統(tǒng)計特征參數(shù)、線性矩、離均系數(shù)等相互關(guān)系與計算方法進行了全面總結(jié)，并與國內(nèi)常用皮爾遜Ⅲ型分布在統(tǒng)計參數(shù)、設(shè)計成果等方面差異性進行了對比分析。

1 廣義極值分布

1.1 函數(shù)形式

1928年，F(xiàn)isher和Tippett證明極值具有漸近分布函數(shù)(即其極限概率分布)，并概括了與原始分布對應(yīng)的通常有三種類型的極限概率分布(漸近的極值分布)模型。其中，第Ⅰ型為指數(shù)原始分布，又稱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、α和ξ分別為形狀參數(shù)、尺度參數(shù)和位置參數(shù)。k=0時為第Ⅰ型分布，k<0時為第Ⅱ型分布，k>0時為第Ⅲ型分布。其概率密度函數(shù)為：

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 分布參數(shù)與統(tǒng)計特征參數(shù)的關(guān)系

(1) 當k≠0且k>-1/3時，分布參數(shù)與統(tǒng)計特征參數(shù)關(guān)系式為[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為離差系數(shù)、Cs為偏態(tài)系數(shù)，sign為符號函數(shù)，Γ(*)為gamma函數(shù)。

(2) 當k=0時，分布參數(shù)與統(tǒng)計特征參數(shù)關(guān)系式為[5]：

(10)

(11)

EX=ξ+αγ

(12)

(13)

Cs=1.139 547 1

(14)

式中：γ為歐拉常數(shù)，約等于0.577 215 7。

2 線性矩與參數(shù)估計

2.1 分布函數(shù)的線性矩

假定序列X服從特定分布函數(shù)，從序列中挑選出一組容量為n的樣本，其次序統(tǒng)計量為x1∶n≤x2∶n≤…xn∶n，則變量所服從分布函數(shù)的前四階線性矩可以表示為[18]：

(15)

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

線性矩法中，定義線性離差系數(shù)τ2(即L-Cv)、線性偏態(tài)系數(shù)τ3(即L-Cs)、線性峰度系數(shù)τ4(即L-Ce)分別表示為[11，18]：

(16)

2.2 樣本線性矩的計算公式

當樣本容量有限，變量X的n觀測值為x1≤x2…≤xn，則λ1，λ2，λ3對應(yīng)的樣本矩l1，l2，l3及線性離差系數(shù)τ2、線性偏態(tài)系數(shù)τ3、線性峰度系數(shù)τ4計算公式如下[12-13]：

(17)

當樣本為連續(xù)序列時，b0、b1、b2的計算公式如下[13]：

(18)

(19)

(20)

當樣本為具有歷史洪水的不連續(xù)序列時，b0、b1、b2的計算公式如下[13]：

(21)

(22)

(23)

式中：N為歷史洪水考證期；a為歷史洪水個數(shù)；l為實測洪水系列中特大洪水個數(shù)；n為實測洪水系列長度，n0=n-l+a為洪水系列總個數(shù)。在此指出，文獻[7]中b1、b2中第一項求和公式m初值有誤，應(yīng)分別為2、3。

2.3 線性矩與分布參數(shù)的關(guān)系

(1) 當k≠0且k>-1/3時，線性矩、分布參數(shù)及統(tǒng)計特征參數(shù)關(guān)系式為[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時，分布參數(shù)與統(tǒng)計特征參數(shù)關(guān)系式為[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 設(shè)計值推求

頻率P的設(shè)計值Xp一般式如下：

Xp=EX(1+ΦpCv)

(35)

(1) 當k≠0且k>-1/3時，離均系數(shù)Φp的計算公式為[5]：

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

(36)

(2) 當k=0時，離均系數(shù)Φp的計算公式為[5]：

(37)

至此，利用式(18)—式(19)或式(20)—式(22)與式(17)可求得樣本矩l1，l2，l3及線性偏態(tài)系數(shù)τ3。利用式(33)或式(34)求得k值，代入式(36)或式(37)求得離均系數(shù)Φp，代入式(9)求得偏態(tài)系數(shù)Cs；利用式(24)—式(27)或式(28)—式(27)，求得均值EX、離差系數(shù)Cv。

4 應(yīng)用實例

某水文站有1940年—2019年年最大洪峰流量系列，調(diào)查歷史洪水有1864年、1867年、1931年，1931年洪水量級與1980年實測洪水相當，不作特大值處理，加入連續(xù)系列進行頻率計算。N=156，a=2；l=0，n=81，n0=n-l+a=83。由式(21)—式(23)求得參數(shù)b0、b1、b2，按線性矩法估計統(tǒng)計特征參數(shù)，結(jié)果見表1。為驗證GEV分布的適用性，采用文獻[11]中線性矩法估計P-Ш型頻率曲線統(tǒng)計特征參數(shù)，結(jié)果見表1。GEV分布與P-Ш分布頻率曲線圖見圖1。

表1 GEV分布與P-Ш分布統(tǒng)計特征參數(shù)估計成果表

圖1 GEV分布與P-Ш分布頻率曲線圖(線性矩法)

從參數(shù)估計看(見表1)，線性矩法與矩法估計均值相同，線性矩法估計Cv與Cs值偏大。對于GEV分布與P-Ш分布，線性矩法估計均值相同與Cv值相近，Cs值相差較大。從計算結(jié)果看(見圖1)，線性矩法求得的GEV分布與P-Ш分布頻率曲線交叉，GEV分布頻率曲線上部上翹，中部下凹，尾部可能會出現(xiàn)負值的不合理情況，即當重現(xiàn)期大于200年時，GEV分布計算設(shè)計成果偏大，且隨著重現(xiàn)期的增大，兩種線型偏離程度越大；當重現(xiàn)期在10年～200年時，GEV分布計算成果偏小，但兩者相差不大，因此，采用GEV分布推求設(shè)計洪水具有一定的合理性，但推求校核洪水的合理性仍有待進一步研究。

5 結(jié) 論

本文對廣義極值分布函數(shù)分布參數(shù)、統(tǒng)計特征參數(shù)、線性矩、離均系數(shù)等相互關(guān)系與計算方法進行了全面總結(jié)，并與國內(nèi)常用皮爾遜Ⅲ型分布在統(tǒng)計參數(shù)、設(shè)計成果等方面的差異性進行了對比分析，得出以下幾點結(jié)論：

(1) 線性矩法與矩法估計均值相同，線性矩法估計Cv與Cs值偏大。對于GEV分布與P-Ш分布，線性矩法估計均值相同與Cv值相近，Cs值相差較大。

(2) 與P-Ш分布相比，GEV分布曲線上部上翹，中部下凹，尾部可能會出現(xiàn)負值的不合理情況，即隨著重現(xiàn)期的增大，兩種線型偏離程度越大；在曲線中部，GEV分布計算成果偏小，但兩者相差不大。

(3) 廣義極值分布推求設(shè)計洪水具有一定的適用性，但推求校核洪水的可靠性仍有待進一步研究。