重力壩抗滑穩(wěn)定可靠度動態(tài)評估模型研究與應(yīng)用

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

（1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室水利水電學(xué)院，成都 610065；2.四川省水利水電勘測設(shè)計研究院有限公司，成都 610072）

0 引 言

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

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

本文結(jié)合壩基巖體蠕變與瞬時變形參數(shù)反分析技術(shù)，構(gòu)建了重力壩抗滑穩(wěn)定響應(yīng)面功能函數(shù)，提出了基于大壩變形監(jiān)測數(shù)據(jù)的重力壩抗滑穩(wěn)定可靠度動態(tài)評估模型，并將其應(yīng)用到GD重力壩抗滑穩(wěn)定安全分析與預(yù)測中，驗證其有效性和適用性。

1 重力壩材料參數(shù)反分析方法

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

在重力壩的變形統(tǒng)計回歸預(yù)測模型中，其因子函數(shù)主要由時效分量、溫度分量與水壓分量構(gòu)成。因此，其統(tǒng)計回歸模型表達(dá)式如式（1）所示。

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

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

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

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

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

水壓分量一般與重力壩變形的轉(zhuǎn)角、扭矩等有關(guān)，本文采用的水壓分量函數(shù)為：

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

綜上所述，經(jīng)過上述對各因子函數(shù)的選取，建立重力壩變形監(jiān)測分析數(shù)學(xué)模型為：

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

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

按最小二乘法可構(gòu)建參數(shù)反分析的目標(biāo)函數(shù)，如式（7）所示：

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

在混凝土壩-地基系統(tǒng)瞬時-蠕變參數(shù)反分析過程中，在不同時刻每個監(jiān)測點都有不同的監(jiān)測值，因此需要構(gòu)建響應(yīng)面方程序列來擬合壩體在時間過程中的變形特征，目標(biāo)函數(shù)可以表示為：

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

結(jié)合式（8）所示目標(biāo)函數(shù)，將遺傳算法運(yùn)用在混凝土壩-地基系統(tǒng)瞬時-蠕變參數(shù)反分析中，通過參數(shù)編碼、基因重組和適應(yīng)性選擇推動族群進(jìn)化從而得到最優(yōu)參數(shù)。本文模擬蠕變采用的本構(gòu)模型為Burgers蠕變本構(gòu)，模擬瞬時變形采用DP本構(gòu)模型。

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

（1）根據(jù)敏感性分析的結(jié)果擇取對大壩變形有高度影響的參數(shù)作為反演目標(biāo)參數(shù)；

（2）建立大壩位移監(jiān)測點位移的統(tǒng)計方程，分離出監(jiān)測點位移的時效分量，作為反分析的目標(biāo)值；

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

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

（5）進(jìn)行每組試驗的分析計算，統(tǒng)計整理每個監(jiān)測點的數(shù)值模擬結(jié)果。

（6）構(gòu)建每個序列典型時刻位移差響應(yīng)面方程，求解響應(yīng)面方程的系數(shù)。

（7）以監(jiān)測點位移的時效為目標(biāo)值，構(gòu)建反分析的目標(biāo)函數(shù)，采用遺傳算法對目標(biāo)函數(shù)進(jìn)行求解，得到反分析參數(shù)的最優(yōu)解。

（8）利用反分析得到的影響大壩變形關(guān)鍵參數(shù)的優(yōu)化解，進(jìn)行大壩變形分析，比較模擬值與實測值的誤差，若誤差精度滿足要求，則可終止反分析循環(huán)，得到的參數(shù)即為目標(biāo)優(yōu)化解；若精度不滿足要求，重復(fù)步驟（3）～（7），至精度滿足要求為止。

2 基于監(jiān)測數(shù)據(jù)的重力壩抗滑穩(wěn)定可靠度預(yù)測模型構(gòu)建

2.1 抗滑穩(wěn)定功能函數(shù)構(gòu)建

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

在重力壩抗滑穩(wěn)定分析中，建基面或者深層滑動面動力大于阻滑力均屬超承載能力極限狀態(tài)。在有限元分析中重力壩抗滑穩(wěn)定可用抗剪斷公式表達(dá)，因此重力壩抗滑失穩(wěn)功能函數(shù)可構(gòu)建為：

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

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

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

因此，壩與地基失效風(fēng)險可以表達(dá)為：

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

綜上，重力壩與地基抗滑穩(wěn)定安全的可靠指標(biāo)為

2.2 基于監(jiān)測數(shù)據(jù)的抗滑穩(wěn)定可靠度預(yù)測模型構(gòu)建

大壩安全狀態(tài)最直觀的反映就是監(jiān)測數(shù)據(jù)，測點測值的大小、變化趨勢等與大壩穩(wěn)定安全狀態(tài)密切相關(guān)。以重力壩的安全監(jiān)測數(shù)據(jù)為基礎(chǔ)，結(jié)合數(shù)值模擬成果，以測點變形監(jiān)測序列為自變量，重力壩抗滑穩(wěn)定可靠度為因變量，構(gòu)建重力壩的穩(wěn)定可靠度評估模型：

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

式中：ai為待定系數(shù)；yi為變形序列的監(jiān)測值；n為變形序列的個數(shù)。

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

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

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

式中各分量定義同上。

本文提出的基于監(jiān)測數(shù)據(jù)的抗滑穩(wěn)定可靠度預(yù)測模型構(gòu)建流程如圖1所示。

3 工程應(yīng)用

3.1 工程概況

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

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

3.2 有限元模型

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

3.3 重力壩壩基巖體參數(shù)反分析

3.3.1 參數(shù)敏感性分析

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

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

3.3.2 壩體與地基參數(shù)反分析

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

表1 反演參數(shù)表Tab.1 Inversion parameter table

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

表2 參數(shù)反演誤差統(tǒng)計表Tab.2 Statistical table of parameter inversion error

3.4 抗滑穩(wěn)定可靠度預(yù)測模型構(gòu)建

3.4.1 失效通道搜尋

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

3.4.2 抗滑穩(wěn)定可靠度分析

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

表3 材料參數(shù)隨機(jī)特征Tab.3 Random characteristics of material parameters

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

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

3.4.3 抗滑穩(wěn)定可靠度預(yù)測模型構(gòu)建

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

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

表4 大壩抗滑穩(wěn)定可靠指標(biāo)預(yù)測精度統(tǒng)計Tab.4 Prediction accuracy statistics of reliability index of dam anti-sliding stability

4 結(jié) 論

本文通過對重力壩與地基材料參數(shù)反分析，抗滑穩(wěn)定可靠度及其動態(tài)預(yù)測模型的深入研究，獲得以下幾點認(rèn)識：

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

（2）采用荷載增量法確定重力壩抗滑穩(wěn)定危險滑面，基于蒙特卡洛抽樣模擬，通過抗滑穩(wěn)定安全時變可靠度分析，構(gòu)建了抗滑穩(wěn)定可靠度-監(jiān)測數(shù)據(jù)耦聯(lián)預(yù)測模型，提出了重力壩抗滑穩(wěn)定實時評估方法。

（3）GD工程應(yīng)用表明，基于反演參數(shù)得到的位移模擬結(jié)果與監(jiān)測數(shù)據(jù)平均誤差小于11%，表明本文提出的方法精度較高。提出的抗滑穩(wěn)定可靠度預(yù)測模型平均誤差為0.79%，最大誤差為1.55%，精度高，適用性好?！?/p>