齊哲嫻，宋松柏

(1 西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100；2 浙江大學(xué) 建筑工程學(xué)院，浙江 杭州 310058)

?

廣義極值分布序列經(jīng)驗概率的計算

齊哲嫻1,2，宋松柏1

(1 西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100；2 浙江大學(xué) 建筑工程學(xué)院，浙江 杭州 310058)

【目的】 研究并建立廣義極值分布無偏經(jīng)驗概率計算公式，為廣義極值分布經(jīng)驗概率的計算提供支持?！痉椒ā?應(yīng)用次序統(tǒng)計量原理，推導(dǎo)廣義極值分布無偏經(jīng)驗概率計算公式，采用統(tǒng)計試驗方法將推導(dǎo)的公式與現(xiàn)有的經(jīng)驗頻率公式進(jìn)行對比，對其進(jìn)行檢驗，最后以陜北地區(qū)12個水文測站的年最大洪峰流量系列為例進(jìn)行模型應(yīng)用?！窘Y(jié)果】 推導(dǎo)出了便于工程設(shè)計應(yīng)用的廣義極值分布經(jīng)驗概率計算公式GEVQG。統(tǒng)計試驗和實例應(yīng)用表明：推導(dǎo)的公式GEVQG和現(xiàn)有經(jīng)驗頻率公式Cunnane公式的相對誤差和偏差均較小，并且對研究區(qū)的擬合效果良好?！窘Y(jié)論】 推導(dǎo)的計算公式GEVQG和Cunnane公式均可以作為廣義極值分布經(jīng)驗概率計算的優(yōu)選公式，為工程水文經(jīng)驗概率計算提供了新的選擇。

無偏經(jīng)驗頻率公式；廣義極值分布；洪水頻率分析；次序統(tǒng)計量；陜北地區(qū)

繪點位置公式也稱經(jīng)驗頻率公式，是目前適線法確定分布參數(shù)和評定各類參數(shù)估算方法的依據(jù)。從20世紀(jì)50年代開始，我國一直沿用數(shù)學(xué)期望公式進(jìn)行P-Ⅲ型分布參數(shù)適線法估算。武漢大學(xué)郭生練等[1]指出：不少人在這種不合適的繪點位置基礎(chǔ)上推敲各種適線準(zhǔn)則的優(yōu)劣，或抽取“理想”樣本作為比較各種參數(shù)估計方法的基準(zhǔn)，這顯然勞而無功或不太合理。

Hazen[2]最早將經(jīng)驗頻率公式應(yīng)用到水文頻率分析中。Weibull[3]以次序統(tǒng)計量頻率的數(shù)學(xué)期望值(E[F(xm)])為依據(jù)，建立了數(shù)學(xué)期望公式。Kimball[4]推薦以次序統(tǒng)計量的數(shù)學(xué)期望頻率(F[E(xm)])為依據(jù)建立經(jīng)驗概率公式。Harter[5]指出，經(jīng)驗頻率公式選擇存在各種混亂和爭議的原因在于F[E(xm)]不等于E[F(xm)]，并研究推薦使用中值公式。當(dāng)樣本長度小于20時，可以利用Johnson[6]給出的表格查算經(jīng)驗頻率；而當(dāng)樣本長度大于20時，可以應(yīng)用Blom公式[7]。

Gringorten[8]建立了適用于極值Ⅰ型分布的以F[E(xm)]為依據(jù)的經(jīng)驗頻率公式。Cunnane[9]提出了一個與分布無關(guān)的無偏經(jīng)驗概率公式。Arnell等[10]利用概率權(quán)重矩法推求了廣義極值分布的無偏繪點位置公式，由于數(shù)值計算問題該公式僅適用于樣本容量小于35的情況，但其也給出了樣本容量更大時的近似求解公式。In-na等[11]在Arnell公式的基礎(chǔ)上，進(jìn)一步發(fā)展了廣義極值分布無偏經(jīng)驗頻率公式。Goel等[12]利用概率權(quán)重矩法提出了具有偏態(tài)系數(shù)的廣義極值分布經(jīng)驗頻率公式。De[13]修正了Gringorten公式。Kim等[14]利用遺傳算法推求了廣義極值分布的經(jīng)驗頻率公式。

郭生練[15]選用湖北省42個測站的年最大流量系列，分析對比了6種線型，發(fā)現(xiàn)廣義極值分布和P-Ⅲ分布一樣都能較好地擬合實測資料。陳興旺[16]以南昌市年汛期最大降水量為例，采用廣義極值分布理論進(jìn)行擬合，認(rèn)為該理論分布合理可靠，可以滿足實際工程設(shè)計的需要。上述研究初步說明，廣義極值分布基本適用于我國的水文情勢[17]，但由于水文情勢影響因素的復(fù)雜性，不同地區(qū)適宜的水文頻率和參數(shù)估計方法具有一定的差異，而目前我國尚缺乏不同地區(qū)水文廣義極值分布無偏經(jīng)驗頻率計算理論的研究報道。

本研究采用水文統(tǒng)計原理和數(shù)值計算方法，研究和推導(dǎo)極值Ⅰ型、極值Ⅱ型和極值Ⅲ型分布無偏經(jīng)驗頻率計算公式和求解技術(shù)，在Goel和De[12]推薦的概化公式基礎(chǔ)上，推導(dǎo)理論上更適合廣義極值分布經(jīng)驗概率的計算公式，并與現(xiàn)有的經(jīng)驗頻率公式進(jìn)行對比，利用統(tǒng)計試驗優(yōu)選最佳經(jīng)驗頻率公式，最后以陜北地區(qū)12個水文測站的洪水資料進(jìn)行洪水分布參數(shù)和設(shè)計值推求，以期為研究區(qū)水利工程建筑物的設(shè)計提供科學(xué)依據(jù)，并為工程水文經(jīng)驗概率計算提供理論支撐。

1 理論與方法

1.1 廣義極值分布[12,18]

廣義極值分布(GEV)的定義為：

F(x)=P(X

exp{-[1-k(x-u)/α]1/k},k≠0;

(1)

F(x)=P(X

exp{-exp[-(x-u)/α]},k=0。

(2)

式中：X為隨機(jī)變量，u、α和k分別為位置、尺度和形狀參數(shù)。

廣義極值分布分為Ⅰ、Ⅱ、Ⅲ型，即有：

1)k=0時，分布為極值Ⅰ型分布(EV1)，也稱耿貝爾分布(Gumbel distribution)，x的取值范圍為-∞

2)k<0時，分布為極值Ⅱ型分布(EV2)，也稱弗雷德分布(Frechet distribution)，x的取值范圍為(u+α/k)≤x<∞；分布曲線的下端有限，上端無限。

3)k>0時，分布為極值Ⅲ型分布(EV3)，也稱威布爾分布(Weibull distribution)，x的取值范圍為-∞

形狀參數(shù)k與偏態(tài)系數(shù)Cs的關(guān)系見圖1。

圖 1 廣義極值分布形狀參數(shù)k和偏態(tài)系數(shù)Cs的關(guān)系

形狀參數(shù)k與偏態(tài)系數(shù)Cs的關(guān)系式為：

(3)

μ2=Γ(1+2k)-Γ2(1+k);

(4)

μ3=Γ(1+3k)-3Γ(1+2k)Γ(1+k)+

2Γ3(1+k),k<0;

(5)

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

2Γ3(1+k)],k>0。

(6)

1.2 廣義極值分布次序統(tǒng)計量的數(shù)學(xué)期望[12,19]

給定樣本容量為n的次序化隨機(jī)觀測樣本y(1)≤y(2)≤…≤y(n)，第m階次序統(tǒng)計量的概率密度函數(shù)為：

(7)

式中：F(·)和f(·)分別為樣本的分布函數(shù)和密度函數(shù)。

則y(m)的數(shù)學(xué)期望為：

(8)

由此可得極值Ⅰ型、Ⅱ型、Ⅲ型分布第m階次序統(tǒng)計量的數(shù)學(xué)期望依次為：

(9)

(10)

(11)

將式(9)、(10)和(11)分別應(yīng)用二項級數(shù)展開，進(jìn)行計算后可得：

(-1)s[γ+ln(m+s)](m+s)-1，

(12)

(13)

E(y(3m))=-E(y(2m))。

(14)

式中：γ為歐拉常數(shù)，γ=0.577 2；s為參數(shù)；其他符號意義同前。

1.3 廣義極值分布經(jīng)驗概率的計算

以F(E(Xm))作為繪點位置不僅具有良好的理論基礎(chǔ)支撐，而且會得到良好的統(tǒng)計性質(zhì)[1]。因此，廣義極值分布第m階次序統(tǒng)計量y(m)的繪點位置計算公式為：

Pm=F[E(y(m))]。

(15)

E(y(m))為上述方法推得廣義極值分布次序統(tǒng)計量的數(shù)學(xué)期望。為了便于工程設(shè)計的使用，通過綜合歸納，多數(shù)經(jīng)驗頻率公式可以概化為下列通式：

(16)

nPm-m=a+bPm。

(17)

式中：a和b為未知系數(shù)。

式(17)右邊是Pm的線性函數(shù)，對于給定的樣本容量n和偏態(tài)系數(shù)Cs(或形狀參數(shù)k)，(nPm-m)與Pm之間存在線性關(guān)系，應(yīng)用式(12)、(13)、(14)計算廣義極值分布第m階次序統(tǒng)計量的數(shù)學(xué)期望E(y(m))?？闪顈=E(y(m))，應(yīng)用下式計算Pm，然后利用MATLAB軟件將(nPm-m)和Pm進(jìn)行回歸分析，即可得到a和b的值。

(18)

本研究中形狀參數(shù)k取-0.3 ～1.5，k計算步長為0.1。n取10～100，n計算步長為5。參照公式(3)中k與Cs的關(guān)系，經(jīng)過回歸分析計算后，廣義極值分布經(jīng)驗概率公式有：

(19)

將式(19)命名為GEVQG公式。

2 統(tǒng)計試驗研究

為了比較得出廣義極值分布連續(xù)序列最佳經(jīng)驗概率公式，下面將用2個試驗來檢驗該公式的描述能力和預(yù)見能力。針對廣義極值分布，在現(xiàn)有的經(jīng)驗頻率公式中選擇Hazen公式、中值(Median)公式、Weilbull公式、Gringorten公式(只適用于極值Ⅰ型分布,k=0)、Cunnane公式(替代Gringorten公式用于極值Ⅱ、Ⅲ型分布)、Arnell公式、In-na公式和GEVAC公式，通過統(tǒng)計試驗方法評價GEVQG與以上經(jīng)驗頻率公式優(yōu)劣性。

2.1 試驗1[1,20]

(20)

然后，計算出j個經(jīng)驗頻率公式[2-3,8-12,20]的P(i),j值，根據(jù)下式計算X(i),j：

X(i),j=F-1(P(i),j),i=1,2,…,N。

(21)

式中：X(i),j為采用j公式計算的第i階次序統(tǒng)計量，F(xiàn)-1代表廣義極值分布的逆函數(shù)，i、j分別代表樣本次序和所采用的公式。

(22)

(23)

表1列出了樣本容量分別為30，50，70時不同經(jīng)驗頻率公式在極值Ⅰ型(k=0，均值EX=100，變差系數(shù)Cv=0.5，Cs=1.139)、Ⅱ型(k=-0.2，均值EX=100，變差系數(shù)Cv=0.5，Cs=3.535)和Ⅲ型(k=0.05，均值EX=100，變差系數(shù)Cv=0.5，Cs=0.868)分布時的平均絕對偏差(AD)和均方根誤差(RMSE)的計算結(jié)果。綜合來看，GEVQG公式誤差較小，In-na公式在Cs值較小時誤差較小，GEVAC公式和Cunnane公式的誤差值比較接近，中值(Median)公式和Hazen公式有較大偏差，Weibull公式和Arnell公式的平均絕對偏差和均方根誤差均較大(表1)。

表 1 不同經(jīng)驗頻率公式平均絕對偏差(AD)和均方根誤差(RMSE)的比較

2.2 試驗2[1,21]

(24)

圖2為k=0.1，k=-0.1，k=0時廣義極值分布的重現(xiàn)期與相對誤差的關(guān)系。從圖2可以看出，k=0.1和k=-0.1時，Arnell公式和Weibull公式相對誤差偏大；k=0時，Weibull公式和中值公式相對誤差也較大；In-na公式在部分重現(xiàn)期表現(xiàn)良好，部分重現(xiàn)期相對誤差較大；Cunnane公式、Gringorten公式和GEVAC公式的相對誤差較小，且三者的相對誤差值也很接近；GEVQG公式在分位數(shù)估計中表現(xiàn)最為良好。整體上來看，重現(xiàn)期的增大使得各個經(jīng)驗頻率公式的相對誤差均隨之增大。

3 實例應(yīng)用

選擇陜西省陜北地區(qū)12個水文測站(神木(二)、趙石窯、綏德、交口河、杏河、棗園、吳旗、劉家河、安塞、黃陵、張村驛和志丹)的年最大洪峰流量系列，在資料審查的基礎(chǔ)上，采用兼具廣泛性和良好統(tǒng)計特性的GEVQG公式、GEVAC公式和Cunnane公式進(jìn)行洪水頻率分析。

利用編制的MATLAB程序?qū)ρ芯繀^(qū)12個水文測站的年最大洪峰流量系列進(jìn)行頻率計算，獲得3種經(jīng)驗公式下的點據(jù)分布，并與理論頻率進(jìn)行相關(guān)關(guān)系比較，其中以代表性水文測站黃陵站和劉家河站為例的計算結(jié)果見圖3。從圖3可以看出，經(jīng)驗點數(shù)據(jù)基本均聚集在45°線附近，由此證明理論分布和樣本系列的擬合效果較為良好，即選擇的廣義極值分布理論與該地區(qū)的適合程度較為良好。

圖 2 不同形狀參數(shù)(k)下廣義極值分布的重現(xiàn)期與相對誤差的關(guān)系

△.GEVAC公式 GEVAC formula；*.Cunnane公式 Cunnane formula；○.GEVQG公式 GEVQG formula

圖 3 陜北地區(qū)2個水文測站年最大洪峰流量經(jīng)驗頻率和理論頻率的相關(guān)關(guān)系

Fig.3 Relationship between empirical and theoretical frequencies for annual maximum flood peak flow series at two stations in Northern Shaanxi

采用離差平方和最小準(zhǔn)則法(OLS)和相對離差平方和最小準(zhǔn)則法(WLS)進(jìn)行擬合效果的定量評價分析，表2中列出了不同評價方法和不同重現(xiàn)期(50年一遇、100年一遇和200年一遇)下擬合效果最優(yōu)的經(jīng)驗概率公式。從表2可以看出，擬合效果的評價標(biāo)準(zhǔn)不同，得到的結(jié)論也不盡相同?？傮w上而言，GEVQG公式和Cunnane公式對于大多數(shù)測站的擬合效果較為良好。選擇代表性水文測站黃陵站和劉家河站為例繪制頻率曲線圖，結(jié)果如圖4所示。由圖4可以看出，經(jīng)驗點據(jù)和理論頻率曲線的擬合效果良好，其中曲線中部和尾部的擬合程度最好，可以較好地兼顧曲線首部。

表 2 不同評價方法和重現(xiàn)期下陜北各水文測站擬合效果最優(yōu)時對應(yīng)的經(jīng)驗概率公式

◇.GEVAC公式 GEVAC formula；○.GEVGQ公式 GEVGQ formula；×.Cunnane公式 Cunnane formula

圖 4 不同公式擬合陜北地區(qū)2個水文測站年最大洪峰流量頻率曲線的比較

Fig.4 Comparison of frequency curves for annual maximum flood peak flow series at two stations in Northern Shaanxi with different formulas

4 結(jié) 論

1)Goel和De[12]推薦的GEVAC公式基于對稱分布的通式進(jìn)行推導(dǎo)，本研究應(yīng)用的通式從理論上更適合廣義極值分布。

Arnell公式[10]未考慮樣本數(shù)據(jù)的性質(zhì)，變量之間未形成一個準(zhǔn)確的聯(lián)系，只羅列出了部分Cs、N、m的值，在實際應(yīng)用中十分不便，而且在Arnell所列的表格中會出現(xiàn)參數(shù)α=β的情況，顯然這不符合實際。在統(tǒng)計試驗中，Arnell公式的偏差很大，因此不建議使用。In-na公式[11]在Cs總體較小時展現(xiàn)出較為良好的統(tǒng)計特性，但在Cs總體較大時會產(chǎn)生較大的偏差，該公式推導(dǎo)時假定的是一組降序排列的數(shù)值，所以In-na公式很難應(yīng)用到水文頻率分析中。Weibull公式[3]與總體分布無關(guān)，在我國得到了較為廣泛的應(yīng)用，但從統(tǒng)計試驗可以看出，Weibull公式誤差較大，不宜繼續(xù)使用。中值公式[20]與總體分布無關(guān)，使用簡便，一般情況下實測資料被認(rèn)為是隨機(jī)獨立的，但是中值公式頭尾2項的頻率間隔與樣本中間部分不同，不符合樣本各項為等可能獨立事件的基本假定，而且在統(tǒng)計試驗中，其抽樣誤差也較大。Hazen公式[2]是純經(jīng)驗公式，理論基礎(chǔ)薄弱，其重現(xiàn)期T=n/(1-0.5)=2n，即n年觀測資料中最大項的重現(xiàn)期為2n，也就是說如果要得到100年一遇的結(jié)果，只要50年的資料就可以，顯然計算結(jié)果偏于不安全。Gringortern公式[8]只適用于極值Ⅰ型分布，在統(tǒng)計試驗中可以看出，其無偏性和有效性都較優(yōu)。Cunnane公式[9]和GEVQG公式表現(xiàn)都較好，且統(tǒng)計特性接近，只是GEVQG公式與偏態(tài)系數(shù)有關(guān)，需初定統(tǒng)計參數(shù)。

2)對研究區(qū)陜北神木(二)、趙石窯、綏德、交口河、杏河、棗園、吳旗、劉家河、安塞、黃陵、張村驛和志丹等12個水文測站洪水頻率的實例分析研究表明，采用廣義極值分布作為該地區(qū)的頻率分布曲線合理可行。

3)在實際應(yīng)用中，推薦將GEVQG公式和Cunnane公式作為廣義極值分布的最佳經(jīng)驗概率公式。

[1]郭生練,葉守澤.論水文計算中的經(jīng)驗頻率公式 [J].武漢水利電力學(xué)院學(xué)報,1992(2):38-45.

Guo S L,Ye S Z.Another look at plotting position formulae in hydrology [J].Journal of Wuhan University of Hydraulic and Electric Engineering,1992(2):38-45.

[2]Hazen A.Storage to be provided in impounding reservoirs for municipal water supply [J].Transactions of the American Society of Civil Engineers,1914,77(1):15-39.

[3]Weibull W.A statistical theory of the strength of materials [M].Stockholm:Generalstabens Litografiska Anstals Forlag,1939.

[4]Kimball B F.Assignment of frequencies to a completely ordered set of sample data [J].Transactions American Geophysical Union,1946,27(6):843-846.

[5]Harter H L.Another look at plotting positions [J].Communications in Statistics-Theory and Methods,1984,13(13):1613-1633.

[6]Johnson L G.The median ranks of sample values in their population with an application to certain fatigue studies [J].Bulletin of the American Mathematical Society,1950,56(3):247.

[7]Blom G.Statistical estimates and transformed beta-variables [M].New York:Wiley,1958.

[8]Gringorten I I.A plotting rule for extreme probability paper [J].Journal of Geophysical Research,1963,68(3):813-814.

[9]Cunnane C.Unbiased plotting positions：a review [J].Journal of Hydrology,1978,37(3):205-222.

[10]Arnell N W,Beran M,Hosking J R M.Unbiased plotting positions for the general extreme value distribution [J].Journal of Hydrology,1986,86:59-69.

[11]In-na N,Nguyen V.An unbiased plotting position formula for the general extreme value distribution [J].Journal of Hydrology,1989,106(89):193-209.

[12]Goel N K,De M.Development of unbiased plotting position formula for general extreme value distributions [J].Stochastic Hydrology & Hydraulics,1993,7(1):1-13.

[13]De M.A new unbiased plotting position formula for Gumbel distribution [J].Stochastic Environmental Research Risk Assessment,2000,14(1):1-7.

[14]Kim S,Shin H,Joo K,et al.Development of plotting position for the general extreme value distribution [J].Journal of Hydrology,2012,475(12):259-269.

[15]郭生練.考慮歷史洪水資料的頻率計算方法 [J].武漢水利電力學(xué)院學(xué)報,1988(1):15-21.

Guo S L.Flood frequency analysis with historic information [J].Journal of Wuhan University of Hydraulic and Electric Engineering,1988(1):15-21.

[16]陳興旺.廣義極值分布理論在重現(xiàn)期計算的應(yīng)用 [J].氣象與減災(zāi)研究,2008,31(4):52-54.

Chen X W.Application of generalized extreme value distribution theory in calculating return periods [J].Journal of Meteorology and Disaster Reduction Research,2008,31(4):52-54.

[17]原秀紅.基于部分概率權(quán)重矩的洪水頻率分布參數(shù)估計方法研究 [D].陜西楊凌：西北農(nóng)林科技大學(xué),2014.

Yuan X H.Partial probability weighted moments for estimating distribution parameters in flood frequency analysis [D].Yangling，Shaanxi：Northwest A&F University,2014.

[18]Jenkinson A F.The frequency distribution of the annual maximum (or minimum) values of meteorological elements [J].Quarterly Journal of the Royal Meteorological Society,1955,81(348):158-171.

[19]陳元芳，黃振平.水文統(tǒng)計學(xué) [M].北京：中國水利水電出版社，2011.

Chen Y F,Huang Z P.Hydrologic statistics [M].Beijing:China Waterpower Press,2011.

[20]Benard A,Bos-Levenbach E C.The plotting of observations on probability paper [J].Statistica Rijswijk,1953:163-173.

[21]Guo S L.A discussion on unbiased plotting positions for the general extreme value distribution [J].Journal of Hydrology,1990,121(90):33-44.

Plotting position formula for general extreme value distribution

QI Zhexian1,2,SONG Songbai1

(1 College of Water Resources and Architectural Engineering,Northwest A&F University,Yangling,Shaanxi 712100,China;2CollegeofCivilEngineeringandArchitecture,ZhejiangUniversity,Hangzhou,Zhejiang310058,China)

【Objective】 This study developed the unbiased plotting position formula for general extreme value (GEV) distribution.【Method】 Based on the theory of order statistics,this paper developed an unbiased plotting position formula for the GEV distribution.The statistical testing was used to compare the developed formula and other existing formulas.At last,the parameters of annual maximum floods peak flow series of twelve stations in Northern Shaanxi were estimated.【Result】 This paper developed a practical and convenient plotting position formula for the GEV distribution.The case study and statistical testing showed that the developed formula and Cunnane formula had small error and bias as well as good fitness in the study area.【Conclusion】 The developed formula and Cunnane formula can be used as better methods for the GEV distribution,which provides new choice for empirical probability calculation in hydrology engineering.

unbiased plotting position formula;general extreme value distribution;flood frequency analysis;order statistics;Northern Shaanxi

時間:2016-10-20 16:37

10.13207/j.cnki.jnwafu.2016.12.030

2015-07-10

國家自然科學(xué)基金項目(51479171，51179160，71273081)；高等學(xué)校博士點基金項目(20110204110017)

齊哲嫻(1993-)，女，河北唐山人，在讀博士，主要從事市政工程研究。E-mail:qzx@zju.edu.cn

宋松柏(1965-)，男，陜西永壽人，教授，博士，博士生導(dǎo)師，主要從事水文與水資源研究。 E-mail:ssb6533@nwsuaf.edu.cn

P333.9

A

1671-9387(2016)12-0219-07

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/61.1390.S.20161020.1637.060.html