基于蒙特卡羅法的緩變型防洪風險分析

周方明， 王 鋒， 谷明晗， 蔡 磊

(1.宿州學院 資源與土木工程學院， 安徽 宿州 234000； 2.江西水利職業學院 水利工程系，江西 南昌 330013； 3.南昌大學 建筑工程學院， 江西 南昌 330031)

1 研究背景

洪澇災害是世界發生頻次最高、損失最大的自然災害，我國大部分地區均存在著不同程度的洪澇災害，對人民生命財產、國家經濟建設、地域生態環境構成了重大威脅[1-3]。世界壩工委員會[4]于1998年總結了水庫防洪安全設計的進展，主要是根據經驗和水庫的共性粗略分級制定標準列表備查，并未考慮到水庫的個性，如水庫大小、壩型、庫容以及垮壩的危害。外國學者Kwon等[5]采用了Monte-Carlo模擬隨機變量分布，用來研究大壩的水文風險概率；國內學者梅亞東等[6]認為水庫大壩在長期服役過程中，影響防洪風險的諸多因素(水力、水文、設計)伴隨時間歷程發生隨機不確定性的變化，致使大壩的防洪能力呈現明顯的時變特性[7-9]。時變特性具體呈現為不確定性洪水的荷載效應與結構的抗力作用是隨機變化的過程。故防洪建筑物的失事因素與時變效應存在著協同作用，導致了諸多破壞形式和失效模式，如洪水漫頂、邊坡失穩、滲透破壞等；其中洪水漫頂是最常發生的失效模式。

因此，防洪風險分析顯得尤為重要，對水庫大壩防洪能力進行評估能最大限度地降低其風險和災害影響，具有重大意義。故本文針對水庫大壩泄流能力不足或遭遇超標準洪水所引起水庫漫壩的風險進行研究。一般的風險分析方法并未考慮長期服役下大壩防洪性能的衰退，也未能對不同服役環境下大壩實際性態的

差異給予定量的分析。據此，本文在時變效應理論的基礎上，試圖構建時變隨機變量量化的函數模型，并提出緩變型防洪風險分析模型，同時運用蒙特卡羅法求解風險率，從而定量分析大壩防洪時變安全度。

2 時變風險率計算方法

2.1 蒙特卡羅法的基本原理

蒙特卡羅法通過構建一個概率模型或隨機過程，使其參數等于所求問題的解，然后通過隨機試驗或觀察等來計算所求參數的統計特征，由此得出問題的近似解，而解的精度可由估計值的標準差來表示，并且據此可得結果的置信度及置信區間，從而解決實際工程問題[10-11]。首先由多個相互獨立的隨機變量X1,X2,…,Xn構成功能函數為Z=g(X1,X2,…,Xn)，再分別對其抽樣獲得各隨機變量的分位值x1,x2,…,xn，接著對功能函數式zi=g(x1,x2,…,xn)進行求解，設抽樣次數為N，再計算每組對應的功能函數值zi，其中zi<0的次數為M，在足夠多的抽樣結束后，其失效概率為Pf=M/N。

首先選用一種特定的方式構造隨機數，在[0,1]上均勻分布的隨機數是最基礎的，其他分布的隨機數均可在[0,1]范圍內的均勻分布基礎上經變換獲得。混合同余法[12]以速度快、周期長、統計性質好的特點而被廣泛運用。

混合同余法的計算公式為：

xi+1=(axi+b)(modm)

(1)

式中：a、b、m都取非負整數，m≠0。

計算時，引入一個參數ki，令：

(2)

式中：符號int表示取整。

則：

xi+1=axi+b-mki

(3)

將xi+1除以m后，即得[0,1]上的均勻隨機數ri+1：

(4)

2.2 隨機變量的抽樣

(5)

隨后將抽樣值依次代入功能函數式中，求解成果唯有g(xi)>0和g(xi)≤0兩種情況，故可定義指標函數：

(6)

依據大數定律計算其失效概率Pf為：

(7)

式中：M為N次模擬計算中g(xi)≤0的總次數。

本文擬采納95%的置信度來確保蒙特卡洛法求解在容許誤差ε內：

(8)

由公式(8)可知，抽樣數目N越大，容許誤差ε隨之減小。若計算結果要滿足既定的精度要求，樣本數目N需滿足一定要求：

(9)

求解風險率的具體形式如下：

(1)依據實測資料定義功能函數中每個隨機變量以及其分布規律；

(2)產生0到1之間的均勻隨機數ri；

(3)由ri產生各個隨機變量的抽樣值；

(4)將隨機變量的抽樣值依次代入功能函數式，求解功能函數值zi；

(5)統計出現zi<0的次數M和總抽樣次數N；

3 緩變型防洪風險分析模型

當壩前水深接近壩頂高程時，大壩防洪風險將顯著標志，建立其極限狀態方程:

Z=H-D=0

(10)

由于大壩服役過程中內外環境等因素的變化，致使各種隨機性因素發生緩慢變異，有學者將該時變效應稱之為緩變型時變效應[13-15]。大壩防洪風險的緩變特性隨時間逐步發生，或許是結構上的變化，或許是非結構上的。在防洪時變風險分析中，各種不確定性因素將發生緩慢變化，D(n)是壩頂高程一個隨機過程，故利用衰減函數模型來表述壩頂高程的沉降緩變歷程：

D(n)=a(n,k)D0

(11)

式中：D(n)為服役n年以后隨機壩頂高程；D0為大壩的壩頂高程初設值；a(n,k)是衰減函數;k為衰減系數，相應的衰減模型用指數函數來表示：

(12)

式中：T為大壩設計基準期;k需要通過工程經驗或試驗資料確定。

壩前最高水位變量Hi對極值Ⅰ型分布、對數正態分布和正態分布均不拒絕[16-17]，依據實測資料對各年的壩前最高水位Hi進行排序(H1

(13)

隨機變量Hi的均值μH、標準差σH及變異系數δH可應用矩法公式求解：

(14)

(15)

(16)

依照公式(13)求解經驗分布F(Hk)，并采用K-S檢驗法[18]依次檢驗相應的經驗分布F(Hk)，建立統計量：

Dn=

(17)

根據顯著性水平(一般取95%)，查K-S檢驗臨界值表得Dn,95%與Dn，通過對比來判辨分布概型的擬合度，從而獲得隨機變量Hi的最優概率分布。假定H與D是相互獨立且服從正態分布的隨機變量，D的均值μD和標準差σD由以下公式確定：

μD(n)=μD0·a(n,k)

(18)

σD(n)=σD0·a(n,k)

(19)

隨機變量H的均值μH和標準差σH由荷載評估期TE內荷載統計分析方法求解：

(20)

(21)

由此可知，大壩第n年的防洪風險率的計算公式如下：

Pf(n)=Φ(-βn)=

(22)

若隨機變量H和D分別服從對數正態分布和正態分布，經公式(4)產生[0,1]上的均勻隨機數ri，再運用坐標變換法獲得正態分布N(0,1)的兩個隨機數ui1，ui2：

(23)

(24)

由于隨機變量H服從對數正態分布，令H′=lnH，可知H′為正態分布，則有：

(25)

(26)

(27)

具體計算流程見圖1。

4 實例分析

4.1 工程概況及計算參數

金溪流域上池潭混凝土壩已正常運行多年，選取該壩的典型壩段進行防洪風險分析，典型壩段剖面示意圖見圖2。

工程主要建筑物為二級永久建筑物，壩頂高程為384.5 m，巖基建基面高程314.5 m，正常蓄水位382.0 m，設計洪水位(P=1%)為382.0 m，校核洪水位(P=0.1%)為384.0 m，死水位354.0 m，為年調節水庫。

圖1 計算流程圖

圖2 典型壩段剖面示意圖

選取1990年至2011年間的壩前實測年最高水位以及對應的壩頂實測年沉陷最大值進行分析，如表1所述。

表1 服役期實測資料

對上文所述的3種分布函數進行K-S檢驗，其檢驗成果見表2。

查K-S檢驗臨界值表獲得允許概率差值D22.95%=0.2807，與表2的計算成果對比分析，3種分布的最大概率差值Dn均小于D22.95%，表明隨機變量H和D均服從這3種分布的假定，且正態分布的Dn值小于對數正態分布和極值I型分布的Dn，說明實測序列值與正態分布假定之間差別最小，隨機變量均為正態分布時的擬合效果最好。

表2 隨機變量分布檢驗表

本工程中的設計基準期為50 a(TN=50)，結構荷載評估基準期TE可采用類比法式求解：

(28)

隨機變量H在繼續使用期TM內的統計參數μH、σH可依據公式(20)、(21)計算，其計算結果見表3。

由于缺乏壩頂沉降衰減模型的資料，故選用公式(12)模擬沉降衰減；根據已服役22 a內壩頂沉降統計的的均值近似確定壩頂高程的衰減函數a(n,k)。

已知μD0=384.5 m，μD(22)=384.495 m，σD(22)=7.43×10-4m，則由公式(18)、(19)可得：

μD(n)=μD0·a(n,k)

(29)

σD(n)=σD0·a(n,k)

(30)

4.2 大壩防洪風險率評估

選用表3所得的荷載評估期統計參數，依據公式(20～22)來推求大壩繼續使用期內的防洪風險率，并計算現行規范設防標準下的防洪風險率，其求解成果如表4及圖3、4所示。

表3 隨機變量H在繼續使用期內的統計參數表

表4 大壩繼續使用期內防洪風險率統計表

圖3 大壩繼續使用期內的防洪風險率

圖4 大壩后續服役期防洪風險率與現行規范設防標準對比圖

由圖3和4可知，大壩在服役過程中其防洪性能伴隨時間不斷的衰退，防洪風險率逐步增大，算例中的大壩在后續服役期內的防洪風險率完全符合現行規范設防標準，且存在一定的安全余度。為進一步確保大壩生命周期內的防洪安全，減小洪水漫頂風險率 ，建議在工程建設前期充分考慮泥沙作用效應，并預設一定的淤積庫容；也可在后續服役期內增加壩頂高程，以增大防洪庫容來提升大壩防洪能力。

5 結 論

(1)本文針對制約大壩防洪能力的不確定性因素，探究了適宜解決隨機問題的蒙特卡羅法(Monte Carlo Method)，并論述了其基本原理和求解過程。

(2)將時變效應引入防洪的風險分析中，構建了時變隨機變量量化的函數模型，據此提出了緩變型防洪風險分析模型。

(3)根據實際工程的實測資料序列，采用隨機風險率計算方法定量地分析大壩后續服役期內的防洪風險率，評估大壩繼續使用期內的防洪能力。同時，本文的風險分析方法亦可為其他水利工程防洪風險評估提供參考。