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

張冬冬，魯 帆，劉少華

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

由式(3)得：

(4)

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

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

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

2 MCMC模擬

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

(1)根據(jù)參數(shù)的后驗分布，通過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ù)經(jīng)過一步迭代以后的邊緣分布與i無關(guān)，這里定義轉(zhuǎn)移核為p(θ,θ′)=q(θi+1=θ′|θi=θ)。并使得后驗分布f(θ|x)為平穩(wěn)分布，即要求：

?θ′∈Θ

(5)

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

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

在使用MCMC方法進(jìn)行模擬的時候，馬爾科夫鏈的穩(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)抽取一個推薦值θ*。

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

(6)

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

(7)

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

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

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

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

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

(8)

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

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

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

4 大渡河流域典型站點分析

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

4.1 參數(shù)估計

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

由于方差ωμ、ωφ、ωξ的設(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ù)的收斂值可以發(fā)現(xiàn)，兩個情景下的參數(shù)收斂值基本上是一致的，這說明選取ωμ、ωφ、ωξ的大小僅僅對于模型收斂速度產(chǎn)生影響，而并不影響對模型參數(shù)的估計。

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

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

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

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

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

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

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

站點參數(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檢驗及OSL、AIC準(zhǔn)則。

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

(9)

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

AIC信息準(zhǔn)則：

(10)

AIC=nln(MSE)+2k

(11)

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

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

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

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

4.3 洪水頻率不確定性分析

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

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

(12)

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

xp=μ-σlog(-logp)

(13)

(14)

(15)

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

(16)

以丹巴站為例，根據(jù)式(16)得出丹巴站年最大洪峰流量的重現(xiàn)水平圖如圖5所示，可以看出隨著重現(xiàn)期的增加，設(shè)計洪峰流量的置信區(qū)間也在增大，洪水頻率的不確定性也在增加。圖6給出了丹巴站各個典型重現(xiàn)期下的設(shè)計洪峰流量的后驗密度估計圖。根據(jù)后驗密度估計，計算出流域各站點各典型重現(xiàn)期下的設(shè)計洪峰流量值及95%置信度下的置信區(qū)間，并與傳統(tǒng)Delta進(jìn)行比較，結(jié)果如表3所示。由表3可以看出，根據(jù)貝葉斯MCMC方法得到的各個站點各個重現(xiàn)期的設(shè)計值總是小于其置信下限和置信上限的平均值，而傳統(tǒng)Delta方法得出的各個站點各個重現(xiàn)期的設(shè)計值是等于其置信下限和置信上限的平均值的，實際中，由于大的洪水資料比較有限，在估計某個洪水頻率設(shè)計值上限時往往不確定性要大于下限的不確定性，由此可以看出根據(jù)貝葉斯MCMC方法得到的置信區(qū)間往往更接近實際的情況。同時，從表3可以看到，重現(xiàn)期越長，重現(xiàn)水平越大，置信區(qū)間的寬度越寬，代表不確定性越大，與Delta方法相比，貝葉斯MCMC方法得到的置信區(qū)間的寬度相對要小一些，說明該方法在評估洪水頻率不確定性時要比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ù)估計的一個有效方法，與傳統(tǒng)的極大似然估計和矩估計相比，貝葉斯MCMC估計方法的擬合效果略優(yōu)于另外兩種方法。貝葉斯MCMC方法通過利用參數(shù)的先驗信息得到參數(shù)以及設(shè)計值的后驗密度函數(shù)，通過參數(shù)和設(shè)計值概率分布的形式表達(dá)洪水頻率中的不確定性信息。同時，將貝葉斯MCMC方法得出的置信區(qū)間與傳統(tǒng)Delta方法得出的置信區(qū)間進(jìn)行比較，發(fā)現(xiàn)貝葉斯MCMC方法得到的置信區(qū)間的寬度相對要小 一些，置信上限和與估計值之間距離大于置信下限與估計值之間距離，這種不對等性與實際更加接近，更能精準(zhǔn)地估計洪水頻率的置信區(qū)間。

圖6 丹巴站各個典型重現(xiàn)期設(shè)計洪峰流量后驗概率密度圖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等線型加入，綜合分析洪水頻率分析中各個環(huán)節(jié)的不確定性。

□

[1] 朱慧明，韓玉啟. 貝葉斯多元統(tǒng)計推斷理論[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é)報，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)工程理論與實踐，2008,28(1):100-108.

[8] 張志華，屈 斐. 軟件可靠性模型的 Bayes 推斷及 Gibbs 算法[J]. 高校應(yīng)用數(shù)學(xué)學(xué)報(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é)報，2013,44(8):942-949.

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

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