基于廣義Pareto分布的洪水序列頻率分析*

趙玲玲，陳子燊，劉昌明，許揚(yáng)生

(1. 廣州地理研究所廣東省地理空間信息技術(shù)與應(yīng)用公共實(shí)驗(yàn)室，廣東 廣州510070； 2. 中國科學(xué)院地理科學(xué)與資源研究所中國科學(xué)院陸地水循環(huán)及 地表過程重點(diǎn)實(shí)驗(yàn)室，北京100101； 3. 中山大學(xué)地理科學(xué)與規(guī)劃學(xué)院，廣東 廣州510275； 4. 北京師范大學(xué)水科學(xué)研究院，北京1008755； 5.廣東省水文局，廣東 廣州 510050)

變化環(huán)境下導(dǎo)致極端事件呈現(xiàn)增加趨勢，其中洪水屬于典型的災(zāi)害事件?；谟邢薜乃挠^測數(shù)據(jù)中發(fā)現(xiàn)盡可能多的洪水的規(guī)律，提高推算洪水重現(xiàn)水平的精度，對防洪工程規(guī)劃設(shè)計(jì)和災(zāi)害風(fēng)險(xiǎn)評估有重要作用。超限頻率分析是極值統(tǒng)計(jì)建模理論的重要組成部分[1]，國內(nèi)外已有不少探索與研究。研究人員對超定量樣本的頻率分布線型的適用性做了較多研究，如指數(shù)分布[2-3]、Gamma分布和Weibull分布[4]。之后研究集中在廣義Pareto分布(GPD)[5-8]。而GPD應(yīng)用于實(shí)際數(shù)據(jù)集的成功與否很大程度上取決于參數(shù)估計(jì)過程。Hosking等[9]使用極大似然法(ML)和概率權(quán)重矩法(PWM)估計(jì)GPD模型的參數(shù)，Zhang[10]提出了基于GPD的參數(shù)估計(jì)似然矩法。Bermudez等[11-12]詳細(xì)介紹了多數(shù)情況下使用的最大似然(ML)、概率加權(quán)矩(PWM)和矩法(MOM)參數(shù)估計(jì)方法及其優(yōu)缺點(diǎn)。王劍峰等[13]用常規(guī)線性矩法和改進(jìn)線性矩法對廣義Pareto分布參數(shù)進(jìn)行估計(jì)，并對比分析了超定量序列頻率。周長讓等[14]采用高階概率權(quán)重矩估計(jì)其分布參數(shù)，統(tǒng)計(jì)試驗(yàn)表明該方法具更高的參數(shù)估計(jì)精度，估參結(jié)果對應(yīng)的GPD曲線能較好地?cái)M合稀遇頻率洪水段的經(jīng)驗(yàn)頻率分布。但在洪水超限頻率分析中如何選取閾值這一關(guān)鍵問題、超定量數(shù)和超定量分布的擬合優(yōu)度檢驗(yàn)等還需要更多的探索與實(shí)踐[7]。

本文以山區(qū)中小流域日流量過程為例，深入探討廣義Pareto分布(GPD)模型的閾值選擇方法，洪水序列超定量數(shù)檢驗(yàn)、擬合優(yōu)度檢驗(yàn)，最后對GPD、GEV分布與P-III型分布推算的設(shè)計(jì)水平加以對比。

1 極值分布模型

單變量極值分布模型有2種常用方法：一是分組區(qū)塊最大值模型(block maximum group of models，BM)。首先是對所得到的數(shù)據(jù)進(jìn)行分塊，常用年最大值方法(annual maximal series，AMS)采用年最大值樣本作為建立模型的觀測值。BM模型的要求數(shù)據(jù)樣本獨(dú)立同分布(IID)。而對洪水序列而言，常常一年內(nèi)出現(xiàn)多次洪水樣本，所以按年最大值抽樣后會出現(xiàn)多個洪峰流量遠(yuǎn)大于枯水年份出現(xiàn)的最大值的現(xiàn)象。按年最大值抽樣顯然會造成洪水信息利用不充分；二為峰量超閾值模型(peaks over threshold models，POT)，其選取超特定閾值抽取樣本，該方法可獲取較多的洪水極值序列建立模型。POT模型要求滿足超定量發(fā)生的時間服從泊松分布，且彼此相互獨(dú)立服從GPD(generalized Pareto distribution)分布)[8]。

1.1 廣義Pareto分布

設(shè)序列{xn} 的分布函數(shù)為F(x)，定義Fμ(y)為隨機(jī)變量X超過閾值μ的條件分布函數(shù)：

Fμ(y)=

P(X-μ≤x|X>μ)=

?F(x)=Fμ(y)(1-F(μ))+F(μ)

(1)

研究表明[18-20]，當(dāng)μ閾值足夠大時，條件超量分布函數(shù)Fμ(y) 收斂于廣義Pareto分布，累積分布函數(shù)為：

Fμ(y)≈G(x,ξ,σ,μ)=

(2)

式中，ξ、σ、μ分別為形態(tài)參數(shù)、尺度參數(shù)和位置參數(shù)。當(dāng)ξ=0時，GPD對應(yīng)于指數(shù)分布,為Pareto I型分布；當(dāng)ξ<0時為Pareto II型分布，x∈[μ,∞)；當(dāng)ξ>0，為Pareto III型分布(短尾型)。有關(guān)研究證明了超定量(X-μ)數(shù)服從泊松分布[17]。

GPD模型T年一遇的分位數(shù)xT為：

(3)

1.2 超定量洪水序列閾值確定

閾值μ的合理確定是GPD模型參數(shù)ξ和σ正確估計(jì)的前提。μ取值過高，超限數(shù)據(jù)量少，使估計(jì)出來的參數(shù)方差很大；相反μ取值過低，則難以保證分布的收斂性，估計(jì)偏差較大。閾值選取基本原則：對于每次洪水過程，選取洪峰流量最大值；不同場次洪峰取樣時，各場次洪峰發(fā)生時間間隔要求大于流域的匯流時間，以保證不同洪水之間相互獨(dú)立。采用以下4種相結(jié)合的方法確定μ：

1)繪制樣本的經(jīng)驗(yàn)平均超過函數(shù)圖

令X(1)>X(2)>>X(n)，樣本的經(jīng)驗(yàn)平均超過函數(shù)定義為[21]：

(4)

其中

k=min{i|Xi>μ}

繪制的平均超過函數(shù)圖即為點(diǎn)(μ,e(μ))構(gòu)成的曲線，選取較大的μ作為閾值，它使得當(dāng)x≥μ時e(x)為近似線性函數(shù)。當(dāng)x≥μ時平均超過函數(shù)圖曲線向上傾斜，表明點(diǎn)據(jù)服從形狀參數(shù)ξ為負(fù)的GPD分布，屬于II型分布。當(dāng)x≥μ時曲線向下傾斜，表明數(shù)據(jù)源自于尾部較短的分布，當(dāng)x≥μ時曲線是水平的，表明該數(shù)據(jù)服從指數(shù)分布。因此，如果某個閾值μ后的e(n)趨向于新的線性變化時，可選取這個值為閾值。

2)Anderson-Darling(AD)檢驗(yàn)

AD檢驗(yàn)屬于平方根類經(jīng)驗(yàn)分布函數(shù)統(tǒng)計(jì)檢驗(yàn)，采用AD檢驗(yàn)GPD模型的超定量樣本的經(jīng)驗(yàn)分布和理論分布的擬合優(yōu)度時，AD檢驗(yàn)在小頻率的上尾部區(qū)域賦予了更多權(quán)重。

(5)

3)泊松分布檢驗(yàn)

為了保證所選樣本的獨(dú)立性，選擇一定的閾值區(qū)間，在顯著水平0.05下，分別對閾值區(qū)間內(nèi)不同閾值的超限數(shù)采用χ2假設(shè)檢驗(yàn)其是否服從泊松分布：

P(x=k)=e-λλk/k！，k=0，1，2，

其中λ為平均發(fā)生超限的頻次。對服從泊松分布樣本的閾值再根據(jù)其它擬合優(yōu)度檢驗(yàn)指標(biāo)作進(jìn)一步的篩選。

4)初始閾值與洪水序列的確定

以枯水年最大洪峰流量值為初始閾值，根據(jù)流域匯流時間形成閾值區(qū)間，采用不同的匯流時間間隔分別抽取不同間隔時段內(nèi)不同序列的洪峰閾值，在顯著水平0.05下，對各超定量樣本數(shù)的做獨(dú)立性檢驗(yàn)。對確定的洪水序列樣本估計(jì)分布參數(shù)值及其擬合優(yōu)度檢驗(yàn)值。

通過以上4種相結(jié)合的方法對于確定最終采用的超閾值洪水序列樣本滿足獨(dú)立性，閾值μ為優(yōu)選值。

2 實(shí)例研究

2.1 研究背景與基本數(shù)據(jù)

選取廣東強(qiáng)降水地區(qū)典型中小流域-曹江流域作為研究對象。曹江是粵西獨(dú)流入海鑒江的一級支流，發(fā)源于高州馬貴鎮(zhèn)山心村海拔1 141 m的藍(lán)蓬嶺，中上游雨量充沛，是廣東的暴雨高區(qū)之一，流域坡降大匯流時間短，洪水陡漲陡落，導(dǎo)致洪災(zāi)頻發(fā)，對當(dāng)?shù)卦斐蓢?yán)重的生命威脅與經(jīng)濟(jì)損失。流域多年平均年雨量2 160 mm，最大平均年雨量為3 150 mm。曹江流域出口斷面大拜水文站集水面積394 km2。

利用1967～2013年共47年大拜站逐日流量觀測序列分析曹江流域洪水的超閾值頻率分布特征。流量序列日最大值778 m3/s，最小值僅0.46 m3/s，日平均流量19.4 m3/s。對大拜站流量序列的洪水過程統(tǒng)計(jì)表明，通過大拜站的洪水平均傳播時間絕大多數(shù)為1～3 d，最長歷時可達(dá)10 d(2010年9月21～30日)。

2.2 洪水閾值

采用以下步驟確定曹江流域洪峰流量閾值：

1)在平均超過函數(shù)圖(圖1)內(nèi)流量序列在100 m3/s左右存在線性變化折點(diǎn)，隨之曲線向上傾斜，表明點(diǎn)據(jù)服從形狀參數(shù)為負(fù)的GPD分布，可考慮選取線性變化折點(diǎn)值為參考閾值。

圖1 超限樣本的經(jīng)驗(yàn)平均超過函數(shù)圖Fig.1 Definite mean excess function chart

2)以枯水年(2007年)最大洪峰流量值74 m3/s為初始閾值μ0，形成閾值區(qū)間：μi=μ0+(i-1)×10，i=1,,15，閾值選擇范圍為74～214 m3/s。

3)按不同洪水序列之間相互獨(dú)立的要求，根據(jù)流域匯流時間短的特點(diǎn)，以3 d時間為初始間隔時段，以1 d為步長增量逐步增加至10 d時間間隔，按照步驟1)分別抽取不同間隔時段內(nèi)不同序列的洪峰閾值，在顯著水平0.05下，對各超定量樣本數(shù)的泊松分布加以χ2檢驗(yàn)(表1中原假設(shè)，H0=0：樣本服從泊松分布)；為節(jié)省版面，文中僅列出采用參數(shù)概率權(quán)重矩估計(jì)和極大似然估計(jì)法4，6，8和10 d時間間隔的洪水序列樣本的參數(shù)估計(jì)值及其擬合優(yōu)度檢驗(yàn)值。

4)超限抽樣系列的AD檢驗(yàn)的P值都遠(yuǎn)大于α=0.05的臨界值，超定量樣本的頻率分布與理論CDF擬合良好，表明該超限樣本符合GPD分布。除去不符合泊松分布的超限樣本和λ<1的超限樣本后發(fā)現(xiàn)，在8種采樣時間間隔的洪峰序列中閾值為104 m3/s的超限樣本的PAD值最大，其中采樣時間間隔中間隔8 d的PAD值和PPCC(概率圖位相關(guān)系數(shù))值最大，RMSE(均方根誤差)最小。因此104 m3/s可作為大拜站洪水序列的日洪峰流量序列的優(yōu)選閾值(表1)。

2.3 GPD參數(shù)估計(jì)與擬合優(yōu)度檢驗(yàn)結(jié)果

極值分布模型的參數(shù)估計(jì)是統(tǒng)計(jì)分析的關(guān)鍵點(diǎn)之一。不同的參數(shù)估計(jì)方法推算的分布參數(shù)直接影響極值重現(xiàn)水平。因此本文使用具有統(tǒng)計(jì)特性良好的概率權(quán)重矩法(PWM)和極大似然法(ML)[6]估計(jì)模型的參數(shù)，對各超限樣本的GPD重現(xiàn)水平的推算結(jié)果進(jìn)一步采用了PPCC和RMSE作為擬合優(yōu)度檢驗(yàn)指標(biāo)，主要結(jié)果如下：

1)PPCC值均大于0.984，表明各個超定量樣本點(diǎn)據(jù)與理論分布值的相關(guān)關(guān)系達(dá)0.984以上。同時PWM參數(shù)估計(jì)方法得到的GPD模型的均方誤差RMSE小于ML擬合的模型誤差，其中λ=2.0-3.0估計(jì)的指標(biāo)更可靠，擬合優(yōu)度檢驗(yàn)結(jié)果見表1。

2)模型的形態(tài)參數(shù)為負(fù)值，表明曹江洪水序列服從Pareto II型分布。

前述結(jié)果顯示，GPD模型確定超定量洪水是根據(jù)多個指標(biāo)綜合分析確定的動態(tài)過程。

2.4 不同概率分布模型對比分析

對比GPD和GEV兩種極值分布及P-III型分布推算參數(shù)的分布函數(shù)擬合指標(biāo)值，P-III型分布參數(shù)估計(jì)使用常規(guī)矩法(OME)和線性矩法(L-M)。三種模型最優(yōu)模型參數(shù)和擬合優(yōu)度檢驗(yàn)指標(biāo)見表2，各分布模型推算的洪水重現(xiàn)水平見表3和圖2～4。

三種概率分布的擬合優(yōu)度指標(biāo)對比顯示，超限量抽樣在滿足超限量數(shù)服從泊松分布，超限彼此相互獨(dú)立條件下構(gòu)建的GPD模型，在介于2.0～3.0之間構(gòu)建的GPD模型精度優(yōu)于GEV和P-III型，圖2顯示，大拜站超定量洪水頻率曲線圖上的樣本點(diǎn)與理論曲線非常吻合，尤以PWM參數(shù)估計(jì)推算的GPD模型最佳。以8 d為抽樣間隔的擬合優(yōu)度檢驗(yàn)指標(biāo)為例，閾值74～124 m3/s范圍內(nèi)，采用PWM參數(shù)估計(jì)方法擬合的GPD模型PPCC值大于等于0.984，以閾值104 m3/s對應(yīng)的PAD和PPCC值最大，RMSE也明顯小于P-III分布和GEV分布。

同頻率設(shè)計(jì)值對比結(jié)果表明，設(shè)計(jì)頻率小于2%(重現(xiàn)期50年)時，GPD模型設(shè)計(jì)值介于GEV和P-III之間，隨著設(shè)計(jì)頻率增大，GPD設(shè)計(jì)值超過GEV和P-III分布的設(shè)計(jì)值(表3)。此或反映了自然流域洪水過程的差異性，其原因需要通過更多的實(shí)證加以歸納說明。

表1 GPD模型的閾值、參數(shù)估計(jì)與擬合優(yōu)度檢驗(yàn)結(jié)果Table 1 The selected threshold values, estimated parameters and goodness-fit test of GPD models

表2 最優(yōu)GPD、GEV、P-III分布參數(shù)與擬合優(yōu)度檢驗(yàn)對比Table 2 Comparison of parameters and indicators of goodness-fit test among the optimal GPD, GEV and Pearson Type III models

表3 三種概率分布函數(shù)洪水重現(xiàn)水平值Table 3 Return levels of flood discharge of three probability distributions m3/s

圖2 大拜站超定量洪水頻率曲線圖Fig.2 POT flood frequency curvesforDabaistation

圖3 大拜站年最大洪水GEV分布頻率曲線Fig.3 Annual maximal flood frequency curves basedGEV distribution forDabaistation

圖4 大拜站年最大洪水P-Ⅲ分布頻率曲線Fig.4 Annual maximal flood frequency curves based Pearson Type III distribution for Dabai station

3 結(jié) 論

1)GPD模型的形態(tài)參數(shù)表明洪水頻率分布屬于II型，與短尾分布III型不同，II型表明密度分布函數(shù)峰值右側(cè)分布曲線與橫坐標(biāo)之間的漸進(jìn)性而無切點(diǎn)，難于確定洪水的上限值。此是否反映了山地流域不同時段不同區(qū)域土壤含水量差異大和洪水的產(chǎn)匯流過程不確定性的自然屬性，有待更多的實(shí)例研究。

2)確定超定量洪水GPD 模型是動態(tài)擇優(yōu)過程，對多個滿足GPD模型要求的閾值，需要通過超定量數(shù)的泊松分布檢驗(yàn)和超限樣本的擬合優(yōu)度等綜合評判后構(gòu)建相對最優(yōu)GPD模型。

3)GPD模型推算的設(shè)計(jì)洪水精度普遍優(yōu)于GEV分布和P-III型分布的推算成果。

4)不同參數(shù)估計(jì)方法對于極值模型參數(shù)推算精度有較大影響。POT樣本由PWM法估計(jì)參數(shù)的極值模型精度高于ML法推算的結(jié)果。