李 揚(yáng)

(太原理工大學(xué)水利科學(xué)與工程學(xué)院，太原 030024)

0 引 言

設(shè)計(jì)洪水是估計(jì)水工建筑物抗御洪水能力和各種工程與非工程措施防洪效益的首要依據(jù)，實(shí)際工程設(shè)計(jì)中常需要估計(jì)重現(xiàn)期較長的洪水設(shè)計(jì)值。在一個年最大流量序列中，較小流量和較大流量往往不能同時擬合在某種概率分布曲線附近。為了照顧多數(shù)點(diǎn)據(jù)，擬合時常采用目估適線等方法，側(cè)重?cái)M合頻率曲線的中下部，難以兼顧頻率曲線的上部，往往導(dǎo)致重現(xiàn)期較長的大洪水估計(jì)精度不足。因此，國外研究者提出了針對重現(xiàn)期較長的大洪水段進(jìn)行擬合的幾種參數(shù)估計(jì)方法，基于Cunnane[1]首次提出的刪失樣本概念，Wang[2-6]推導(dǎo)了部分概率權(quán)重矩法(Partial Probability Weighted Moments, PPWM)、高階概率權(quán)重矩法(higher-order probability weighted moments, HPWM)和高階線性矩法(high-order L-moments, LH)的無偏估計(jì)量公式，其中，PPWM法與HPWM法均是基于概率權(quán)重矩法的改進(jìn)參數(shù)估計(jì)方法。

近年來，國內(nèi)研究者開始引入PPWM法、HPWM法及LH矩法進(jìn)行洪水頻率分析的研究。王怡璇等[7]研究了低刪失樣本下廣義極值分布的部分概率權(quán)重矩，并應(yīng)用于陜北地區(qū)主要水文測站的實(shí)測洪水序列擬合，取得了優(yōu)于普通概率權(quán)重矩法的效果；原秀紅等[8]對比了部分概率權(quán)重矩法與矩法在廣義極值分布參數(shù)估計(jì)中的統(tǒng)計(jì)特性，結(jié)果表明基于部分概率權(quán)重矩法的估計(jì)量統(tǒng)計(jì)特性優(yōu)于矩法。Deng[9,10]基于最小交互熵原理將PPWM法應(yīng)用于刪失樣本極值分位數(shù)估計(jì)亦取得了理想的效果。由于部分概率權(quán)重矩的公式推導(dǎo)過程較為復(fù)雜，而高階概率權(quán)重矩的推導(dǎo)和計(jì)算相對簡便，筆者研究了廣義極值分布的高階概率權(quán)重矩[11]，而后推導(dǎo)了P-III分布的高階概率權(quán)重矩及參數(shù)估計(jì)量公式[12]，應(yīng)用于陜北地區(qū)年最大洪峰流量序列的擬合，結(jié)果表明高階概率權(quán)重矩的統(tǒng)計(jì)特性良好，并且可以有效改善大洪水段的擬合效果；肖玲等[13]、胡素端等[14]應(yīng)用高階概率權(quán)重矩進(jìn)行洪水頻率分析的結(jié)果亦顯示：高階概率權(quán)重矩可較好地?cái)M合序列中大洪水段的經(jīng)驗(yàn)點(diǎn)據(jù)。周長讓等[15]推求了廣義Pareto分布的高階概率權(quán)重矩并給出了估參方法，結(jié)果表明高階概率權(quán)重矩法對超定量洪水系列的稀遇頻率洪水段擬合效果較矩法好。此外，王俊珍等[16]和楊惠等[17]分別研究了高階線性矩法在洪水序列和降水序列頻率分析中的應(yīng)用，結(jié)果表明高階線性矩法對序列高尾部分?jǐn)M合效果較好，可以降低設(shè)計(jì)值的估計(jì)偏差。

現(xiàn)有研究表明，這些新型參數(shù)估計(jì)方法在洪水頻率分析中均取得了較為理想的效果，且優(yōu)于矩法、概率權(quán)重矩法等常規(guī)估參方法，計(jì)算雖較常規(guī)估參方法復(fù)雜，但可以借助計(jì)算機(jī)程序簡化計(jì)算過程。然而，現(xiàn)有研究多針對其中一種方法與常規(guī)估參方法進(jìn)行對比，而對不同新型估參方法的橫向比較則不多見。換言之，對于某實(shí)測洪水序列，究竟哪種新型估參方法的統(tǒng)計(jì)特性更穩(wěn)定，估計(jì)精度更高，對序列較大值的擬合效果更為理想？同時，計(jì)算的復(fù)雜程度也是實(shí)際應(yīng)用中選擇估參方法時需要考量的因素。

因此，本文擬對同樣基于概率權(quán)重矩法改進(jìn)的PPWM法和HPWM法進(jìn)行統(tǒng)計(jì)特性和實(shí)際應(yīng)用的比較研究。采用蒙特卡洛試驗(yàn)，對這兩種參數(shù)估計(jì)方法的統(tǒng)計(jì)特性進(jìn)行研究和比較，并通過實(shí)例對其適用性及擬合效果進(jìn)行檢驗(yàn)，以期為實(shí)際應(yīng)用中優(yōu)選參數(shù)估計(jì)方法提供依據(jù)。

1 部分概率權(quán)重矩和高階概率權(quán)重矩

設(shè)隨機(jī)變量X的分布函數(shù)為F(x)=P(X≤x)，其概率權(quán)重矩(Probability Weighted Moments，PWM)的表達(dá)式為[18]：

(1)

式中：r=0,1,2,…表示概率權(quán)重矩的階數(shù)。若給定一個長度為n、服從F分布的樣本，且x(1)≤x(2)≤…≤x(n)，其概率權(quán)重矩βr的無偏估計(jì)br可由下式計(jì)算：

(2)

廣義極值分布(Generalized Extreme Value Distribution, GEV)是國內(nèi)外洪水頻率分析中廣泛應(yīng)用的一種概率分布[19]，其分布函數(shù)為：

(3)

式中：k為形狀參數(shù)，k∈R；α為尺度參數(shù)，α>0；ξ為位置參數(shù)，ξ∈R。

當(dāng)k≠0時，GEV分布的PWM為：

(4)

其無偏估計(jì)量br為：

(5)

1.1 部分概率權(quán)重矩(PPWM)

為有針對性地研究洪水序列較大值的擬合效果，在已知水文序列中選擇大于某一門限值的觀測值組成新的樣本序列，該樣本序列稱為刪失樣本[20]，該門限值稱為刪失門限值。PPWM法的基本思想是給定不同的刪失水平F0和與之對應(yīng)的刪失門限值x0，將小于x0的實(shí)測流量值從序列中去除，剩余較大實(shí)測流量值參與頻率計(jì)算，刪失水平F0=F(x0)，相當(dāng)于序列中樣本點(diǎn)的去除比例。因此，刪失樣本的部分概率權(quán)重矩可以在普通概率權(quán)重矩法的基礎(chǔ)上進(jìn)行擴(kuò)充。將式(4)中F的積分區(qū)間由[0,1]變?yōu)閇F0,1]，即可得到PPWM的表達(dá)式。當(dāng)k≠0時，GEV分布的PPWM為

(6)

式中：P(-(r+1)lnF0，1+k)為不完全gamma函數(shù)，其表達(dá)式為：

(7)

可推出：

(8)

令：

(9)

給定不同k值，計(jì)算擬合z-k曲線k=a0+a1z+a2z2并求得曲線系數(shù)如表1。

可推求出各參數(shù)的估計(jì)量表達(dá)式為：

(10)

表1 不同刪失水平下z-k曲線系數(shù)Tab. 1 Coefficients of z-k curves under different censoring levels

(11)

(12)

(13)

1.2 高階概率權(quán)重矩(HPWM)

PPWM法中，通過定義刪失樣本，小于刪失門限值的樣本點(diǎn)被舍棄，從而不再影響分布的擬合效果。實(shí)際上，較小的流量值對擬合結(jié)果的影響十分有限，擬合效果的好壞仍主要取決于序列的較大值部分，因此，從概率權(quán)重矩的定義出發(fā)，通過提高PWM的階數(shù)r為隨機(jī)變量x賦予更大的權(quán)重，可以更加容易地實(shí)現(xiàn)側(cè)重?cái)M合序列較大值的目的，這也是HPWM法的基本思想。對k≠0，令r=η,η+1,η+2，可推出：

(14)

令：

(15)

給定不同k值，計(jì)算擬合z-k曲線，求得曲線k=a0+a1z+a2z2系數(shù)見表2。

表2 不同階數(shù)下z-k曲線系數(shù)Tab. 2 Coefficients of z-k curves under different orders

可推求出各參數(shù)的估計(jì)量表達(dá)式為：

(16)

(17)

(18)

(19)

2 蒙特卡洛試驗(yàn)

2.1 試驗(yàn)方案

設(shè)置樣本容量為100，模擬次數(shù)為N=10 000，給定位置參數(shù)的初值為ξ=0，尺寸參數(shù)的初值為α=1.0，形狀參數(shù)為k=-0.5和k=0.5，部分概率權(quán)重矩的刪失水平F0=0,0.1,0.2,0.3,0.4,0.5，高階概率權(quán)重矩的階數(shù)η=0，1，2，3，4。

2.2 評價標(biāo)準(zhǔn)

采用偏差(Bias)、標(biāo)準(zhǔn)誤差(Standard Error,SE)和均方根誤差(Root Mean Square Error,RMSE)3個指標(biāo)，評價分布參數(shù)及設(shè)計(jì)值的估計(jì)值的無偏性和有效性。設(shè)計(jì)值的重現(xiàn)期選擇了50年一遇、100年一遇、200年一遇和500年一遇。

偏差Bias反映了估計(jì)量與其數(shù)學(xué)期望的差距，Bias的絕對值越小，表示估計(jì)量的偏差越小，其無偏性也越好，其計(jì)算公式為：

(20)

標(biāo)準(zhǔn)誤差SE反映了估計(jì)量的有效性，對同一總體參數(shù)的估計(jì)量中，SE更小的估計(jì)量更有效，即SE的值越小，估計(jì)量的有效性越好。若統(tǒng)計(jì)試驗(yàn)次數(shù)為N，則SE的計(jì)算公式為：

(21)

均方根誤差RMSE是估計(jì)值與實(shí)際值偏差平方和與試驗(yàn)次數(shù)比值的平方根，用來衡量估計(jì)值與實(shí)際值之間的偏差，RMSE的值越小，表明估計(jì)的偏差越小，有效性越好。若統(tǒng)計(jì)試驗(yàn)次數(shù)為N，則RMSE的計(jì)算公式為：

(22)

2.3 試驗(yàn)結(jié)果

按照以上試驗(yàn)方案，用Matlab R2018a編程進(jìn)行蒙特卡洛試驗(yàn)，結(jié)果如表3～表4及圖1～圖6所示。

可以發(fā)現(xiàn)如下規(guī)律：

(1)對于GEV分布的PPWM法(圖1～圖3)：①參數(shù)和各重現(xiàn)期設(shè)計(jì)值的估計(jì)偏差隨刪失水平F0的提高而增大，參數(shù)的估計(jì)偏差增幅較大，而設(shè)計(jì)值的估計(jì)偏差增幅相對較??；②各重現(xiàn)期設(shè)計(jì)值的估計(jì)偏差不隨重現(xiàn)期的長度而改變；③無論k=-0.5還是k=0.5，參數(shù)估計(jì)的偏差均隨刪失水平的提高而顯著增大，而設(shè)計(jì)值的估計(jì)偏差較小，k=-0.5時設(shè)計(jì)值估計(jì)的無偏性較好，而k=0.5時設(shè)計(jì)值估計(jì)的有效性較好；

(2)對于GEV分布的HPWM法(圖4～圖6)：①k=0.5時，隨著概率權(quán)重矩階數(shù)的升高，參數(shù)及各重現(xiàn)期設(shè)計(jì)值的估計(jì)值均接近無偏估計(jì)，顯著優(yōu)于k=-0.5時的情形；②k=0.5時，參數(shù)估計(jì)的無偏性和有效性均顯著優(yōu)于k=-0.5時的情形；

(3)兩種估計(jì)方法對比：當(dāng)k=-0.5時，GEV分布的HPWM法無論在參數(shù)估計(jì)的無偏性還是有效性方面都顯著優(yōu)于PPWM法，設(shè)計(jì)值的無偏性和有效性差距不顯著，但HPWM法仍優(yōu)于PPWM法。

表3 PPWM法蒙特卡洛試驗(yàn)結(jié)果Tab.3 Monte Carlo test results of PPWM method

表4 HPWM法蒙特卡洛試驗(yàn)結(jié)果Tab.4 Monte Carlo test results of HPWM method

圖1 不同刪失水平下的PPWM法估計(jì)值BiasFig.1 Bias of PPWM estimators in different censoring levels

圖2 不同刪失水平下的PPWM法估計(jì)值SEFig.2 SE of PPWM estimators in different censoring levels

圖3 不同刪失水平下的PPWM法估計(jì)值RMSEFig.3 RMSE of PPWM estimators in different censoring levels

圖4 不同階數(shù)下的HPWM法估計(jì)值BiasFig.4 Bias of HPWM estimators in different orders

圖5 不同階數(shù)下的HPWM法估計(jì)值SEFig.5 SE of HPWM estimators in different orders

圖6 不同階數(shù)下的HPWM法估計(jì)值RMSEFig.6 RMSE of HPWM estimators in different orders

3 實(shí)例應(yīng)用

3.1 資料及研究區(qū)概況

選用陜西省神木水文站1956-2003年的實(shí)測年最大洪峰流量系列資料，系列長度為48年。神木水文站位于陜西省神木縣南郊，系黃河一級支流窟野河中游干流控制站，屬于國家重要水文站，多年平均徑流量4.888 億m3，實(shí)測最大流量高達(dá)13 800 m3/s。窟野河流域地處黃河中游暴雨中心區(qū)，且水土流失嚴(yán)重，流域內(nèi)洪水峰高量小、洪峰尖瘦、暴漲暴落，洪枯流量的變幅也很大，洪災(zāi)頻發(fā)，對當(dāng)?shù)厣a(chǎn)生活造成較大威脅，因此，在該站的洪水頻率計(jì)算中，提高大洪水部分的擬合精度對當(dāng)?shù)氐姆罏?zāi)減災(zāi)工作具有重要意義。

3.2 分布參數(shù)估計(jì)與頻率曲線繪制

分別使用PPWM法及HPWM法，借助Matlab R2018a編程對GEV分布的參數(shù)進(jìn)行估計(jì)，結(jié)果見表5。分別繪制PPWM法和HPWM法的頻率曲線，如圖7和圖8所示，其中F0=0及r=0,1,2相當(dāng)于普通PWM法。

3.3 擬合效果評價

由圖7可知，采用PPWM法估計(jì)參數(shù)并擬合頻率曲線時，隨著刪失水平F0的提高，實(shí)測序列較大值的擬合效果有明顯改善，其中F0=0.4時，擬合效果最佳，F(xiàn)0=0.5時，擬合效果有變差的趨勢。

由圖8可知，采用HPWM法估計(jì)參數(shù)并擬合頻率曲線時，隨著概率權(quán)重矩階數(shù)的提高，在序列較大值部分，理論頻率曲線逐漸接近實(shí)測序列，在階數(shù)為r=4,5,6時，擬合效果最優(yōu)。

表5 GEV分布參數(shù)估計(jì)結(jié)果Tab.5 Parameter estimates of GEV distribution

圖7 神木站PPWM法頻率曲線擬合結(jié)果Fig.7 Curve fittings for different censoring levels of PWMs

圖8 神木站HPWM法頻率曲線擬合結(jié)果Fig.8 Curve fittings for different orders of PWMs

(23)

表6 神木站不同估參方法的設(shè)計(jì)值誤差Tab. 6 Quantile errors using different methods in Shenmu station

為直觀比較兩種方法對同一序列較大值部分?jǐn)M合效果的優(yōu)劣，選擇每種方法中擬合效果最優(yōu)的參數(shù)組合，即F0=0.4和r=4,5,6的情況，單獨(dú)繪制頻率曲線(圖9)進(jìn)行對比。由圖9看出，兩種方法的曲線在P>50%時幾乎重合，擬合效果不相上下，但HPWM法的曲線更接近最大實(shí)測值。

圖9 神木站HPWM法和PPWM法擬合結(jié)果對比Fig.9 Comparison of Curve fittings of HPWMs and PPWMs

綜合蒙特卡洛試驗(yàn)和實(shí)例應(yīng)用的結(jié)果可知，HPWM方法的統(tǒng)計(jì)特性更加優(yōu)良且穩(wěn)定，實(shí)例應(yīng)用中HPWM法和PPWM法的擬合效果基本相當(dāng)，PPWM方法在P=60%～98%時的設(shè)計(jì)值的累積誤差更小，可根據(jù)實(shí)際情況靈活選用。

4 結(jié) 論

現(xiàn)有研究表明，HPWM法、PPWM法、LH矩法等新型參數(shù)估計(jì)方法在洪水頻率分析中均取得了較為理想的效果，且優(yōu)于矩法、概率權(quán)重矩法等常規(guī)估參方法，其復(fù)雜計(jì)算可借助計(jì)算機(jī)程序得以簡化。目前，國內(nèi)對于這些新型估參方法的研究相對較少，且多為單一種類的新型估參方法與傳統(tǒng)方法的對比研究，缺乏不同新型估參方法之間的橫向比較。而在實(shí)際工作中，統(tǒng)計(jì)特性穩(wěn)定、估計(jì)精度高、計(jì)算方便、能切實(shí)提高設(shè)計(jì)值估計(jì)精度、改善擬合效果的新型估參方法顯然更容易實(shí)現(xiàn)推廣應(yīng)用。

PPWM法和HPWM法均為基于概率權(quán)重矩法改進(jìn)的新型參數(shù)估計(jì)方法，均減小了實(shí)測洪水序列中較小流量值對擬合結(jié)果的影響。本文通過蒙特卡洛試驗(yàn)，對PPWM法和HPWM法的統(tǒng)計(jì)特性進(jìn)行了模擬分析與比較，以國家重點(diǎn)水文站神木水文站多年實(shí)測年最大流量系列為實(shí)例，分別對兩種方法的適用性進(jìn)行檢驗(yàn)，并比較兩種方法的優(yōu)劣。結(jié)果表明：在統(tǒng)計(jì)特性方面，HPWM法的無偏性和有效性均顯著優(yōu)于PPWM法，且估計(jì)結(jié)果穩(wěn)定；在擬合效果方面，與普通概率權(quán)重矩相比，提高PWM的階數(shù)(HPWM法)或增加刪失水平(PPWM法)均可顯著改善實(shí)測洪水序列的擬合效果，PPWM方法在P=60%～98%時的設(shè)計(jì)值的累積誤差更?。辉谟?jì)算難度方面，PPWM法估計(jì)量公式推導(dǎo)及計(jì)算難度較大，HPWM法計(jì)算時僅需改變階數(shù)，計(jì)算相對簡便，而這兩種方法對實(shí)例的擬合效果相當(dāng)。

此外，筆者在此前的研究中發(fā)現(xiàn)[21]：HPWM法與PPWM法對陜北地區(qū)其他水文站年徑流序列的頻率分析結(jié)果均遜色于普通概率權(quán)重矩法，且PWM階數(shù)越高/刪失水平越高，擬合效果越差；洪峰流量偏小的站點(diǎn)，HPWM法和PPWM法對擬合效果的改善相對不顯著。

綜上所述，HPWM法與PPWM法估參均適宜用于洪水頻率分析計(jì)算，且擬合效果基本相當(dāng)，適宜用于洪峰流量較大的站點(diǎn)，在實(shí)際工作中更推薦使用計(jì)算相對簡便的HPWM法進(jìn)行頻率分析。