基于貝葉斯MCMC方法的洪水頻率分析及不確定性評估

吳云標，遲藝俠

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

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

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

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

1 研究方法

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

1)貝葉斯公式

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

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

2)先驗分布

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

3)似然函數(shù)

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

4)后驗分布

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

5)參數(shù)估計

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

1.2 洪水設計值計算

假設洪水觀測值x已知，觀測值為x時θ的后驗分布用f(θ|x)表示，如果用z表示未來洪水設計值，則z的預測密度可表示為

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

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

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

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

由式(4)得

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

1.3 MCMC方法

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

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

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

3)計算接受概率αi

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

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

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

1.4 GEV分布

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

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

2 實例分析

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

2.1 洪水頻率分布選擇

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

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

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

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

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

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

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

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

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

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

2.3 擬合檢驗

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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