基于貝葉斯參數(shù)更新的高土石壩壩頂開裂風(fēng)險(xiǎn)動(dòng)態(tài)評(píng)估與預(yù)警

黃會(huì)寶，熊敏，高志良，吳震宇，李俊儒，尹川

（1.四川大學(xué)水利水電學(xué)院，四川成都 610065；2.四川大學(xué)水力學(xué)及山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川成都 610065；3.國(guó)能大渡河流域水電開發(fā)有限公司，四川成都 610041）

0 引言

土石壩對(duì)地形地質(zhì)條件適應(yīng)性強(qiáng)，可以充分利用當(dāng)?shù)夭牧?，造價(jià)較低，是水電開發(fā)中最具應(yīng)用前景的壩型之一［1，2］。在我國(guó)已建土石壩中，潰壩事故主要集中在中低高度大壩，高土石壩主要表現(xiàn)為壩頂開裂、壩體與壩肩（基）滲漏、壩坡滑坡塌坑等非潰壩故障模式，根據(jù)相關(guān)文獻(xiàn)統(tǒng)計(jì)，壩頂開裂約占高心墻堆石壩常見故障模式的40%以上［3-5］。

眾所周知，壩頂開裂與土石壩的變形特性密切相關(guān)，特別是變形傾度指標(biāo)［6-8］。通過大壩原觀監(jiān)測(cè)數(shù)據(jù)，結(jié)合參數(shù)反演算法和數(shù)值模擬方法可以對(duì)高心墻堆石壩運(yùn)行形態(tài)進(jìn)行準(zhǔn)確預(yù)測(cè)［9，10］。當(dāng)基于數(shù)值模擬結(jié)果計(jì)算出的壩頂變形傾度大于某一閾值時(shí)，壩頂存在開裂可能性，但是并不意味著壩頂會(huì)立即開裂［11，12］，若壩頂變形傾度超過閾值的異常狀態(tài)持續(xù)時(shí)間較長(zhǎng)，壩頂開裂的可能性就會(huì)隨著異常狀態(tài)的持續(xù)時(shí)間而增大。由于復(fù)雜的物理機(jī)制和隨機(jī)因素，壩頂開裂前的異常狀態(tài)持續(xù)時(shí)間是不確定的，可以看作是一個(gè)隨機(jī)變量，并可以基于此對(duì)異常狀態(tài)期間發(fā)生壩頂開裂的風(fēng)險(xiǎn)進(jìn)行評(píng)估。

由于缺乏相關(guān)的歷史資料，特別是對(duì)于某一具體的土石壩，無法通過經(jīng)典的頻率統(tǒng)計(jì)分析得到壩頂開裂前異常狀態(tài)持續(xù)時(shí)間的概率分布參數(shù)（PDP）。根據(jù)異常狀態(tài)開始到發(fā)生峰值之間的觀測(cè)時(shí)間間隔，可以利用貝葉斯理論估計(jì)PDP 的概率分布［13］。此外，PDP 的概率分布可以隨著觀察到的證據(jù)的積累過程而動(dòng)態(tài)更新。因此，在貝葉斯框架下，土石壩的壩頂開裂風(fēng)險(xiǎn)可以根據(jù)PDP 的逐步更新概率分布進(jìn)行動(dòng)態(tài)評(píng)估。此外，還可以根據(jù)概率估計(jì)的開裂時(shí)間范圍對(duì)壩頂開裂進(jìn)行預(yù)警。

本文研究的主要問題是：①如何用壩頂沉降異常狀態(tài)持續(xù)時(shí)間來評(píng)價(jià)土石壩壩頂開裂風(fēng)險(xiǎn)；②如何估計(jì)給定概率區(qū)間對(duì)應(yīng)的壩頂開裂次數(shù)以進(jìn)行壩頂開裂預(yù)警；③如何更新壩頂開裂風(fēng)險(xiǎn)公式中參數(shù)λ的概率分布。

1 壩頂開裂風(fēng)險(xiǎn)動(dòng)態(tài)評(píng)估與預(yù)警方法

1.1 總體思路與方法架構(gòu)

首先確定高心墻堆石壩壩頂異常狀態(tài)判定準(zhǔn)則。如引言所述，壩頂開裂與壩頂變形傾度γS密切相關(guān)。當(dāng)γS超過某一閾值時(shí)，壩頂易發(fā)生開裂故障。變形傾度的定義是［14］：

式中：γS為變形傾度；ΔS為同一高程上兩點(diǎn)的沉降差；ΔL為兩點(diǎn)之間的水平距離。

采用數(shù)值模擬計(jì)算壩頂裂縫兩側(cè)變形傾度，為使計(jì)算結(jié)果可以準(zhǔn)確描述壩體變形，基于大壩原觀變形監(jiān)測(cè)數(shù)據(jù)對(duì)數(shù)值模型的材料參數(shù)進(jìn)行反演分析。壩頂開裂前的異常變形狀態(tài)定義為γS＞γC，γC為異常變形的閾值（即臨界傾度）。有學(xué)者根據(jù)土梁撓曲試驗(yàn)，提出γC的取值范圍為0.5%～0.65%［15］；或根據(jù)部分土石壩開裂前變形傾度的統(tǒng)計(jì)數(shù)據(jù)，提出γC可采用0.6%［16］、0.7%［17］或1%［18］等值。從偏保守角度考慮，本文取γC=0.5%。

壩頂開裂前異常狀態(tài)的持續(xù)時(shí)間TAS是不確定的，因此本文將TAS視為隨機(jī)變量，并假定TAS服從指數(shù)分布，故而壩頂開裂的概率可以表示為異常狀態(tài)持續(xù)時(shí)間TAS的函數(shù)。一旦出現(xiàn)異常狀態(tài)（如γS＞γC），就可以用TAS來估計(jì)壩頂開裂的概率。此外，還可以推導(dǎo)出與給定的開裂概率相對(duì)應(yīng)的開裂時(shí)間。由于不具備與壩頂開裂概率相對(duì)應(yīng)的先驗(yàn)知識(shí)，因此通過對(duì)觀察到的開裂次數(shù)的計(jì)算概率的統(tǒng)計(jì)分析來選擇和更新估算開裂時(shí)間的概率。

壩頂開裂概率計(jì)算中的參數(shù)是不確定的。如何估計(jì)參數(shù)的概率分布是該方法應(yīng)用中的一個(gè)主要問題。應(yīng)用貝葉斯理論，根據(jù)觀測(cè)到的從出現(xiàn)異常狀態(tài)到出現(xiàn)壩頂開裂的時(shí)間間隔，估計(jì)和更新參數(shù)的概率分布。

上述方法的總體流程如圖1所示。

圖1 壩頂開裂動(dòng)態(tài)風(fēng)險(xiǎn)評(píng)價(jià)與預(yù)警流程圖Fig.1 Flowchart of the dynamic risk evaluation and early warning of crest cracking

1.2 壩頂開裂概率與預(yù)警

如前所述，壩頂開裂前異常狀態(tài)的持續(xù)時(shí)間TAS是不確定的，其不確定性可以用概率密度函數(shù)（PDF）定量描述。指數(shù)分布可以用來描述服從泊松分布的獨(dú)立隨機(jī)事件的時(shí)間間隔的概率分布［19］。因此，采用指數(shù)分布來表示壩頂開裂前異常狀態(tài)持續(xù)時(shí)間TAS的概率分布，其PDF如下：

式中：λ是TAS的概率密度函數(shù)參數(shù)，λ=E（TAS），即λ是壩頂開裂前異常狀態(tài)持續(xù)時(shí)間的期望值。

當(dāng)壩頂沉降的異常狀態(tài)（即γS＞γC）發(fā)生并持續(xù)時(shí)，在給定的時(shí)間內(nèi)出現(xiàn)壩頂開裂的概率為：

從公式（3）中可以看出，壩頂開裂風(fēng)險(xiǎn)評(píng)估的合理性與準(zhǔn)確性取決于參數(shù)λ的取值。由于缺乏高心墻壩壩頂開裂前異常狀態(tài)持續(xù)時(shí)間的歷史數(shù)據(jù)，使得參數(shù)λ值存在很大的認(rèn)識(shí)不確定性，特別是對(duì)于具體的大壩而言。在壩頂開裂風(fēng)險(xiǎn)評(píng)估中，應(yīng)考慮到這種認(rèn)識(shí)不確定性的影響。在貝葉斯框架下，認(rèn)知不確定性可以用主觀概率來量化，本研究應(yīng)用貝葉斯理論，根據(jù)觀測(cè)到的壩頂開裂前的異常狀態(tài)持續(xù)時(shí)間，對(duì)λ的主觀概率分布進(jìn)行估計(jì)和更新。

1.3 信息驅(qū)動(dòng)的風(fēng)險(xiǎn)模型參數(shù)更新

參數(shù)λ中的認(rèn)知不確定性源于分析人員對(duì)λ的統(tǒng)計(jì)信息的缺乏。這種認(rèn)知不確定性可以用貝葉斯框架中的主觀概率定量地表示。當(dāng)獲得與參數(shù)λ有關(guān)的信息時(shí)，可以根據(jù)貝葉斯定理估計(jì)和更新λ的概率分布如下：

式中：πi(λ|Ei)是在獲得第i個(gè)信息的情況下λ的后驗(yàn)概率分布；πi-1(λ)是λ在獲得第i個(gè)信息之前的先驗(yàn)概率分布；L(Ei|λ)是似然函數(shù)（即給定λ值下觀測(cè)信息的條件概率）；Θ 是λ的定義域。

如何建立似然函數(shù)是用公式（5）更新λ概率分布的關(guān)鍵步驟，在土石壩的長(zhǎng)期運(yùn)行過程中，可以觀察到與λ估算直接相關(guān)的兩類信息。根據(jù)這兩類信息，似然函數(shù)可構(gòu)造如下：

（2）第二類信息：從異常狀態(tài)出現(xiàn)到異常狀態(tài)結(jié)束，未出現(xiàn)壩頂開裂的持續(xù)時(shí)間，相應(yīng)的似然函數(shù)表述為：

根據(jù)新觀測(cè)信息的類型，將相應(yīng)的似然函數(shù)和λ的先驗(yàn)概率分布代入公式（5），即可更新λ的概率分布。當(dāng)?shù)谝淮喂烙?jì)λ的概率分布時(shí)，沒有任何關(guān)于λ的先驗(yàn)信息，此時(shí)采用均勻分布來表示π0(λ)：

式中：λmin和λmax分別是λ的可能最小值和最大值。

從λ的后驗(yàn)概率分布中抽樣，將λ的樣本代入公式（3）和（4），可以計(jì)算出壩頂開裂概率和開裂時(shí)間的均值和置信區(qū)間。

2 工程案例

2.1 工程概況

某心墻堆石壩最大壩高為186 m，壩頂寬度為14 m，主要防滲結(jié)構(gòu)為礫石土心墻及壩基防滲墻，心墻上、下游側(cè)與壩殼堆石間設(shè)反濾層及過渡層，其典型橫剖面見圖2，其中內(nèi)部水平變形測(cè)點(diǎn)（CH17～CH31）與沉降位移測(cè)點(diǎn)（SG17～SG31）同點(diǎn)布置。

圖2 某心墻堆石壩典型橫剖面（單位：m）Fig.2 Typical cross section of the core rockfill dam

大壩始建于2004 年3月，2009 年9 月竣工，該水庫的第一次蓄水始于2009 年11月，水庫水位于2010 年11 月上升到正常水位（850.00 m），此后水庫水位在死水位（790.00 m）和正常水位之間周期性波動(dòng)。大壩填筑過程和庫水位變化過程如圖3所示。

圖3 大壩填筑過程及庫水位變化過程Fig.3 Dam filling process and fluctuation of reservoir water level

2.2 大壩壩頂開裂情況

該大壩在首次達(dá)到正常蓄水位后，于2010 年12 月21 日首次在壩頂壩軸線位置發(fā)現(xiàn)縱向裂縫，采用瀝青封堵進(jìn)行處理，此后歷年高水位運(yùn)行期間均出現(xiàn)壩頂開裂現(xiàn)象。壩頂首次發(fā)生裂縫后進(jìn)行了挖探坑檢查，查明裂縫深度約為1.5～2 m，長(zhǎng)度約230 m。由于壩頂與礫石土心墻頂部的垂直距離為2 m，因此裂縫未延伸至礫石土心墻，此外，心墻滲壓監(jiān)測(cè)數(shù)據(jù)也沒有出現(xiàn)異常。雖然壩頂裂縫并沒有直接危害礫石土心墻的安全，但由于降雨入滲可能侵蝕大壩內(nèi)部材料，需要對(duì)裂縫及時(shí)進(jìn)行封堵，因此，在大壩運(yùn)行過程中有必要對(duì)壩頂開裂風(fēng)險(xiǎn)進(jìn)行評(píng)估和預(yù)警，提前做好壩頂開裂后的修復(fù)準(zhǔn)備。

3 案例分析結(jié)果

3.1 壩頂時(shí)序變形傾度及異常狀態(tài)識(shí)別

采用本文提出的方法對(duì)多次觀測(cè)到壩頂開裂的大壩典型橫剖面進(jìn)行分析，圖4 為該剖面的有限元模型。采用鄧肯-張E-B 模型計(jì)算大壩瞬時(shí)變形，采用七參數(shù)模型計(jì)算大壩蠕變變形［20］，采用三參數(shù)模型計(jì)算大壩首次蓄水產(chǎn)生的濕化變形［21］。實(shí)測(cè)水平及沉降位移與計(jì)算值對(duì)比如圖5所示。

圖4 大壩典型橫剖面的有限元模型Fig.4 Finite element model of typical cross section of the dam

圖5 實(shí)測(cè)水平及沉降位移與計(jì)算值的對(duì)比Fig.5 Comparison of the computed and measured horizontal and settlement displacements

在本研究中壩頂裂縫處變形傾度是通過有限元計(jì)算裂縫兩側(cè)節(jié)點(diǎn)沉降位移后代入式（1）得到的，為保證壩頂變形梯度的計(jì)算精度，基于SG18、SG25、CH19和TP14四個(gè)測(cè)點(diǎn)的監(jiān)測(cè)數(shù)據(jù)進(jìn)行參數(shù)反演，采用響應(yīng)面法代替耗時(shí)的有限元計(jì)算［22］，采用粒子群優(yōu)化算法進(jìn)行目標(biāo)函數(shù)尋優(yōu)［23］。圖5為各測(cè)點(diǎn)水平及沉降位移的實(shí)測(cè)值和反演計(jì)算值的對(duì)比，可以看出兩者的變化規(guī)律一致，量值接近，運(yùn)行期和實(shí)測(cè)水平及沉降位移與計(jì)算值的相對(duì)誤差在6%以內(nèi)，表明本文數(shù)值模型可以比較準(zhǔn)確地預(yù)測(cè)大壩運(yùn)行過程中的變形狀態(tài)。

計(jì)算得到的壩頂順河向沉降變形傾度過程線如圖6 所示，圖6 中還標(biāo)出了歷次壩頂開裂時(shí)間。由圖6 可見，壩頂順河向沉降變形傾度與庫水位具有明顯正相關(guān)關(guān)系，庫水位上升時(shí)變形傾度相應(yīng)增大，最大變形傾度為0.98%。由圖6 可見，壩頂開裂都發(fā)生在變形傾度較大時(shí)，此外，壩頂出現(xiàn)異常變形狀態(tài)到壩頂開裂之間的時(shí)間間隔具有明顯的隨機(jī)性。

圖6 壩頂變形傾度過程線及開裂時(shí)刻Fig.6 The time series of crest settlement inclination and the cracking time

3.2 開裂風(fēng)險(xiǎn)模型參數(shù)λ概率分布的貝葉斯估計(jì)

由于壩頂開裂前的異常狀態(tài)持續(xù)時(shí)間TAS是不確定的，本文假定其為服從指數(shù)分布的隨機(jī)變量。由于其參數(shù)λ也是不確定的，因此，根據(jù)TAS的觀測(cè)值，應(yīng)用2.3 節(jié)所述的貝葉斯方法估計(jì)λ的概率分布。對(duì)于第一次壩頂開裂，λ的先驗(yàn)分布為均勻分布，其中λmin和λmax分別取0 d 和300 d。圖7 為歷次壩頂開裂后λ的后驗(yàn)分布概率密度曲線，表1給出了λ的統(tǒng)計(jì)特征值。

由圖7 可見，在有限的觀測(cè)信息下，λ在第一、第二、第三和第四次壩頂破裂后的后驗(yàn)分布具有明顯的重尾特征。隨著信息（即壩頂開裂前的異常狀態(tài)持續(xù)時(shí)間）的積累，λ的概率分布的離散程度逐漸減小。從表1也可以看出，隨著觀察到的TAS的增加，λ的標(biāo)準(zhǔn)差從72 d 下降到41.4 d。因此，觀察信息的增加可以顯著降低λ的不確定性，這與統(tǒng)計(jì)學(xué)原理是一致的。

表1 參數(shù)λ的統(tǒng)計(jì)特征值Tab.1 Statistical characteristics of λ

圖7 參數(shù)λ概率分布的動(dòng)態(tài)更新結(jié)果Fig.7 Progressively updated distributions of λ

3.3 壩頂開裂風(fēng)險(xiǎn)的動(dòng)態(tài)評(píng)估和預(yù)測(cè)預(yù)警

為了對(duì)壩頂開裂風(fēng)險(xiǎn)進(jìn)行動(dòng)態(tài)評(píng)估，基于更新后的λ的概率分布，根據(jù)式（3），采用蒙特卡羅法對(duì)壩頂時(shí)變開裂概率進(jìn)行計(jì)算。圖8給出了歷次壩頂開裂的時(shí)變風(fēng)險(xiǎn)概率的均值、5%和95%分位值，以及觀察到的實(shí)際開裂時(shí)間和相應(yīng)的開裂概率均值，其中時(shí)間軸的原點(diǎn)對(duì)應(yīng)壩頂開裂前變形異常狀態(tài)的開始時(shí)刻。由圖8 可見，開裂概率隨異常變形狀態(tài)持續(xù)時(shí)間的延長(zhǎng)而增大。隨著觀測(cè)數(shù)據(jù)的積累，由于λ的不確定性逐漸減小，開裂概率估計(jì)值的置信區(qū)間逐漸縮小，開裂風(fēng)險(xiǎn)評(píng)估的不確定性降低。值得注意的是，與實(shí)際開裂時(shí)間對(duì)應(yīng)的開裂概率的均值不是一個(gè)確定值。

圖8 壩頂歷次開裂前估計(jì)的開裂時(shí)變概率Fig.8 Estimated time-varying cracking probability

在給定開裂概率的情況下，可以利用公式（4）估計(jì)壩頂開裂時(shí)間。因此，開裂時(shí)間的預(yù)測(cè)精度取決于預(yù)警所采用的開裂概率閾值。由于與實(shí)際開裂時(shí)間對(duì)應(yīng)的開裂概率具有隨機(jī)性，因此，根據(jù)λ的概率分布，計(jì)算開裂概率閾值的均值和置信區(qū)間。計(jì)算結(jié)果表明，壩頂開裂時(shí)刻對(duì)應(yīng)的開裂概率的均值約為50%，5%和95%分位值分別約為25%和80%。

根據(jù)開裂概率閾值的均值和5%、95%分位值，結(jié)合逐步更新的λ均值（見表1），采用公式（4）預(yù)測(cè)壩頂開裂時(shí)間。開裂時(shí)間預(yù)測(cè)值的均值、上限和下限如圖9所示，并與實(shí)際開裂時(shí)間進(jìn)行了比較。結(jié)果表明，壩頂實(shí)際開裂時(shí)間介于預(yù)測(cè)開裂時(shí)間的上限和下限之間，與開裂時(shí)間預(yù)測(cè)值的均值較為接近。利用預(yù)測(cè)的開裂時(shí)間的下限進(jìn)行壩頂開裂預(yù)警，可以看出，在壩頂發(fā)生開裂前20～70 d（即實(shí)際開裂時(shí)間與預(yù)測(cè)的開裂時(shí)間之間的間隔）可發(fā)出預(yù)警。因此，大壩安全管理單位可以根據(jù)壩頂開裂的早期預(yù)警，提前做好裂縫修復(fù)的準(zhǔn)備工作。

圖9 壩頂開裂時(shí)間的預(yù)測(cè)Fig.9 Prediction of dam crest cracking time

4 結(jié)論

本文針對(duì)傳統(tǒng)大壩安全風(fēng)險(xiǎn)評(píng)估方法存在定量化不足，不能實(shí)時(shí)跟蹤大壩安全風(fēng)險(xiǎn)的動(dòng)態(tài)變化，可能會(huì)錯(cuò)過故障發(fā)現(xiàn)及解決的最佳時(shí)間等問題，提出了集時(shí)變風(fēng)險(xiǎn)分析和參數(shù)動(dòng)態(tài)辨識(shí)于一體的高土石壩壩頂開裂風(fēng)險(xiǎn)動(dòng)態(tài)評(píng)估與預(yù)警方法，實(shí)現(xiàn)了壩頂開裂時(shí)變概率的動(dòng)態(tài)分析和基于概率的壩頂開裂風(fēng)險(xiǎn)提前預(yù)警，為高土石壩壩頂開裂風(fēng)險(xiǎn)預(yù)測(cè)預(yù)警和主動(dòng)調(diào)控提供了科學(xué)依據(jù)。

通過對(duì)某高心墻堆石壩壩頂開裂的風(fēng)險(xiǎn)評(píng)估和預(yù)警，對(duì)該方法進(jìn)行了驗(yàn)證。結(jié)果表明，隨著觀察到的壩頂開裂信息的積累，風(fēng)險(xiǎn)模型參數(shù)λ后驗(yàn)分布的離散程度逐漸減小。壩頂開裂概率隨異常變形狀態(tài)持續(xù)時(shí)間的延長(zhǎng)而增大，開裂概率估計(jì)值的置信區(qū)間隨λ分布的更新而逐漸縮小。通過對(duì)壩頂開裂時(shí)間下限的預(yù)測(cè)，可以在壩頂開裂發(fā)生前至少20 天發(fā)出預(yù)警，有利于大壩管理單位提前做好裂縫修復(fù)的準(zhǔn)備工作，提高大壩安全保障水平。