基于貝葉斯參數更新的高土石壩壩頂開裂風險動態評估與預警

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

（1.四川大學水利水電學院，四川成都 610065；2.四川大學水力學及山區河流開發保護國家重點實驗室，四川成都 610065；3.國能大渡河流域水電開發有限公司，四川成都 610041）

0 引言

土石壩對地形地質條件適應性強，可以充分利用當地材料，造價較低，是水電開發中最具應用前景的壩型之一［1，2］。在我國已建土石壩中，潰壩事故主要集中在中低高度大壩，高土石壩主要表現為壩頂開裂、壩體與壩肩（基）滲漏、壩坡滑坡塌坑等非潰壩故障模式，根據相關文獻統計，壩頂開裂約占高心墻堆石壩常見故障模式的40%以上［3-5］。

眾所周知，壩頂開裂與土石壩的變形特性密切相關，特別是變形傾度指標［6-8］。通過大壩原觀監測數據，結合參數反演算法和數值模擬方法可以對高心墻堆石壩運行形態進行準確預測［9，10］。當基于數值模擬結果計算出的壩頂變形傾度大于某一閾值時，壩頂存在開裂可能性，但是并不意味著壩頂會立即開裂［11，12］，若壩頂變形傾度超過閾值的異常狀態持續時間較長，壩頂開裂的可能性就會隨著異常狀態的持續時間而增大。由于復雜的物理機制和隨機因素，壩頂開裂前的異常狀態持續時間是不確定的，可以看作是一個隨機變量，并可以基于此對異常狀態期間發生壩頂開裂的風險進行評估。

由于缺乏相關的歷史資料，特別是對于某一具體的土石壩，無法通過經典的頻率統計分析得到壩頂開裂前異常狀態持續時間的概率分布參數（PDP）。根據異常狀態開始到發生峰值之間的觀測時間間隔，可以利用貝葉斯理論估計PDP 的概率分布［13］。此外，PDP 的概率分布可以隨著觀察到的證據的積累過程而動態更新。因此，在貝葉斯框架下，土石壩的壩頂開裂風險可以根據PDP 的逐步更新概率分布進行動態評估。此外，還可以根據概率估計的開裂時間范圍對壩頂開裂進行預警。

本文研究的主要問題是：①如何用壩頂沉降異常狀態持續時間來評價土石壩壩頂開裂風險；②如何估計給定概率區間對應的壩頂開裂次數以進行壩頂開裂預警；③如何更新壩頂開裂風險公式中參數λ的概率分布。

1 壩頂開裂風險動態評估與預警方法

1.1 總體思路與方法架構

首先確定高心墻堆石壩壩頂異常狀態判定準則。如引言所述，壩頂開裂與壩頂變形傾度γS密切相關。當γS超過某一閾值時，壩頂易發生開裂故障。變形傾度的定義是［14］：

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

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

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

壩頂開裂概率計算中的參數是不確定的。如何估計參數的概率分布是該方法應用中的一個主要問題。應用貝葉斯理論，根據觀測到的從出現異常狀態到出現壩頂開裂的時間間隔，估計和更新參數的概率分布。

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

圖1 壩頂開裂動態風險評價與預警流程圖Fig.1 Flowchart of the dynamic risk evaluation and early warning of crest cracking

1.2 壩頂開裂概率與預警

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

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

當壩頂沉降的異常狀態（即γS＞γC）發生并持續時，在給定的時間內出現壩頂開裂的概率為：

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

1.3 信息驅動的風險模型參數更新

參數λ中的認知不確定性源于分析人員對λ的統計信息的缺乏。這種認知不確定性可以用貝葉斯框架中的主觀概率定量地表示。當獲得與參數λ有關的信息時，可以根據貝葉斯定理估計和更新λ的概率分布如下：

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

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

（2）第二類信息：從異常狀態出現到異常狀態結束，未出現壩頂開裂的持續時間，相應的似然函數表述為：

根據新觀測信息的類型，將相應的似然函數和λ的先驗概率分布代入公式（5），即可更新λ的概率分布。當第一次估計λ的概率分布時，沒有任何關于λ的先驗信息，此時采用均勻分布來表示π0(λ)：

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

從λ的后驗概率分布中抽樣，將λ的樣本代入公式（3）和（4），可以計算出壩頂開裂概率和開裂時間的均值和置信區間。

2 工程案例

2.1 工程概況

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

圖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）和正常水位之間周期性波動。大壩填筑過程和庫水位變化過程如圖3所示。

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

2.2 大壩壩頂開裂情況

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

3 案例分析結果

3.1 壩頂時序變形傾度及異常狀態識別

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

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

圖5 實測水平及沉降位移與計算值的對比Fig.5 Comparison of the computed and measured horizontal and settlement displacements

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

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

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

3.2 開裂風險模型參數λ概率分布的貝葉斯估計

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

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

表1 參數λ的統計特征值Tab.1 Statistical characteristics of λ

圖7 參數λ概率分布的動態更新結果Fig.7 Progressively updated distributions of λ

3.3 壩頂開裂風險的動態評估和預測預警

為了對壩頂開裂風險進行動態評估，基于更新后的λ的概率分布，根據式（3），采用蒙特卡羅法對壩頂時變開裂概率進行計算。圖8給出了歷次壩頂開裂的時變風險概率的均值、5%和95%分位值，以及觀察到的實際開裂時間和相應的開裂概率均值，其中時間軸的原點對應壩頂開裂前變形異常狀態的開始時刻。由圖8 可見，開裂概率隨異常變形狀態持續時間的延長而增大。隨著觀測數據的積累，由于λ的不確定性逐漸減小，開裂概率估計值的置信區間逐漸縮小，開裂風險評估的不確定性降低。值得注意的是，與實際開裂時間對應的開裂概率的均值不是一個確定值。

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

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

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

圖9 壩頂開裂時間的預測Fig.9 Prediction of dam crest cracking time

4 結論

本文針對傳統大壩安全風險評估方法存在定量化不足，不能實時跟蹤大壩安全風險的動態變化，可能會錯過故障發現及解決的最佳時間等問題，提出了集時變風險分析和參數動態辨識于一體的高土石壩壩頂開裂風險動態評估與預警方法，實現了壩頂開裂時變概率的動態分析和基于概率的壩頂開裂風險提前預警，為高土石壩壩頂開裂風險預測預警和主動調控提供了科學依據。

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