吳云標(biāo)，遲藝俠

(1.河海大學(xué)水文水資源學(xué)院，江蘇南京210098；2.河海大學(xué)文天學(xué)院，安徽馬鞍山243031)

近幾十年來，洪水事件被認(rèn)為是世界上最普遍、分布最廣泛的自然災(zāi)害之一，洪水災(zāi)害是影響經(jīng)濟(jì)社會(huì)可持續(xù)發(fā)展的障礙性因素之一[1]。因此，洪水頻率分析已成為水文學(xué)中的一個(gè)研究熱點(diǎn)。在洪水頻率分析中，洪水設(shè)計(jì)值的計(jì)算和不確定性評(píng)估是水利工程規(guī)劃和水資源管理的一個(gè)重要課題[2]。

極值理論(EVT)是研究極端事件最為成功的方法之一，在水文、氣象、地震、保險(xiǎn)、金融等領(lǐng)域得到廣泛的應(yīng)用[3]。在極值理論中，獨(dú)立同分布隨機(jī)變量的區(qū)組最大值分布的漸近分布為廣義極值分布(GEV)[3-5]。傳統(tǒng)的GEV分布參數(shù)估計(jì)有極大似然估計(jì)(MLE)法、矩法、L-矩法等，其中極大似然估計(jì)法最為常用[6-7]。Martins等[8]研究表明，小樣本的極大似然估計(jì)不穩(wěn)定，形狀參數(shù)的估計(jì)有時(shí)會(huì)偏小，導(dǎo)致分位數(shù)與觀測(cè)值的均方根誤差偏大。Coles[4]也指出極大似然估計(jì)法雖然是從考慮不確定性出發(fā)，可以用來計(jì)算參數(shù)的置信區(qū)間，但由于極大似然估計(jì)量以漸近正態(tài)為基礎(chǔ)，除非樣本足夠大，否則不足以說明這種不確定性。近年來，國內(nèi)外圍繞小樣本參數(shù)估計(jì)和不確定性評(píng)估的研究取得了較大的進(jìn)展，其中貝葉斯估計(jì)法最具有代表性[9]。與傳統(tǒng)的參數(shù)估計(jì)法中將分布參數(shù)作為確定性未知常量不同，貝葉斯估計(jì)法將未知參數(shù)看成隨機(jī)變量，在對(duì)總體分布參數(shù)進(jìn)行統(tǒng)計(jì)推斷時(shí)，除了樣本提供的信息外，還需在抽樣前確定總體分布參數(shù)的先驗(yàn)分布，即使樣本序列較短，也能獲得比較精確的推斷結(jié)果[9-10]。近年來，貝葉斯估計(jì)法已逐步被引入到洪水頻率分析中。Kuczera[11]根據(jù)貝葉斯理論，運(yùn)用重要性抽樣法搜索參數(shù)的后驗(yàn)狀態(tài)空間計(jì)算洪水的置信區(qū)間。O'Connell等[12]利用貝葉斯理論將歷史和古水文信息等考慮到洪水頻率分析中，發(fā)現(xiàn)古水文信息能有效減少分位數(shù)估計(jì)偏差且能縮小設(shè)計(jì)洪水估計(jì)的不確定性。Liang等[2]基于貝葉斯理論，將模型不確定性和參數(shù)不確定性同時(shí)納入洪水頻率分析，對(duì)分位數(shù)的不確定性進(jìn)行了有效估計(jì)。

本文以GEV分布為洪水頻率分布線型，利用基于貝葉斯理論的MCMC模擬方法，研究洪水高分位數(shù)設(shè)計(jì)值的推求方法，得到相應(yīng)的置信區(qū)間，從而對(duì)估計(jì)值的不確定性進(jìn)行定量評(píng)價(jià)，并以洞庭湖流域四個(gè)水文站點(diǎn)洪水頻率分析為實(shí)例，驗(yàn)證方法的有效性。

1 研究方法

1.1 貝葉斯參數(shù)估計(jì)

1)貝葉斯公式

假設(shè)樣本x=(X1,…,Xn)的密度函數(shù)屬于參數(shù)族F={f(x;θ):θ∈Θ}，其中參數(shù)θ是一個(gè)未知的常數(shù)。在經(jīng)典的統(tǒng)計(jì)推斷中，在得到樣本x之前，對(duì)參數(shù)θ是一無所知的。在貝葉斯框架下，總體參數(shù)θ被當(dāng)作隨機(jī)變量，根據(jù)歷史資料或經(jīng)驗(yàn)，給出參數(shù)的概率密度函數(shù)π(θ)，稱為θ的先驗(yàn)密度。貝葉斯公式將先驗(yàn)信息與樣本信息進(jìn)行耦合，如式(1)。

其中：π(θ)和f(θ|x)分別為參數(shù)θ的先驗(yàn)和后驗(yàn)密度函數(shù)；f(x|θ)為樣本x=(X1,…,Xn)的似然函數(shù)；是θ的參數(shù)空間。

2)先驗(yàn)分布

先驗(yàn)分布π(θ)的確定是貝葉斯分析的關(guān)鍵，目前主要采用一些經(jīng)驗(yàn)的方法。如：在無歷史數(shù)據(jù)以及經(jīng)驗(yàn)的情況下，按照貝葉斯理論的無信息先驗(yàn)分布的原則，可認(rèn)為θ服從均勻分布；若參數(shù)θ傾向于取較小的實(shí)數(shù)值，則可用具有較大方差的正態(tài)分布表示[4]。

3)似然函數(shù)

Xi=(i=1,…,n)相互獨(dú)立，樣本x=(X1,…,Xn)的似然函數(shù)f(x|θ)可用式(2)計(jì)算。

4)后驗(yàn)分布

后驗(yàn)分布f(θ|x)由先驗(yàn)分布π(θ)和似然函數(shù)f(x|θ)通過式(1)確定。

5)參數(shù)估計(jì)

在貝葉斯理論中，總體分布參數(shù)θ的統(tǒng)計(jì)完全由后驗(yàn)分布推斷，不需要借助極大似然估計(jì)的漸近正態(tài)性來得到參數(shù)估計(jì)量的漸近分布。取后驗(yàn)分布的均值作為θ的估計(jì)值。

1.2 洪水設(shè)計(jì)值計(jì)算

假設(shè)洪水觀測(cè)值x已知，觀測(cè)值為x時(shí)θ的后驗(yàn)分布用f(θ|x)表示，如果用z表示未來洪水設(shè)計(jì)值，則z的預(yù)測(cè)密度可表示為

與其他的預(yù)測(cè)方法相比，其優(yōu)點(diǎn)在于它包含了反映模型不確定的f(θ|x)，以及反映未來觀測(cè)值變異性的f(z|θ)。

由式(3)可得洪水的未來分布

其中含有參數(shù)的不確定性和未來觀測(cè)值的隨機(jī)性。解方程可得m年重現(xiàn)水平的洪水設(shè)計(jì)值（即分位數(shù)）。

方程(5)中后驗(yàn)分布的計(jì)算較復(fù)雜，即便是用數(shù)值積分方法計(jì)算也比較困難。近年來模擬方法的快速發(fā)展為該問題的解決提供了新途徑，Markov Chain Monte Carlo（MCMC）法為其中有效方法之一[3]。本文采用MCMC法模擬產(chǎn)生服從后驗(yàn)分布的隨機(jī)樣本，去除前k個(gè)不穩(wěn)定的樣本，將剩余序列θk+1,θk+2,…,θn看作f(θ|x)的觀測(cè)值。

由式(4)得

再由數(shù)值方法求解得式(5)的解。

1.3 MCMC方法

MCMC方法的思想是通過建立Markov鏈模擬產(chǎn)生服從后驗(yàn)分布的隨機(jī)樣本，從而模擬樣本估計(jì)的后驗(yàn)分布[3]。該方法通過迭代產(chǎn)生模擬序列θ0,θ1,θ2,…，其中θ0為任意初始值，θi+1由條件分布q(·|θi)產(chǎn)生，即θi+1只依賴于當(dāng)前的θi，與前面的θ0,θ1,θ2,…,θi-1無關(guān)。不同的抽樣方法形成不同的Markov鏈，其中Metropolis--Hasting算法是目前應(yīng)用最為廣泛的MCMC抽樣方法之一，其基本步驟如下：

1)確定參數(shù)的初始值θ0，選定建議分布q(·|θi)；

2)由q(·|θi)產(chǎn)生一個(gè)新的建議值θ*；

3)計(jì)算接受概率αi

4)以概率αi接受θ*為下一個(gè)θi+1，即

其中μ為[0,1]均勻分布隨機(jī)數(shù)。

重復(fù)步驟2)~4)直到產(chǎn)生足夠多的樣本為止(樣本數(shù)為n)，去除前k個(gè)不穩(wěn)定的樣本，使得剩余的序列θk+1,θk+2,…,θn達(dá)到平穩(wěn)狀態(tài)，則此序列可認(rèn)為是后驗(yàn)分布的抽樣，并用來估計(jì)后驗(yàn)分布的數(shù)字特征。

1.4 GEV分布

GEV分布是Gumbel、Fréchet和Weibull分布3種極值分布的統(tǒng)一形式[3-5]，其分布函數(shù)為

其中：μ,σ,ξ分別表示GEV分布的位置、尺度、形狀參數(shù)，滿足

2 實(shí)例分析

洞庭湖流域位于長(zhǎng)江流域中下游，總面積約為26萬km2。流域內(nèi)徑流年際變化大，旱澇災(zāi)害發(fā)生幾率高[13]。本文選取流域內(nèi)四水(湘江、資水、沅江、澧水)主要水文控制站點(diǎn)(湘潭站、桃江站、桃源站、石門站)年最大洪水流量作為研究對(duì)象，其中：石門、桃江站流量資料為1951—2014年，桃源站為1953—2014年，湘潭站為1951—2012年。

2.1 洪水頻率分布選擇

在水文頻率分析中，洪水頻率通常采用GEV分布、皮爾遜III型(P-III)分布來分析。為從上述分布中選擇一種最佳概率分布，分別對(duì)洞庭湖流域內(nèi)四個(gè)站點(diǎn)年最大洪水流量用上述待選分布進(jìn)行擬合。為便于計(jì)算，年最大洪水流量單位統(tǒng)一為103m3·s-1。

采用均方根誤差(RMSE)檢驗(yàn)洪水頻率曲線的擬合效果，RMSE值越小表示擬合效果越好。計(jì)算如式(10)。

分別采用兩種待選分布對(duì)四個(gè)站點(diǎn)年最大洪水流量進(jìn)行擬合，其擬合效果檢驗(yàn)的RMSE值見表1，其中GEV分布參數(shù)估計(jì)采用極大似然法(MLE)，P-III分布參數(shù)估計(jì)采用矩法。由表1可知，除石門站P-III分布擬合的RMSE值略小于GEV外，其余3個(gè)站點(diǎn)GEV分布擬合的RMSE值均小于P-III分布，說明GEV分布對(duì)4個(gè)站點(diǎn)年最大洪水流量擬合有良好的適應(yīng)性。由于石門站P-III分布和GEV分布擬合結(jié)果接近，為便于分析，石門站仍采用GEV分布。

表1 各站點(diǎn)年最大洪水流量概率分布擬合效果RMSE檢驗(yàn)值Tab.1 RMSE values of probability distribution fitting effect of annual maximum flood flow distributions at each station

2.2 GEV分布參數(shù)的貝葉斯估計(jì)

分別對(duì)石門、桃源、桃江、湘潭4個(gè)水文站的年最大洪水流量建立GEV分布模型。采用Metropolis-Hastings算法產(chǎn)生隨機(jī)樣本，用去除前k個(gè)不穩(wěn)定樣本后的序列對(duì)GEV分布參數(shù)進(jìn)行統(tǒng)計(jì)推斷。選擇先驗(yàn)密度函數(shù)為π(μ,σ,ξ)=πμ(μ)πσ(σ)πξ(ξ)，其中：參數(shù)μ,σ,ξ相互獨(dú)立；πμ(μ),πσ(σ)和πξ(ξ)為均值為，方差分別為vμ,vσ,vξ正態(tài)分布的概率密度函數(shù)。為使密度函數(shù)平坦，選擇足夠大的方差，按照文獻(xiàn)[4]的建議，選取vμ=vσ=104,vξ=102。參數(shù)μ,σ,ξ的建議分布為各自坐標(biāo)軸上的隨機(jī)游動(dòng)，即μ*=μ+εμ,σ*=σ+εσ,ξ*=ξ+εξ,其中εμ,εσ,εξ為均值為0，方差分別為ωμ,ωσ,ωξ的正態(tài)隨機(jī)數(shù)。在實(shí)驗(yàn)中，通過調(diào)整ωμ,ωσ,ωξ的值，使接受概率αi落在0.2~0.5之間[14]。為使MCMC快速收斂，本文選擇GEV模型參數(shù)的極大似然估計(jì)值作為參數(shù)μ,σ,ξ的Metropolis-Hastings抽樣初始值。

以石門站為例。圖1為石門站年最大洪水流 量GEV模型參數(shù)經(jīng)過10 000次迭代產(chǎn)生的MCMC序列圖。從圖中可以看出，由極大似然估計(jì)值作為初始值生成的序列收斂很快，均在初始值附近就趨于穩(wěn)定，說明極大似然估計(jì)值作為抽樣的初始值合理，且算法效率高。為確保序列的平穩(wěn)性，去除前500個(gè)樣本，將剩余序列的模擬值作為后驗(yàn)分布的觀測(cè)值。

圖1 石門站GEV分布參數(shù)MCMC模擬Fig.1 MCMC simulation for parameters of GEV model at Shimen station

表2為4個(gè)站點(diǎn)基于貝葉斯理論MCMC法的GEV分布參數(shù)后驗(yàn)分布的統(tǒng)計(jì)特征(置信水平為95%)和參數(shù)抽樣的初始值(MLE值)。與傳統(tǒng)參數(shù)的估計(jì)方法相比，貝葉斯法不僅給出了參數(shù)的估計(jì)值，還給出了參數(shù)估計(jì)值的置信區(qū)間，從而能量化模型參數(shù)估計(jì)的不確定性。

表2 各站點(diǎn)GEV分布參數(shù)估計(jì)Tab.2 Estimation parameters of GEV distribution at each station

2.3 擬合檢驗(yàn)

以石門站為例。圖2給出了石門站年最大洪水流量的GEV模型貝葉斯估計(jì)分位數(shù)擬合效果圖。從圖中可以看出，樣本分位數(shù)和理論分位數(shù)均分布在45゜線附近，表明GEV模型與實(shí)測(cè)樣本序列擬合效果較好。

對(duì)4個(gè)站點(diǎn)擬合結(jié)果分別采用均方根誤差(RMSE)、Kolmogrov--Smirnov檢驗(yàn)(K-S檢驗(yàn))[4]進(jìn)行擬合優(yōu)度檢驗(yàn)。表3給出的是四個(gè)站點(diǎn)擬合優(yōu)度檢驗(yàn)結(jié)果。從表3中可以看出，貝葉斯估計(jì)法和極大似然估計(jì)法均通過K-S檢驗(yàn)(顯著水平α=0.05)。說明兩種估計(jì)方法均適用于GEV模型參數(shù)估計(jì)。對(duì)比各站點(diǎn)采用兩種參數(shù)估計(jì)方法所得的RMSE可以看出，石門、桃源、桃江站由貝葉斯估計(jì)得到的RMSE略大于極大似然估計(jì)，但是相差非常小，相差最大的石門站僅為0.065 1，說明貝葉斯估計(jì)和極大似然估計(jì)結(jié)果近似相同。而在湘潭站，貝葉斯估計(jì)得到的RMSE略小于極大似然估計(jì)，說明貝葉斯估計(jì)擬合略好于極大似然估計(jì)。總體來說，兩者擬合效果基本一致。

圖2 石門站GEV模型貝葉斯估計(jì)分位數(shù)圖Fig.2 Quantile plots for the GEV model at Shimen station

表3 各站點(diǎn)擬合優(yōu)度檢驗(yàn)結(jié)果Tab.3 Results of goodness-of-fit test of each station

2.4 重現(xiàn)水平估計(jì)

當(dāng)GEV分布模型建立后，T=1/p年重現(xiàn)期的重現(xiàn)水平x1-p可由下式計(jì)算

分別將各站點(diǎn)模型參數(shù)μ，σ，ξ的模擬值代入式(11)，可得到相應(yīng)于T=1/p年重現(xiàn)期的重現(xiàn)水平(設(shè)計(jì)洪水)后驗(yàn)分布樣本。

圖3給出了石門站各典型重現(xiàn)期下的年最大洪水流量后驗(yàn)密度估計(jì)。根據(jù)后驗(yàn)密度估計(jì)，計(jì)算流域內(nèi)四個(gè)站點(diǎn)各重現(xiàn)期下的設(shè)計(jì)洪水流量及95%的置信區(qū)間。其結(jié)果如表4。

表4位采用貝葉斯法對(duì)4個(gè)站點(diǎn)年最大洪水流量不同重現(xiàn)水平的估計(jì)。由表4可以看出，貝葉斯法估計(jì)的洪水設(shè)計(jì)值均小于置信區(qū)間的平均值，說明置信區(qū)間不關(guān)于設(shè)計(jì)值對(duì)稱。在實(shí)際中，由于大洪水資料有限，洪水設(shè)計(jì)值上限的不確定性通常大于下限的不確定性[10]。這說明貝葉斯法估計(jì)的結(jié)果與實(shí)際相符。此外，從表4還可以看出，重現(xiàn)期越大，對(duì)應(yīng)的設(shè)計(jì)值越大，相應(yīng)的置信區(qū)間越寬，表明其不確定性越大。

圖3 石門站各重現(xiàn)期年最大洪水流量的后驗(yàn)密度估計(jì)Fig.3 Posterior density plots of annual maximum flood flow under different return periods at Shimen station

表4 各站點(diǎn)年最大洪水流量的不同重現(xiàn)水平估計(jì)Tab.4 Different return level estimates of annual maximum flood flow at each station

3 結(jié) 論

以GEV分布作為洪水分布線型，利用基于Metropolis-Hastings抽樣的貝葉斯MCMC方法估計(jì)GEV分布參數(shù)，推求洪水設(shè)計(jì)值的點(diǎn)估計(jì)和區(qū)間估計(jì)，應(yīng)用于洞庭湖流域內(nèi)四個(gè)水文站點(diǎn)洪水頻率分析，得以下結(jié)論。

1)以GEV分布參數(shù)的極大似然估計(jì)值作為Metropolis-Hastings抽樣初始值可有效提高M(jìn)CMC方法收斂速度。

2)通過貝葉斯MCMC方法，可獲得洪水分布參數(shù)的后驗(yàn)分布和洪水設(shè)計(jì)值的后驗(yàn)分布。與傳統(tǒng)方法相比，本文方法不僅能得到洪水設(shè)計(jì)值的估計(jì)值，也可得到設(shè)計(jì)值的置信區(qū)間，從而可對(duì)估計(jì)結(jié)果的不確定性進(jìn)行定量評(píng)價(jià)。

3)基于貝葉斯法估計(jì)的洪水設(shè)計(jì)值小于置信區(qū)間的平均值，置信區(qū)間上限與估計(jì)值的距離大于置信區(qū)間下限與估計(jì)值的距離。這種不對(duì)稱性比傳統(tǒng)方法更貼近于實(shí)際，說明貝葉斯法估計(jì)進(jìn)一步提高了洪水頻率分析結(jié)果的可靠性。