重力壩抗滑穩定可靠度動態評估模型研究與應用

甘珩佚，陳建康，高 山，吳震宇，周靖人

（1.四川大學水力學與山區河流開發保護國家重點實驗室水利水電學院，成都 610065；2.四川省水利水電勘測設計研究院有限公司，成都 610072）

0 引 言

混凝土重力壩因其施工工藝簡單，耐久性和抗滲性好等優點，在水利工程中得到廣泛應用［1］，如位于廣西的龍灘重力壩和位于四川的官地重力壩等［2］。據統計，重力壩失事事故中，因地基破壞導致失穩的比例高達40%［3］。受復雜運行環境和材料劣化等因素的影響，重力壩抗滑穩定安全具有時變特征，結合實時反映大壩安全性態的安全監測數據，實現監測-模擬數據的同化與融合，建立重力壩抗滑穩定安全風險-監測數據耦聯模型，實時跟蹤大壩抗滑穩定性態的動態演變，實現大壩穩定安全風險實時評估與預警，具有重要的理論意義與應用價值［4］。

重力壩抗滑穩定及其實時評估預測問題是近來的關注焦點［5］。楊志勇等改進了重力壩設計規范中的壩基抗滑穩定計算，以雙滑面深層抗滑穩定導出簡化的分項系數多滑面公式，取得了較好的預測效果［6］；彭圣軍利用優化算法融合混沌理論對大壩位移進行擬合和預測，構建了考慮殘差餛飩因子的混凝土壩位移混合監控模型，實時監控和預測大壩的安全性［7］；谷艷昌等以大壩潰決模式與壩體結構特性轉異為紐帶，利用有限元計算獲得相應的變形臨界值，進而建立了基于風險管理的大壩變形預警指標［8］。沈振中［9］等將基于塊體理論的非連續變形分析引入重力壩安全因素和預警指標研究中，根據失穩判據和壩體位移與強度折減系數之間的關系曲線確定了大壩的結構性態，并估計出大壩彈性狀態和承載極限狀態的變形預警指標。但結合實時監測數據的抗滑穩定安全評估方法鮮見報道，因此，如何利用監測數據實時對大壩安全進行評估，盡早對潛在的危險進行預警，將災害風險和損失降為最低是我們急需解決的問題。

本文結合壩基巖體蠕變與瞬時變形參數反分析技術，構建了重力壩抗滑穩定響應面功能函數，提出了基于大壩變形監測數據的重力壩抗滑穩定可靠度動態評估模型，并將其應用到GD重力壩抗滑穩定安全分析與預測中，驗證其有效性和適用性。

1 重力壩材料參數反分析方法

1.1 基于統計回歸模型的重力壩變形時效位移分量提取

在重力壩的變形統計回歸預測模型中，其因子函數主要由時效分量、溫度分量與水壓分量構成。因此，其統計回歸模型表達式如式（1）所示。

式中：Y(t)為壩體變形監測值在時間t的統計估計值；F1[θ(t)]為壩體變形的時效分量；F2[θ(t)]為壩體變形的溫度分量；F3[H(t)]為壩體變形的水壓分量。

時效分量是隨時間推移而朝某一方向發展的不可逆的分量，一般采用指數式、對數式、雙曲線式等表示，本文采用的時效因子函數為對數線性式［10］，即：

式中：ti為觀測的總天數；c1、c2為待回歸系數。

溫度對重力壩位移的影響主要表現在氣溫的季節性變化上，本文采用的溫度因子函數為多段平均氣溫的線性組合［11］。

式中：Ti為觀測日當天、前3 d、前5 d、前10 d、前15 d、前30 d 的平均大氣氣溫。

水壓分量一般與重力壩變形的轉角、扭矩等有關，本文采用的水壓分量函數為：

式中：H(t)為變形觀測日的平均水位；H(t0)為基準水位；Hi(t)為觀測日當天、前1～3 d、前4～7 d、前8～15 d、前16～30 d、前31～60 d平均水位；ai為待回歸系數［12］。

綜上所述，經過上述對各因子函數的選取，建立重力壩變形監測分析數學模型為：

基于上述模型，重力壩變形的時效位移分量可表達為：

1.2 基于時效位移的壩與地基材料參數反分析

按最小二乘法可構建參數反分析的目標函數，如式（7）所示：

式中：E為反分析目標函數值；x為待求參數，其上下限由統計資料得到；δ*i為測點i的實測位移；δi(x)為測點i的響應面函數輸出位移；k為測點總數［13］。

在混凝土壩-地基系統瞬時-蠕變參數反分析過程中，在不同時刻每個監測點都有不同的監測值，因此需要構建響應面方程序列來擬合壩體在時間過程中的變形特征，目標函數可以表示為：

式中：F為目標函數；δ*i(tj)為監測點i在不同時刻下的實測位移；δi(x1，x2，…，xn，tj)為監測點i的響應面函數輸出位移歷時過程；x1，x2，…，xn為E-B 模型和蠕變模型中待反演的參數；k為選取的目標監測點數量；m為選取的時間節點總數［14］。

結合式（8）所示目標函數，將遺傳算法運用在混凝土壩-地基系統瞬時-蠕變參數反分析中，通過參數編碼、基因重組和適應性選擇推動族群進化從而得到最優參數。本文模擬蠕變采用的本構模型為Burgers蠕變本構，模擬瞬時變形采用DP本構模型。

影響混凝土重力壩變形的材料參數眾多，傳統的位移反分析采用對所有參數進行反分析，導致計算效率低，且參數之間的相互影響可能導致次要影響因素遮蔽主要影響因素。為此，本文針對影響大壩變形的巖體材料物理力學參數采用極差分析法和方差分析法確定影響大壩變形的關鍵參數。本文基于巖石蠕變特性與監測時效變形的物理意義，利用大壩位移統計回歸模型分離出的時效分量，建立基于時效變形的重力壩參數反分析方程，其具體步驟如下：

（1）根據敏感性分析的結果擇取對大壩變形有高度影響的參數作為反演目標參數；

（2）建立大壩位移監測點位移的統計方程，分離出監測點位移的時效分量，作為反分析的目標值；

（3）根據設計資料中的材料設計參數來確定反演參數的合理取值范圍。

（4）采用正交試驗法進行試驗設計，通過參數分布類型抽樣或根據建議取值范圍確定試驗中參數水平。

（5）進行每組試驗的分析計算，統計整理每個監測點的數值模擬結果。

（6）構建每個序列典型時刻位移差響應面方程，求解響應面方程的系數。

（7）以監測點位移的時效為目標值，構建反分析的目標函數，采用遺傳算法對目標函數進行求解，得到反分析參數的最優解。

（8）利用反分析得到的影響大壩變形關鍵參數的優化解，進行大壩變形分析，比較模擬值與實測值的誤差，若誤差精度滿足要求，則可終止反分析循環，得到的參數即為目標優化解；若精度不滿足要求，重復步驟（3）～（7），至精度滿足要求為止。

2 基于監測數據的重力壩抗滑穩定可靠度預測模型構建

2.1 抗滑穩定功能函數構建

對重力壩與地基結構體系而言，設其隨機變量為x1、x2、…、xn，定義功能函數G=g(x1，x2，…，xn)來表征結構與地基的工作狀態，則其極限狀態可以表述為

在重力壩抗滑穩定分析中，建基面或者深層滑動面動力大于阻滑力均屬超承載能力極限狀態。在有限元分析中重力壩抗滑穩定可用抗剪斷公式表達，因此重力壩抗滑失穩功能函數可構建為：

式中：f′i、ci為滑裂面第i個單元的抗剪強度指標；σi為第i個單元滑裂面方向受到的正應力；τi為第i個單元滑裂面方向受到的切應力；li為第i個單元沿滑裂面方向的單元長度；n為整個滑裂面上所有單元的總數［15］。

基于式（10），采用蒙特卡洛法計算抗滑穩定可靠度，即隨機抽取隨機變量X的樣本值{x1，x2，…，xn}T，以功能函數值G(X1，X2，…，Xn)≤0 的次數占總抽樣次數比獲得失效風險Pf。為便于分析，此處可以構造一個取值為0和1的函數，即：

式中：Df表征結構的失效域，即G(x)∈Df時，結構失效；反之，結構可靠。

因此，壩與地基失效風險可以表達為：

式中：為壩與地基失效風險的抽樣估計值；N為總抽樣數；G(x)i為第i次抽樣對應的功能函數值。

綜上，重力壩與地基抗滑穩定安全的可靠指標為

2.2 基于監測數據的抗滑穩定可靠度預測模型構建

大壩安全狀態最直觀的反映就是監測數據，測點測值的大小、變化趨勢等與大壩穩定安全狀態密切相關。以重力壩的安全監測數據為基礎，結合數值模擬成果，以測點變形監測序列為自變量，重力壩抗滑穩定可靠度為因變量，構建重力壩的穩定可靠度評估模型：

式中：F[Y(t)]為測點變形分量；βs(t)為重力抗滑穩定可靠指標在時間t的統計估計值；C為待定常數項。

式中：ai為待定系數；yi為變形序列的監測值；n為變形序列的個數。

實際工程中，大壩變形測點測值一般包括四部分，即水壓分量、時效分量、環境量影響及系統誤差等，因此：

式中：F1[H(t)]為監測點水壓分量；F2[T(t)]為溫度分量；F3[θ(t)]為時效分量；C為待定常數項。

因本次有限元計算中未考慮溫度因素，僅利用時效和水壓分量構建可靠度預測模型，故將上述方程中變形值的時效和水壓分量分離出來構建大壩抗滑穩定可靠度預測模型。

式中各分量定義同上。

本文提出的基于監測數據的抗滑穩定可靠度預測模型構建流程如圖1所示。

3 工程應用

3.1 工程概況

GD 重力壩位于四川省涼山彝族自治州西昌市和鹽源縣交界的打羅村境內。主體樞紐主要由攔河碾壓混凝土重力壩、泄洪消能建筑物、引水發電建筑物等組成，工程為一等大（1）型工程，壩軸線方向N12°E，攔河壩壩頂高程1 334 m，壩頂長度516 m。水庫正常蓄水位1 330.00 m，設計洪水位1 330.18 m，校核洪水位1 330.44 m，死水位1 328.00 m，極限死水位1 321.00 m，總庫容7.6 億m3，屬日調節水庫。電站裝機容量2 400 MW，多年平均發電量118.7 億kWh。

本次反演分析選用GD 重力壩的9#壩段，該壩段為一擋水壩段，建基面高程1 186 m，壩頂高程1 334 m，最大壩高148 m，工程地質條件較復雜。倒垂測點IP3 和雙金屬標測點DS3 埋設在1 205 m 的大壩下部，最低端位于1 155 m，與基巖緊密結合；真空激光測點LA9、LA1分別布置于高程1 334 m 的壩頂和高程1 205 m的大壩底部廊道內。GD水電站相關監測儀器布置圖如圖2所示。

3.2 有限元模型

有限元模型的模擬范圍為：自建基面向下延伸2.5 倍壩高（約370 m），上下游方向同樣延伸2.5 倍壩高（約370 m）。為對比計算位移與實測位移，根據真空激光測點LA9、倒垂線測點IP3、雙金屬標DS3的埋設位置，在有限元模型中布料，其中壩體由C15、C20、C25 三種標號的碾壓混凝土和C20 標號的常態混凝土組成，壩基由Ⅱ類、Ⅲ類、Ⅳ類、Ⅴ類巖體、斷層和破碎帶組成。有限元模型共劃分為65 190 個單元和73 964 個結點。大壩三維模型順河向X軸以指向下游為正，鉛直向Y軸以垂直向上為正，橫河向Z軸以指向右岸為正，模型及材料分區如圖3所示。

3.3 重力壩壩基巖體參數反分析

3.3.1 參數敏感性分析

本文模擬蠕變采用的本構模型為Burgers蠕變本構，該模型主要有Maxwell 剪切模量mshear、Kelvin 剪切模量kshear、Maxwell 黏滯系數mviscosity、Kelvin 黏滯系數kviscosity 四個參數，模擬瞬時變形采用DP 本構模型，該模型主要有E彈模，c凝聚力，φ摩擦角3個參數。而GD重力壩9號壩段共4類巖體8個蠕變參數，15個力學參數。根據GD重力壩監測點的布置情況，選取9號壩段的壩頂激光位移測點LA9以及底部廊道中的激光測點LA1 作為此次分析的特征點（如圖4）。通過數值模擬，采用極差分析法和方差分析法分別對特征點的順河向位移、垂直位移和綜合位移三個目標量進行敏感性分析。

極差分析法和方差分析法顯示，壩體C25 碾壓混凝土的E；C20 碾壓混凝土的E、c；II 類巖體的E、c、φ、mshear、kshear、mviscosity、kviscosity；Ⅲ類巖體的E、c對大壩變形均有影響或有顯著影響。因此，選擇以上參數作為此次反分析參數。

3.3.2 壩體與地基參數反分析

選取9 號擋水壩段的壩頂激光位移測點LA9、底部廊道中的倒垂測點IP3 和雙金屬標測點DS3 作為參數反演分析的特征點，選取目標序列2014.8.1、2015.8.1、2016.8.1 和2017.8.1 的測值為目標值進行參數反分析，蠕變參數與DP 模型參數反演結果見表1。

表1 反演參數表Tab.1 Inversion parameter table

對所有序列進行誤差對比分析，其反分析參數模擬值與實測值最大誤差均在19%以內；選擇時間序列的中間節點2016.12.31 和末尾端點2018.12.31 作為誤差分析的結果展示，根據表1中大壩變形參數的優化解和初始設計值，分別對GD重力壩的變形進行數值模擬分析，如表2所示，通過有限元模擬與測值比較分析可見，采用設計參數的變形模擬值與實測值的最大誤差50.21%，平均誤差為25.49%。采用反分析參數模擬值與實測值的最大誤差17.43%，平均誤差10.56%，誤差在合理范圍內，曲線吻合較好，表明參數反分析效果良好，可以反映大壩的變形性態。

表2 參數反演誤差統計表Tab.2 Statistical table of parameter inversion error

3.4 抗滑穩定可靠度預測模型構建

3.4.1 失效通道搜尋

本文采用荷載增量法搜索最危險的滑移通道。GD 重力壩壩基滑移通道如圖5所示。

3.4.2 抗滑穩定可靠度分析

本文確定材料的特征值如表5，壩基巖體的彈模E、凝聚力c和內摩擦角φ以及壩體C25 混凝土、C20 混凝土的彈模E考慮為隨機變量。

表3 材料參數隨機特征Tab.3 Random characteristics of material parameters

采用不含交叉項的二次多項式，構建GD 重力壩抗滑穩定安全系數響應面方程。采用正交試驗法，獲得不同時刻大壩抗滑穩定響應面方程系數。本次分析的時間序列為2012.4.1-2018.12.31，在該時間段內每隔3 個月建立一個響應面方程，共計建立28 個響應面方程來評估其運行期可靠度，本次28 個響應面方程擬合效果較好，復相關系數均大于0.88，響應面方程具有良好的顯著性。

進行大壩抗滑穩定可靠度分析時需要對參數進行截尾處理，參數在均值的3σ范圍內取值有99.7%的保證率，故選擇μ-3σ、μ+ 3σ分別為隨機變量的左右截尾點，根據關鍵滑移通道響應面方程，采用蒙特卡洛法獲得大壩抗滑穩定的可靠指標β與失效概率Pf歷時曲線見圖6。由圖6可見，大壩初期蓄水階段可靠指標相對較低，2014年后水位僅在小范圍內波動，可靠指標也在小范圍內呈現相似的波動趨勢?，F階段可靠指標達到4.27，其變化趨于穩定。

3.4.3 抗滑穩定可靠度預測模型構建

基于本文提出的重力壩抗滑穩定可靠度動態評估模型和式（13），集合2012.4.1-2018.12.31 間的可靠度模擬值，構建GD重力壩可靠度預測模型，如圖7所示。結果顯示，復相關系數大于0.93，可靠指標計算與擬合值平均誤差0.79%，最大誤差1.55%，精度較高，預測模型如式（14）所示。通過對后180 d 和360 d 大壩抗滑穩定可靠指標預測分析，其平均誤差和最大誤差如表6所示都在5%以內。工程應用表明，本文提出的基于監測數據的重力壩抗滑穩定可靠度動態評估模型效率和精度較高，工程應用價值明顯。

式中：y1為IP3 順河向位移時效分量；y2為DS3 垂直位移時效分量；y3為LA9 順河向位移時效分量；y4為LA9 垂直位移時效分量；y5為LA9 順河向位移水壓分量；y6為LA9 垂直位移水壓分量；y7為DS3 垂直位移水壓分量；y8為IP3 順河向位移水壓分量。

表4 大壩抗滑穩定可靠指標預測精度統計Tab.4 Prediction accuracy statistics of reliability index of dam anti-sliding stability

4 結 論

本文通過對重力壩與地基材料參數反分析，抗滑穩定可靠度及其動態預測模型的深入研究，獲得以下幾點認識：

（1）針對壩與地基材料參數眾多，參數設計值難以準確表征材料性能等問題，采用極差分析和方差分析法，提取影響重力壩變形的關鍵影響參數，通過時效位移分離、目標函數構建和參數優化求解，提出了基于時效位移的混凝土與巖體參數反分析方法。

（2）采用荷載增量法確定重力壩抗滑穩定危險滑面，基于蒙特卡洛抽樣模擬，通過抗滑穩定安全時變可靠度分析，構建了抗滑穩定可靠度-監測數據耦聯預測模型，提出了重力壩抗滑穩定實時評估方法。

（3）GD工程應用表明，基于反演參數得到的位移模擬結果與監測數據平均誤差小于11%，表明本文提出的方法精度較高。提出的抗滑穩定可靠度預測模型平均誤差為0.79%，最大誤差為1.55%，精度高，適用性好。□