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

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

1 研究背景

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

因此，防洪風(fēng)險(xiǎn)分析顯得尤為重要，對(duì)水庫(kù)大壩防洪能力進(jìn)行評(píng)估能最大限度地降低其風(fēng)險(xiǎn)和災(zāi)害影響，具有重大意義。故本文針對(duì)水庫(kù)大壩泄流能力不足或遭遇超標(biāo)準(zhǔn)洪水所引起水庫(kù)漫壩的風(fēng)險(xiǎn)進(jìn)行研究。一般的風(fēng)險(xiǎn)分析方法并未考慮長(zhǎng)期服役下大壩防洪性能的衰退，也未能對(duì)不同服役環(huán)境下大壩實(shí)際性態(tài)的

差異給予定量的分析。據(jù)此，本文在時(shí)變效應(yīng)理論的基礎(chǔ)上，試圖構(gòu)建時(shí)變隨機(jī)變量量化的函數(shù)模型，并提出緩變型防洪風(fēng)險(xiǎn)分析模型，同時(shí)運(yùn)用蒙特卡羅法求解風(fēng)險(xiǎn)率，從而定量分析大壩防洪時(shí)變安全度。

2 時(shí)變風(fēng)險(xiǎn)率計(jì)算方法

2.1 蒙特卡羅法的基本原理

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

首先選用一種特定的方式構(gòu)造隨機(jī)數(shù)，在[0,1]上均勻分布的隨機(jī)數(shù)是最基礎(chǔ)的，其他分布的隨機(jī)數(shù)均可在[0,1]范圍內(nèi)的均勻分布基礎(chǔ)上經(jīng)變換獲得?；旌贤喾╗12]以速度快、周期長(zhǎng)、統(tǒng)計(jì)性質(zhì)好的特點(diǎn)而被廣泛運(yùn)用。

混合同余法的計(jì)算公式為：

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

(1)

式中：a、b、m都取非負(fù)整數(shù)，m≠0。

計(jì)算時(shí)，引入一個(gè)參數(shù)ki，令：

(2)

式中：符號(hào)int表示取整。

則：

xi+1=axi+b-mki

(3)

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

(4)

2.2 隨機(jī)變量的抽樣

(5)

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

(6)

依據(jù)大數(shù)定律計(jì)算其失效概率Pf為：

(7)

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

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

(8)

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

(9)

求解風(fēng)險(xiǎn)率的具體形式如下：

(1)依據(jù)實(shí)測(cè)資料定義功能函數(shù)中每個(gè)隨機(jī)變量以及其分布規(guī)律；

(2)產(chǎn)生0到1之間的均勻隨機(jī)數(shù)ri；

(3)由ri產(chǎn)生各個(gè)隨機(jī)變量的抽樣值；

(4)將隨機(jī)變量的抽樣值依次代入功能函數(shù)式，求解功能函數(shù)值z(mì)i；

(5)統(tǒng)計(jì)出現(xiàn)zi<0的次數(shù)M和總抽樣次數(shù)N；

3 緩變型防洪風(fēng)險(xiǎn)分析模型

當(dāng)壩前水深接近壩頂高程時(shí)，大壩防洪風(fēng)險(xiǎn)將顯著標(biāo)志，建立其極限狀態(tài)方程:

Z=H-D=0

(10)

由于大壩服役過(guò)程中內(nèi)外環(huán)境等因素的變化，致使各種隨機(jī)性因素發(fā)生緩慢變異，有學(xué)者將該時(shí)變效應(yīng)稱之為緩變型時(shí)變效應(yīng)[13-15]。大壩防洪風(fēng)險(xiǎn)的緩變特性隨時(shí)間逐步發(fā)生，或許是結(jié)構(gòu)上的變化，或許是非結(jié)構(gòu)上的。在防洪時(shí)變風(fēng)險(xiǎn)分析中，各種不確定性因素將發(fā)生緩慢變化，D(n)是壩頂高程一個(gè)隨機(jī)過(guò)程，故利用衰減函數(shù)模型來(lái)表述壩頂高程的沉降緩變歷程：

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

(11)

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

(12)

式中：T為大壩設(shè)計(jì)基準(zhǔn)期;k需要通過(guò)工程經(jīng)驗(yàn)或試驗(yàn)資料確定。

壩前最高水位變量Hi對(duì)極值Ⅰ型分布、對(duì)數(shù)正態(tài)分布和正態(tài)分布均不拒絕[16-17]，依據(jù)實(shí)測(cè)資料對(duì)各年的壩前最高水位Hi進(jìn)行排序(H1

(13)

隨機(jī)變量Hi的均值μH、標(biāo)準(zhǔn)差σH及變異系數(shù)δH可應(yīng)用矩法公式求解：

(14)

(15)

(16)

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

Dn=

(17)

根據(jù)顯著性水平(一般取95%)，查K-S檢驗(yàn)臨界值表得Dn,95%與Dn，通過(guò)對(duì)比來(lái)判辨分布概型的擬合度，從而獲得隨機(jī)變量Hi的最優(yōu)概率分布。假定H與D是相互獨(dú)立且服從正態(tài)分布的隨機(jī)變量，D的均值μD和標(biāo)準(zhǔn)差σD由以下公式確定：

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

(18)

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

(19)

隨機(jī)變量H的均值μH和標(biāo)準(zhǔn)差σH由荷載評(píng)估期TE內(nèi)荷載統(tǒng)計(jì)分析方法求解：

(20)

(21)

由此可知，大壩第n年的防洪風(fēng)險(xiǎn)率的計(jì)算公式如下：

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

(22)

若隨機(jī)變量H和D分別服從對(duì)數(shù)正態(tài)分布和正態(tài)分布，經(jīng)公式(4)產(chǎn)生[0,1]上的均勻隨機(jī)數(shù)ri，再運(yùn)用坐標(biāo)變換法獲得正態(tài)分布N(0,1)的兩個(gè)隨機(jī)數(shù)ui1，ui2：

(23)

(24)

由于隨機(jī)變量H服從對(duì)數(shù)正態(tài)分布，令H′=lnH，可知H′為正態(tài)分布，則有：

(25)

(26)

(27)

具體計(jì)算流程見(jiàn)圖1。

4 實(shí)例分析

4.1 工程概況及計(jì)算參數(shù)

金溪流域上池潭混凝土壩已正常運(yùn)行多年，選取該壩的典型壩段進(jìn)行防洪風(fēng)險(xiǎn)分析，典型壩段剖面示意圖見(jiàn)圖2。

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

圖1 計(jì)算流程圖

圖2 典型壩段剖面示意圖

選取1990年至2011年間的壩前實(shí)測(cè)年最高水位以及對(duì)應(yīng)的壩頂實(shí)測(cè)年沉陷最大值進(jìn)行分析，如表1所述。

表1 服役期實(shí)測(cè)資料

對(duì)上文所述的3種分布函數(shù)進(jìn)行K-S檢驗(yàn)，其檢驗(yàn)成果見(jiàn)表2。

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

表2 隨機(jī)變量分布檢驗(yàn)表

本工程中的設(shè)計(jì)基準(zhǔn)期為50 a(TN=50)，結(jié)構(gòu)荷載評(píng)估基準(zhǔn)期TE可采用類比法式求解：

(28)

隨機(jī)變量H在繼續(xù)使用期TM內(nèi)的統(tǒng)計(jì)參數(shù)μH、σH可依據(jù)公式(20)、(21)計(jì)算，其計(jì)算結(jié)果見(jiàn)表3。

由于缺乏壩頂沉降衰減模型的資料，故選用公式(12)模擬沉降衰減；根據(jù)已服役22 a內(nèi)壩頂沉降統(tǒng)計(jì)的的均值近似確定壩頂高程的衰減函數(shù)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 大壩防洪風(fēng)險(xiǎn)率評(píng)估

選用表3所得的荷載評(píng)估期統(tǒng)計(jì)參數(shù)，依據(jù)公式(20～22)來(lái)推求大壩繼續(xù)使用期內(nèi)的防洪風(fēng)險(xiǎn)率，并計(jì)算現(xiàn)行規(guī)范設(shè)防標(biāo)準(zhǔn)下的防洪風(fēng)險(xiǎn)率，其求解成果如表4及圖3、4所示。

表3 隨機(jī)變量H在繼續(xù)使用期內(nèi)的統(tǒng)計(jì)參數(shù)表

表4 大壩繼續(xù)使用期內(nèi)防洪風(fēng)險(xiǎn)率統(tǒng)計(jì)表

圖3 大壩繼續(xù)使用期內(nèi)的防洪風(fēng)險(xiǎn)率

圖4 大壩后續(xù)服役期防洪風(fēng)險(xiǎn)率與現(xiàn)行規(guī)范設(shè)防標(biāo)準(zhǔn)對(duì)比圖

由圖3和4可知，大壩在服役過(guò)程中其防洪性能伴隨時(shí)間不斷的衰退，防洪風(fēng)險(xiǎn)率逐步增大，算例中的大壩在后續(xù)服役期內(nèi)的防洪風(fēng)險(xiǎn)率完全符合現(xiàn)行規(guī)范設(shè)防標(biāo)準(zhǔn)，且存在一定的安全余度。為進(jìn)一步確保大壩生命周期內(nèi)的防洪安全，減小洪水漫頂風(fēng)險(xiǎn)率 ，建議在工程建設(shè)前期充分考慮泥沙作用效應(yīng)，并預(yù)設(shè)一定的淤積庫(kù)容；也可在后續(xù)服役期內(nèi)增加壩頂高程，以增大防洪庫(kù)容來(lái)提升大壩防洪能力。

5 結(jié) 論

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

(2)將時(shí)變效應(yīng)引入防洪的風(fēng)險(xiǎn)分析中，構(gòu)建了時(shí)變隨機(jī)變量量化的函數(shù)模型，據(jù)此提出了緩變型防洪風(fēng)險(xiǎn)分析模型。

(3)根據(jù)實(shí)際工程的實(shí)測(cè)資料序列，采用隨機(jī)風(fēng)險(xiǎn)率計(jì)算方法定量地分析大壩后續(xù)服役期內(nèi)的防洪風(fēng)險(xiǎn)率，評(píng)估大壩繼續(xù)使用期內(nèi)的防洪能力。同時(shí)，本文的風(fēng)險(xiǎn)分析方法亦可為其他水利工程防洪風(fēng)險(xiǎn)評(píng)估提供參考。