部分概率權重矩法和高階概率權重矩法在洪水頻率分析中的比較研究?剻?

李 揚

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

0 引 言

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

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

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

因此，本文擬對同樣基于概率權重矩法改進的PPWM法和HPWM法進行統計特性和實際應用的比較研究。采用蒙特卡洛試驗，對這兩種參數估計方法的統計特性進行研究和比較，并通過實例對其適用性及擬合效果進行檢驗，以期為實際應用中優選參數估計方法提供依據。

1 部分概率權重矩和高階概率權重矩

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

(1)

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

(2)

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

(3)

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

當k≠0時，GEV分布的PWM為：

(4)

其無偏估計量br為：

(5)

1.1 部分概率權重矩(PPWM)

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

(6)

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

(7)

可推出：

(8)

令：

(9)

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

可推求出各參數的估計量表達式為：

(10)

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

(11)

(12)

(13)

1.2 高階概率權重矩(HPWM)

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

(14)

令：

(15)

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

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

可推求出各參數的估計量表達式為：

(16)

(17)

(18)

(19)

2 蒙特卡洛試驗

2.1 試驗方案

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

2.2 評價標準

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

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

(20)

標準誤差SE反映了估計量的有效性，對同一總體參數的估計量中，SE更小的估計量更有效，即SE的值越小，估計量的有效性越好。若統計試驗次數為N，則SE的計算公式為：

(21)

均方根誤差RMSE是估計值與實際值偏差平方和與試驗次數比值的平方根，用來衡量估計值與實際值之間的偏差，RMSE的值越小，表明估計的偏差越小，有效性越好。若統計試驗次數為N，則RMSE的計算公式為：

(22)

2.3 試驗結果

按照以上試驗方案，用Matlab R2018a編程進行蒙特卡洛試驗，結果如表3～表4及圖1～圖6所示。

可以發現如下規律：

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

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

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

表3 PPWM法蒙特卡洛試驗結果Tab.3 Monte Carlo test results of PPWM method

表4 HPWM法蒙特卡洛試驗結果Tab.4 Monte Carlo test results of HPWM method

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

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

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

圖4 不同階數下的HPWM法估計值BiasFig.4 Bias of HPWM estimators in different orders

圖5 不同階數下的HPWM法估計值SEFig.5 SE of HPWM estimators in different orders

圖6 不同階數下的HPWM法估計值RMSEFig.6 RMSE of HPWM estimators in different orders

3 實例應用

3.1 資料及研究區概況

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

3.2 分布參數估計與頻率曲線繪制

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

3.3 擬合效果評價

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

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

表5 GEV分布參數估計結果Tab.5 Parameter estimates of GEV distribution

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

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

(23)

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

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

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

綜合蒙特卡洛試驗和實例應用的結果可知，HPWM方法的統計特性更加優良且穩定，實例應用中HPWM法和PPWM法的擬合效果基本相當，PPWM方法在P=60%～98%時的設計值的累積誤差更小，可根據實際情況靈活選用。

4 結 論

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

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

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

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