基于貝葉斯理論的洪水頻率不確定性分析

張冬冬，魯 帆，劉少華

(1.長江水利委員會水文局，武漢 430010；2. 中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100038 )

洪水頻率分析主要是利用已有的實(shí)測洪水?dāng)?shù)據(jù)構(gòu)造洪水分布函數(shù)，通過分布函數(shù)計(jì)算一定重現(xiàn)期下的洪水設(shè)計(jì)值。洪水頻率分析中的難點(diǎn)在于極值事件的頻率計(jì)算，由于極值洪水相對于一般的洪水樣本，發(fā)生頻率較小，而在實(shí)際計(jì)算的過程中，極值洪水以及可以利用的歷史特大洪水觀測數(shù)據(jù)往往較為缺乏，影響分布函數(shù)參數(shù)的估計(jì)和高分位數(shù)計(jì)算，從而導(dǎo)致洪水頻率分析中存在較大的不確定性。如何定量的評估洪水頻率計(jì)算中線型參數(shù)以及推求設(shè)計(jì)值的不確定性，成為洪水頻率分析中的熱點(diǎn)問題。

本文選取大渡河流域典型站點(diǎn)洪水資料作為分析對象，選取廣義極值分布作為洪水分布的線型，利用貝葉斯原理，通過先驗(yàn)分布和似然函數(shù)構(gòu)造后驗(yàn)分布；利用MCMC模擬方法，根據(jù)抽樣值構(gòu)造一定頻率下的洪水設(shè)計(jì)值，并得到其相應(yīng)的置信區(qū)間，以定量描述洪水頻率分析的不確定性，并與傳統(tǒng)方法的置信區(qū)間進(jìn)行對比，從而對該方法的有效性和優(yōu)良性進(jìn)行評估。

1 貝葉斯洪水頻率分析的基本原理

貝葉斯估計(jì)的理論基礎(chǔ)是貝葉斯原理，主要是綜合考慮未知總體分布和樣本分布中的參數(shù)先驗(yàn)信息，根據(jù)信息求出樣本的后驗(yàn)分布，再根據(jù)后驗(yàn)分布推斷未知總體分布的參數(shù)[1]。貝葉斯估計(jì)與傳統(tǒng)的極大似然估計(jì)法、矩法和權(quán)函數(shù)法等參數(shù)估計(jì)方法的主要區(qū)別在于：貝葉斯估計(jì)將未知參數(shù)視為隨機(jī)變量，在利用樣本信息推斷總體分布參數(shù)之前，需要規(guī)定樣本參數(shù)的先驗(yàn)分布，先驗(yàn)分布的選取通常是利用已有的歷史洪水資料或利用決策者的經(jīng)驗(yàn)進(jìn)行設(shè)定。在洪水頻率計(jì)算中，能夠運(yùn)用一些有效的信息(歷史洪水等)，表示成先驗(yàn)分布的形式，通過最大化的利用有效信息，可以獲得較其他傳統(tǒng)參數(shù)估計(jì)方法更精確的推斷。與傳統(tǒng)估計(jì)方法相比，貝葉斯估計(jì)方法優(yōu)勢在于對模型不確定性的評估，貝葉斯方法利用大量隨機(jī)樣本模擬，概率分布的形式輸出模型的參數(shù)以及一定重現(xiàn)期下的設(shè)計(jì)值，概率分布相比與單一估計(jì)值包含了額外的信息，可以反映洪水頻率分析中的不確定性。近年來，許多學(xué)者嘗試在水文頻率分析研究中采用貝葉斯的估計(jì)方法。Wood等[2]最早提出了基于貝葉斯理論的水文模型參數(shù)估計(jì)方法的技術(shù)框架，Kuczera[3]根據(jù)貝葉斯理論同時(shí)探索了應(yīng)用分布參數(shù)不確定性問題，主要是利用重要抽樣方法分析了皮爾遜III型和對數(shù)皮爾遜III型曲線的分布參數(shù)。W Femandes等[4]認(rèn)為PMF的不確定性和先驗(yàn)分布線型有關(guān)，針對洪水頻率線型多數(shù)沒有考慮上下界問題，他采用貝葉斯方法和具有上界的分布函數(shù)相結(jié)合，用以福爾瑟姆水庫(美國)為例進(jìn)行洪水頻率分析。后驗(yàn)統(tǒng)計(jì)推斷中分母積分的計(jì)算是應(yīng)用貝葉斯方法的一個(gè)難點(diǎn)。為了簡化后驗(yàn)分布中的積分計(jì)算，可以通過選擇合理的先驗(yàn)分布族的方法來實(shí)現(xiàn)，但某些情況下，數(shù)值積分的方法無法計(jì)算參數(shù)較多且結(jié)構(gòu)復(fù)雜的后驗(yàn)分布。馬爾科夫鏈蒙特卡羅(MCMC)方法能夠解決貝葉斯后驗(yàn)分布高維積分的問題[5]。

假設(shè)洪水序列以及分布線型是已知的情況下，利用貝葉斯估計(jì)方法進(jìn)行洪水頻率分析的基本框架如圖1所示。

圖1 貝葉斯參數(shù)估計(jì)方法的基本原理Fig.1 Basic principle of bayesian parameters estimation

(1)先驗(yàn)分布。依據(jù)一定的原則確定總體分布參數(shù)的先驗(yàn)分布。在先驗(yàn)分布的選擇上，目前沒有一個(gè)統(tǒng)一的規(guī)范。但是有一些推薦的原則，例如，在沒有歷史數(shù)據(jù)以及經(jīng)驗(yàn)可以作為參考的情況下，可以根據(jù)參數(shù)的物理意義，采用取值范圍內(nèi)的均勻分布作為總體分布參數(shù)的先驗(yàn)分布。如果認(rèn)為參數(shù)在取值范圍內(nèi)更傾向取較小的值，推薦應(yīng)用方差較大的正態(tài)分布，相反，可以應(yīng)用方差較小的正態(tài)分布。

(2)似然函數(shù)。以θ表示總體分布參數(shù)向量，可以用π(θ)表示其先驗(yàn)分布的密度。由于Xi是相互獨(dú)立的，所以樣本x=(X1,X2,…,Xn)的似然函數(shù)f(x|θ)可用下式計(jì)算：

(1)

(3)后驗(yàn)分布。由貝葉斯原理，后驗(yàn)分布可以由先驗(yàn)密度π(θ)以及式(1)確定的似然函數(shù)計(jì)算得到，如下式：

(2)

(4)總體分布的參數(shù)估計(jì)。根據(jù)式(2)求得總體分布參數(shù)θ的后驗(yàn)分布，貝葉斯統(tǒng)計(jì)推斷與傳統(tǒng)的參數(shù)估計(jì)相比，最大的不同是貝葉斯給出的總體分布參數(shù)θ的分布函數(shù)。基于分布函數(shù)可以對θ的進(jìn)行統(tǒng)計(jì)推斷，一般情況下，θ的估計(jì)值有兩種選取方法，第一種方法是將后驗(yàn)分布的50%分位數(shù)作為 的估計(jì)值，第二種方法是將使得后驗(yàn)密度最大的θ值作為其估計(jì)值。

(5)洪水頻率設(shè)計(jì)值計(jì)算。假設(shè)已知洪水觀測值x的情況下，θ的后驗(yàn)分布可以用f(θ|x)來表示，設(shè)定 表示需要計(jì)算的一定重現(xiàn)期下的洪水設(shè)計(jì)值，那么在已知x的條件下，z的密度函數(shù)可以表示為：

(3)

由式(3)得：

(4)

Pr(Z≤z|x)可以理解為，已知洪水資料條件下，某一重現(xiàn)期洪水頻率設(shè)計(jì)值z的分布函數(shù)。

如果假設(shè)洪水重現(xiàn)期為m年一遇，即Pr(Z≤z|x)=1-1/m，通過求解方程即可求得設(shè)計(jì)值。由于在輸入的時(shí)候使用了樣本參數(shù)的先驗(yàn)分布以及似然函數(shù)，因此貝葉斯方法可以將模型參數(shù)的不確定性考慮進(jìn)去，得出的設(shè)計(jì)值也是概率分布的形式，體現(xiàn)了由于參數(shù)不確定性引起的頻率設(shè)計(jì)的不確定性。

在求解方程的過程中，由于后驗(yàn)概率密度較為復(fù)雜，需要借助特定的抽樣模擬方法估計(jì)式(4)后驗(yàn)分布的估計(jì)值。

2 MCMC模擬

MCMC方法[6]的基本流程如下。

(1)根據(jù)參數(shù)的后驗(yàn)分布，通過MCMC方法產(chǎn)生服從式(2)要求的一階馬爾科夫鏈的模擬序列θ0,θ1,θ2,…，其中θ0為任意初始值，θi+1只依賴于當(dāng)前的θi，與之前的序列θ0,θ1,…,θi-1無關(guān)，即θi+1由條件分布q(·|θi)產(chǎn)生。

(2)設(shè)定轉(zhuǎn)移核。為了保證馬爾科夫鏈?zhǔn)菚r(shí)間齊次的，即任意參數(shù)經(jīng)過一步迭代以后的邊緣分布與i無關(guān)，這里定義轉(zhuǎn)移核為p(θ,θ′)=q(θi+1=θ′|θi=θ)。并使得后驗(yàn)分布f(θ|x)為平穩(wěn)分布，即要求：

?θ′∈Θ

(5)

(3)通過MCMC方法產(chǎn)生隨機(jī)樣本。設(shè)定任意一個(gè)初始的θ0，由(1)和(2)中定義的馬爾克夫鏈進(jìn)行抽樣模擬，根據(jù)訓(xùn)練步長產(chǎn)生隨機(jī)序列θ1,θ2,…,θn。文獻(xiàn)[7]指出，由于初始設(shè)定θ0的邊緣分布并不一定是f(θ|x)，需要一定的訓(xùn)練時(shí)期使得序列達(dá)到平穩(wěn)，因此在利用MCMC模擬的時(shí)候，通常去掉非平穩(wěn)序列的θ1,θ2,…,θk，用剩下的θk+1,θk+2,…,θn作為后驗(yàn)分布f(θ|x)的抽樣即可。

(4)利用(3)中產(chǎn)生的序列進(jìn)行蒙特卡洛積分，本文選擇模擬樣本的50%分位數(shù)作為后驗(yàn)分布的均值，而后驗(yàn)分布密度由模擬序列的頻率直方圖生成。

在使用MCMC方法進(jìn)行模擬的時(shí)候，馬爾科夫鏈的穩(wěn)定性依賴于轉(zhuǎn)移核的選取。目前，常用于轉(zhuǎn)移核構(gòu)造的方法有兩種，分別是Metropolis-Hastings算法[7]和Gibbs抽樣[8]。Metropolis-Hastings算法(簡稱M-H算法)是最早提出是在1953年，經(jīng)過后人改進(jìn)成為目前應(yīng)用較為廣泛的MCMC方法。

Metropolis-Hastings算法基本步驟如下。

①確定參數(shù)的初值θ0以及轉(zhuǎn)移核q(·|θi)。

②進(jìn)行迭代，從q(·|θi)抽取一個(gè)推薦值θ*。

③計(jì)算相應(yīng)的接受概率αi。

(6)

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

(7)

式中：μ是均勻分布在一定區(qū)間內(nèi)的隨機(jī)數(shù)，其取值范圍在0～1之間。

(5)重復(fù)①～④的步驟n次，去除前面k個(gè)不穩(wěn)定序列后，認(rèn)為θk+1,θk+2,…為滿足抽樣要求的平穩(wěn)序列，根據(jù)這些平穩(wěn)序列進(jìn)行后驗(yàn)分布的各種統(tǒng)計(jì)推斷(后驗(yàn)均值及各個(gè)分位數(shù)等)。從極值理論角度來講，通過式(7)給出的概率進(jìn)行優(yōu)選可以使得序列達(dá)到平穩(wěn)狀態(tài)，并且保證模擬序列具備所求的邊緣分布的統(tǒng)計(jì)特征。

3 GEV模型分布參數(shù)與洪水的估計(jì)

作為搭建極值變量的一個(gè)重要分布，GEV分布在水文、氣象、環(huán)境、保險(xiǎn)和金融等領(lǐng)域得到廣泛的應(yīng)用[9, 10]。對于GEV分布，可以采用極大似然估計(jì)法和矩估計(jì)等對其參數(shù)進(jìn)行估計(jì)，其中極大似然估計(jì)方法是目前采用較多的方法[11]。盡管極大似然估計(jì)法可以計(jì)算極值變量分位數(shù)的置信區(qū)間，以此表達(dá)相應(yīng)估計(jì)值的不確定性，然而，極大似然估計(jì)方法要求待估計(jì)的樣本具有漸近正態(tài)性的特征，對于樣本長度要求較高，當(dāng)樣本系列較小時(shí)，極大似然估計(jì)方法不能表達(dá)分位數(shù)估計(jì)的中所帶來不確定性。Martins和Stedinger[12]指出，序列長度對極大似然估計(jì)有一定的影響，較小的樣本會造成極大似然估計(jì)方法估計(jì)出的形狀參數(shù)偏小，從而影響對于總體分位數(shù)估計(jì)的準(zhǔn)確性。

GEV分布函數(shù)可以表示為：

(8)

式中：μ、σ、ξ分別表示GEV分布模型的位置參數(shù)(Location parameter)、尺度參數(shù)(Scale parameter)和形狀參數(shù)(Shape parameter)，且同時(shí)滿足條件為μ∈R，σ>0，ξ∈R，1+ξ(x-μ)σ>0。

分別選取6個(gè)站點(diǎn)的年最大洪峰流量系列建立GEV分布模型。為方便計(jì)算和表達(dá)，洪峰流量系列的單位統(tǒng)一為103m3/s，以下類同。令φ=logσ，以保證σ>0。先驗(yàn)密度參數(shù)確定如下：選擇先驗(yàn)密度函數(shù)為π(μ,φ,ξ)=πμ(μ)πφ(φ)πξ(ξ)。其中πμ(·)、πφ(·)和πξ(·)為正態(tài)密度函數(shù)，其統(tǒng)計(jì)特征是均值為零，方差分別為υμ、υφ和υξ，同時(shí)保證參數(shù)μ、φ和ξ是相互獨(dú)立的。根據(jù)文獻(xiàn)[11]的推薦，初步設(shè)定υμ=υφ=104，υξ=100，從而幫助密度函數(shù)變得更加平滑。

參考上文中M-H算法的步驟[13]，利用概率αi對向量(μ,σ,ξ)的每個(gè)分量進(jìn)行篩選，同時(shí)定義每個(gè)分量在坐標(biāo)軸上隨機(jī)移動的步長為：μ*=μ+εμ，φ*=φ+εφ，ξ*=ξ+εξ，其中，εμ、εφ、εξ是均值為0，方差分別為ωμ、ωφ、ωξ的正態(tài)變量。ωμ、ωφ、ωξ的數(shù)值根據(jù)試驗(yàn)確定。

4 大渡河流域典型站點(diǎn)分析

選取大渡河流域6個(gè)典型站點(diǎn)年最大洪峰系列進(jìn)行分析，分別為丹巴、大金、猴子巖、瀘定、瀑布溝和雙江口，6個(gè)站點(diǎn)覆蓋流域大部分地區(qū)，具有一定代表性，將各個(gè)站點(diǎn)年最大洪峰流量系列作為洪水頻率計(jì)算的樣本數(shù)據(jù)。

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

設(shè)定每個(gè)參數(shù)的訓(xùn)練次數(shù)為20 000次，通過MCMC模擬產(chǎn)生丹巴站GEV分布的3個(gè)參數(shù)序列如圖2所示。其中，位置參數(shù)和形狀參數(shù)為原序列，而尺度參數(shù)是經(jīng)過σi=eφi轉(zhuǎn)化以后得到的模擬序列。由于初始值的邊緣分布并不一定滿足后驗(yàn)分布要求，因此去掉前面不平穩(wěn)的序列之后，可以認(rèn)為將其余比較穩(wěn)定的模擬值當(dāng)做符合要求的后驗(yàn)分布的觀測值。

由于方差ωμ、ωφ、ωξ的設(shè)定具有一定的主觀性，為了進(jìn)一步明晰不同的ωμ、ωφ、ωξ對參數(shù)收斂速度的影響，設(shè)置了2組不同的數(shù)值進(jìn)行對比分析如圖2所示，在圖2(a)的MCMC模擬中，設(shè)置ωμ、ωφ、ωξ的取值分別為0.02、0.01、0.1，從圖2(a)中可以看出，經(jīng)過2500次迭代后序列趨于穩(wěn)定；而圖2(b)的MCMC模擬中，設(shè)置ωμ、ωφ、ωξ的取值分別為0.1、0.05、0.5，為圖2(a)的5倍，從圖2(b)中可以看出，經(jīng)過1000次迭代后序列趨于穩(wěn)定。通過對比可以看出，設(shè)置較大的方差，可以一定程度上加快MCMC序列的收斂速度，同時(shí)，對比兩個(gè)情景下參數(shù)的收斂值可以發(fā)現(xiàn)，兩個(gè)情景下的參數(shù)收斂值基本上是一致的，這說明選取ωμ、ωφ、ωξ的大小僅僅對于模型收斂速度產(chǎn)生影響，而并不影響對模型參數(shù)的估計(jì)。

按照上述的方法，選取ωμ=0.1、ωφ=0.05、ωξ=0.5的值得出其余5個(gè)站點(diǎn)的MCMC序列如圖3所示。

根據(jù)20 000組參數(shù)，得到每個(gè)站點(diǎn)GEV各個(gè)參數(shù)后驗(yàn)分布統(tǒng)計(jì)特征如表1所示。相比傳統(tǒng)的參數(shù)估計(jì)方法，貝葉斯分布不僅給出參數(shù)的估計(jì)值，同時(shí)也給出參數(shù)置信區(qū)間，通過置信區(qū)間表示參數(shù)估計(jì)的不確定性。

4.2 擬合優(yōu)度檢驗(yàn)

(1)圖像分析法。圖像分析法一種較為直觀的利用圖形來描述擬合的優(yōu)劣程度的方法。主要是通過點(diǎn)繪理論聯(lián)合概率值和經(jīng)驗(yàn)聯(lián)合概率值，如果得到的點(diǎn)距較均勻地分布在45°線附近，則說明建立的概率分布模型是合理的[14]。圖4給出了大渡河流域各個(gè)站點(diǎn)最大洪峰流量系列GEV模型經(jīng)驗(yàn)頻率與理論頻率擬合效果圖，從圖4可以看出，樣本系列基本都落在了經(jīng)驗(yàn)頻率與理論頻率所在的45°線上，說明通過MCMC序列估計(jì)出來的GEV分布模型與實(shí)際樣本系列擬合程度較好。

圖2 丹巴站Bayes參數(shù)估計(jì)的MCMC模擬圖Fig.2 Bayes MCMC simulation for parameters at Danba station

圖3 各站Bayes參數(shù)估計(jì)的MCMC模擬圖 Fig.3 Bayes MCMC simulation for parameters of GEV model

圖4 大渡河流域各個(gè)站點(diǎn)GEV模型經(jīng)驗(yàn)頻率與理論頻率擬合效果圖Fig.4 The fitting effect chart for empirical frequency and theoretical frequency of GEV model at different stations in Dadu River basin

站點(diǎn)參數(shù)置信下限期望值置信上限位置參數(shù)3．2253．3963．568丹巴尺度參數(shù)0．3910．4980．644形狀參數(shù)-0．0900．0570．282位置參數(shù)1．8482．0642．277大金尺度參數(shù)0．4700．6120．802形狀參數(shù)-0．262-0．0560．203位置參數(shù)2．8603．0343．218猴子巖尺度參數(shù)0．4350．5430．682形狀參數(shù)-0．1260．0120．212位置參數(shù)3．1073．3273．556瀘定尺度參數(shù)0．5440．6870．868形狀參數(shù)-0．157-0．0240．154位置參數(shù)4．3744．5634．787瀑布溝尺度參數(shù)0．4890．6310．829形狀參數(shù)-0．0210．1700．400位置參數(shù)2．0302．2202．419雙江口尺度參數(shù)0．4380．5670．727形狀參數(shù)-0．0910．0690．284

(2)K-S檢驗(yàn)及OSL、AIC準(zhǔn)則。

離差平方和最小準(zhǔn)則(OSL)：

(9)

OSL值越小，說明GEV模型擬合得越好。

AIC信息準(zhǔn)則：

(10)

AIC=nln(MSE)+2k

(11)

式中：Fi表示經(jīng)驗(yàn)頻率；Ci表示理論頻率；n為分布樣本的個(gè)數(shù)；k為模型參數(shù)的個(gè)數(shù)，在GEV模型中k=3。

與OSL相同，AIC值越小，說明GEV模型擬合得越好。

采用OSL和AIC、Kolmogrov-Smirnow(簡稱K-S)法[15]3種指標(biāo)，對傳統(tǒng)的兩種參數(shù)估計(jì)方法和貝葉斯MCMC方法得到模型擬合效果進(jìn)行了對比分析，評定結(jié)果如表2所示。3種參數(shù)估計(jì)方法均通過K-S方法的檢驗(yàn)(顯著水平5%)，說明3種方法均適用于估計(jì)GEV模型的參數(shù)。對比各個(gè)參數(shù)估計(jì)方法的OSL和AIC可以看出，對于丹巴、大金、猴子巖、瀑布溝和瀘定站點(diǎn)，由貝葉斯估計(jì)參數(shù)得出的OSL和AIC的值是最小的，也就是說線型擬合的效果最好，而對于雙江口來說，矩估計(jì)得出的OSL和AIC值最小，但是矩估計(jì)的OSL和AIC值與貝葉斯估計(jì)相差較小，可以認(rèn)為矩估計(jì)和貝葉斯估計(jì)擬合出效果近似相同，極大似然估計(jì)相對比其他兩種估計(jì)方法，擬合效果相對較差。與傳統(tǒng)的參數(shù)估計(jì)方法相比，貝葉斯估計(jì)方法計(jì)算的OSL和AIC值總體上更小，這說明GEV模型擬合的效果更好，貝葉斯估計(jì)方法具有一定的優(yōu)越性。

表2 各個(gè)站點(diǎn)擬合優(yōu)度檢驗(yàn)結(jié)果Tab.2 Results of goodness of fit test

4.3 洪水頻率不確定性分析

估計(jì)一定重現(xiàn)期下的洪水設(shè)計(jì)值是水文頻率分析的主要目的之一，廣義極值分布的p(0

xp=μ-σ[1-(-logp)-ξ]/ξ

(12)

當(dāng)(μi,σi,ξi)時(shí)，即Gumbel分布的分位數(shù)為：

xp=μ-σlog(-logp)

(13)

(14)

(15)

貝葉斯方法將向量(μi,σi,ξi)的每個(gè)模擬值代入式(16)就得到相應(yīng)的1/(1-p)年重現(xiàn)期的設(shè)計(jì)洪峰流量值的后驗(yàn)分布樣本及置信區(qū)間。

(16)

以丹巴站為例，根據(jù)式(16)得出丹巴站年最大洪峰流量的重現(xiàn)水平圖如圖5所示，可以看出隨著重現(xiàn)期的增加，設(shè)計(jì)洪峰流量的置信區(qū)間也在增大，洪水頻率的不確定性也在增加。圖6給出了丹巴站各個(gè)典型重現(xiàn)期下的設(shè)計(jì)洪峰流量的后驗(yàn)密度估計(jì)圖。根據(jù)后驗(yàn)密度估計(jì)，計(jì)算出流域各站點(diǎn)各典型重現(xiàn)期下的設(shè)計(jì)洪峰流量值及95%置信度下的置信區(qū)間，并與傳統(tǒng)Delta進(jìn)行比較，結(jié)果如表3所示。由表3可以看出，根據(jù)貝葉斯MCMC方法得到的各個(gè)站點(diǎn)各個(gè)重現(xiàn)期的設(shè)計(jì)值總是小于其置信下限和置信上限的平均值，而傳統(tǒng)Delta方法得出的各個(gè)站點(diǎn)各個(gè)重現(xiàn)期的設(shè)計(jì)值是等于其置信下限和置信上限的平均值的，實(shí)際中，由于大的洪水資料比較有限，在估計(jì)某個(gè)洪水頻率設(shè)計(jì)值上限時(shí)往往不確定性要大于下限的不確定性，由此可以看出根據(jù)貝葉斯MCMC方法得到的置信區(qū)間往往更接近實(shí)際的情況。同時(shí)，從表3可以看到，重現(xiàn)期越長，重現(xiàn)水平越大，置信區(qū)間的寬度越寬，代表不確定性越大，與Delta方法相比，貝葉斯MCMC方法得到的置信區(qū)間的寬度相對要小一些，說明該方法在評估洪水頻率不確定性時(shí)要比Delta方法得到的不確定性要小。綜上所述，貝葉斯MCMC方法在評估洪水頻率較傳統(tǒng)的Delta方法具有一定的優(yōu)越性。

圖5 丹巴站Bayes重現(xiàn)水平圖Fig.5 Return level diagram of GEV model at Danba station

5 結(jié) 語

以GEV模型為洪水頻率分析線型，應(yīng)用貝葉斯MCMC算法對大渡河流域洪水頻率分析不確定性進(jìn)行研究，研究結(jié)果表明，Metropolis-Hastings抽樣的貝葉斯MCMC方法是洪水頻率參數(shù)估計(jì)的一個(gè)有效方法，與傳統(tǒng)的極大似然估計(jì)和矩估計(jì)相比，貝葉斯MCMC估計(jì)方法的擬合效果略優(yōu)于另外兩種方法。貝葉斯MCMC方法通過利用參數(shù)的先驗(yàn)信息得到參數(shù)以及設(shè)計(jì)值的后驗(yàn)密度函數(shù)，通過參數(shù)和設(shè)計(jì)值概率分布的形式表達(dá)洪水頻率中的不確定性信息。同時(shí)，將貝葉斯MCMC方法得出的置信區(qū)間與傳統(tǒng)Delta方法得出的置信區(qū)間進(jìn)行比較，發(fā)現(xiàn)貝葉斯MCMC方法得到的置信區(qū)間的寬度相對要小 一些，置信上限和與估計(jì)值之間距離大于置信下限與估計(jì)值之間距離，這種不對等性與實(shí)際更加接近，更能精準(zhǔn)地估計(jì)洪水頻率的置信區(qū)間。

圖6 丹巴站各個(gè)典型重現(xiàn)期設(shè)計(jì)洪峰流量后驗(yàn)概率密度圖Fig.6 Posterior probability density map of flood peak for different return periods at Danba station

m3/s

對于洪水頻率不確定性研究，本文僅僅采用GEV一種線型進(jìn)行研究，沒有考慮線型選擇對于洪水頻率分析不確定性造成的影響，未來工作需要進(jìn)一步將國內(nèi)較為常見的P-Ⅲ、Lognormal等線型加入，綜合分析洪水頻率分析中各個(gè)環(huán)節(jié)的不確定性。

□

[1] 朱慧明，韓玉啟. 貝葉斯多元統(tǒng)計(jì)推斷理論[M]. 上海：科學(xué)出版社, 2006.

[2] Wood E F, Rodríguez Iturbe I. Bayesian inference and decision making for extreme hydrologic events[J]. Water Resources Research, 1975,11(4):533-542.

[3] Claps P, Laio F. Can continuous streamflow data support flood frequency analysis? An alternative to the partial duration series approach[J]. Water Resources Research, 2003,39(8):1 216-1 225.

[4] Fernandes W, Naghettini M, Loschi R. A Bayesian approach for estimating extreme flood probabilities with upper-bounded distribution functions[J]. Stochastic Environmental Research and Risk Assessment, 2010,24(8):1 127-1 143.

[5] 朱 嵩，毛根海，劉國華，等. 改進(jìn)的 MCMC 方法及其應(yīng)用[J]. 水利學(xué)報(bào)，2009,(8)：1 019-1 023.

[6] 梁忠民，戴 榮，李彬權(quán). 基于貝葉斯理論的水文不確定性分析研究進(jìn)展[J]. 水科學(xué)進(jìn)展，2010,21(2):274-281.

[7] 陳 平，徐若曦. Metropolis-Hastings 自適應(yīng)算法及其應(yīng)用[J]. 系統(tǒng)工程理論與實(shí)踐，2008,28(1):100-108.

[8] 張志華，屈 斐. 軟件可靠性模型的 Bayes 推斷及 Gibbs 算法[J]. 高校應(yīng)用數(shù)學(xué)學(xué)報(bào)(A 輯)，2005,20(4):401-408.

[9] Coles S G, Tawn J A. A Bayesian analysis of extreme rainfall data[J]. Applied Statistics, 1996:463-478.

[10] Walshaw D. Modelling extreme wind speeds in regions prone to hurricanes[J]. Journal of the Royal Statistical Society (Series C Applied Statistics), 2000,49(1):51-62.

[11] Huang W, Xu S, Nnaji S. Evaluation of GEV model for frequency analysis of annual maximum water levels in the coast of United States[J]. Ocean Engineering, 2008,35(11):1 132-1 147.

[12] Reis Jr D S, Stedinger J R. Bayesian MCMC flood frequency analysis with historical information[J]. Journal of Hydrology, 2005,313(1):97-116.

[13] 魯 帆，嚴(yán)登華. 基于廣義極值分布和 Metropolis-Hastings 抽樣算法的貝葉斯 MCMC 洪水頻率分析方法[J]. 水利學(xué)報(bào)，2013,44(8):942-949.

[14] 曾惠芳，朱慧明，李素芳，等. 基于 MH 算法的貝葉斯分位自回歸模型[J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版)，2010,(2):88-92.

[15] 張冬冬，魯 帆，嚴(yán)登華，等. 基于 Archimedean Copula 函數(shù)的洪水多要素聯(lián)合概率分布研究[J]. 中國農(nóng)村水利水電，2015,(1):68-74.