心墻堆石壩非飽和滲流模型識別方法及應用

許云鵬，吳震宇，尹 川

(四川大學水利水電學院 水力學與山區河流開發保護國家重點實驗室，四川 成都 610065)

1 概述

土質心墻堆石壩是我國西部高壩廣泛采用的壩型之一，已建成的長河壩(240m)、糯扎渡壩(261.5m)、兩河口壩(293m)等心墻堆石壩工程運行良好，在建的古水壩(305m)、雙江口壩(314m)是超過300m的世界級水平大型心墻堆石壩工程[1]。防滲心墻的滲透穩定性是影響大壩安全的關鍵因素。庫水位的變動會影響心墻堆石壩的滲透穩定性。一般來說，堆石壩心墻由浸潤線分成飽和區和非飽和區，滲透系數在飽和區域是一個常數[2- 3]。庫水位變化時，總存在一個非飽和區域，其滲透系數隨孔隙水壓力和給水度的變化而變化。但在通常的心墻堆石壩滲流數值模擬計算中，視滲透系數為常數，上述因素對計算的影響往往被忽略[4- 5]。實踐證明，在水庫水位變化情況下計算土壩滲流時，如果不考慮心墻非飽和特性因素，計算結果和實際情況相比會有很大的誤差[2，6]。

非飽和滲透系數可根據土壤水分特征曲線(SWCC)進行推定[6]。土壤水分特征曲線成為非飽和滲流數值模擬的重要基礎。直接測定土壤水分特征曲線成本較高且十分復雜，經驗公式法預測SWCC逐漸應用開來。經驗公式法是利用代數關系式來描述SWCC，代表性的模型主要有Brooks-Corey模型(BC模型)、van Genuchten模型(VG模型)、Campbell模型、Fredlund模型等[7]。在眾多SWCC模型中，因為BC模型和VG模型適用性范圍更廣，準確度更高，故而相關研究中得到了廣泛的應用[8]。同時BC模型和VG模型提供滲透系數函數模型，用于表示滲透系數隨基質吸力的變化情況。Gardner模型[9]僅提供滲透系數函數，可與VG模型的SWCC結合，用于非飽和滲流分析，稱為VG-G模型。

土-水特征曲線和滲透系數函數可以準確反映材料飽和度隨孔隙水壓力的變化情況，進而反映材料滲透系數隨孔隙水壓力的變化情況，直接影響非飽和非穩定滲流模型的正確性。因此基于滲流實測數據如何合理確定SWCC和滲透系數函數模型，是構建非飽和滲流模型、開展后續滲流分析研究的關鍵。模型選擇不僅要考慮模型的復雜程度，也要考慮模型的預測效果，過于簡單的模型不能充分描述數據的內在結構變化，造成預測偏差，模型復雜雖然能充分描述數據的內在關聯，但過于復雜的模型也會導致預測的方差偏大，因此選擇合適的模型就是要找到預測偏差與方差的平衡點。1974年日本學者Akaike在信息理論研究中，為平衡估計模型的復雜性和擬合數據優良性，提出了AIC準則[10]。1978年Schwarz使用蒙特卡洛法研究AIC準則時發現，當樣本數量增加時，AIC準則計算精度逐漸下降，進而在AIC準則中引入了包含樣本數量的懲罰項，提出BIC準則以滿足大樣本條件下的模型識別[11]。隨著模型識別方法的逐漸完善，越來越多的國內外學者圍繞模型識別開展了深入研究。汪建均[12]在識別非正態響應部分因子試驗的顯著性因子中，結合經驗貝葉斯先驗和變量后驗，提出了基于廣義線性模型的模型選擇方法。段偲默[13]根據AIC準則和BIC準則從基于Copula函數的風光聯合出力的靜態模型和基于Copula函數的風光聯合出力的動態模型中獲取最優模型。本文以PBG礫石土心墻堆石壩為例，基于原觀數據進行非飽和滲流參數反演然后基于AIC、BIC準則識別非飽和滲流模型。

2 非飽和滲流模型

達西定律被證明既適合飽和土體，也適合非飽和土體。與飽和土體相比，非飽和土體滲透系數不是一個常數，它隨著土體飽和度或基質吸力的變化而變化。根據達西定律及流體運動連續性方程可導出二維非穩定飽和—非飽和滲流控制方程為：

(1)

(2)

式中，kx、ky—x、y方向的滲透系數；H—總水頭；Q—邊界流量；γw—水的重度；mw—比水重度，定義為土水特征曲線斜率；θw—土體含水率；ua—孔隙氣壓力；uw—孔隙水壓力。

對于穩定流條件，基質吸力剖面的時間依賴性消失，等式(1)的右側變為零。以下介紹了本工作中用于描述飽和和非飽和流動的主要本構模型和參數。Brooks和Corey基于土體含水量試驗提出了有效飽和度和基質吸力之間關系的土水特征曲線模型：

(3)

式中，h—為吸力水頭；hb—進氣壓力水頭。

Van Genuchten提出了一種形式平滑、封閉的3參數土水特征曲線模型如下：

Se=[1+(-ah)n]-m

(4)

式中，a—約等于進氣壓力值的倒數，kPa-1；m、n—無量綱擬合參數，m=1-1/n。

滲透系數隨基質吸力變化而變化。BC模型和VG模型分別給出了滲透系數函數關系方程如下：

(5)

(6)

式中，Ks—飽和土滲透系數。

除上述2種模型外，Gardner提出了冪函數形式的雙參數滲透系數模型：

K(h)=Ks[1+(-αh)β]-1

(7)

式中，α—滲透系數開始下降的起點；β—滲透系數隨基質吸力下降的斜率。

將Brooks-Corey表達土水特征曲線的公式(3)和滲透系數的公式(5)稱為Brooks-Corey-Burdine模型，簡稱為Brooks-Corey模型(BC模型)；將van Genuchten表達土水特征曲線的公式(4)和滲透系數的公式(6)稱為van Genuchten-Mualem模型，簡稱為van Genuchten模型(VG模型)。Gardner模型僅提供滲透系數函數，可與VG模型的土水特征曲線結合，用于非飽和滲流分析，稱為VG-G模型。

3 模型識別方法

3.1 模型識別準則

最優模型準則包括似然函數最大化和模型未知參數最少化。似然函數值越大模型擬合效果越好，但易導致未知參數越來越多，模型過于復雜，因此最優模型應該綜合擬合精度和參數個數最優化。AIC準則全稱是赤池信息準則(Akaike Information Criterion)，主要思想是平衡模型的復雜度和數據的擬合精度，其表達式如下：

AIC=2k-2lnL

(8)

式中，k—參數個數，L—極大似然函數。

當樣本容量較大時，AIC準則會選擇參數較多的模型，造成過擬合現象。針對該問題，Schwartz提出BIC準則，全稱是貝葉斯信息準則(Bayesian Information Criterion)，其表達式如下：

BIC=klnN-2lnL

(9)

式中，N—樣本數量，其他參數意義同上式。

BIC準則考慮了樣本數量，避免精度過高選擇復雜模型。使用AIC、BIC準則進行模型識別時，值越小模型越優。

3.2 模型識別方法流程

首先選擇待反演的堆石壩心墻非飽和區滲透系數、飽和區滲透系數、非飽和滲流模型特性參數作為待反演參數。根據滲透系數和非飽和模型特性參數取值區間，在材料正常的取值范圍內上下波動參數進行正交試驗，得到多組試驗參數組合。根據壩體結構輪廓尺寸、材料分區和壩址地形地質條件等資料，構建二維心墻堆石壩有限元模型，輸入實測庫水位過程線，計算各組試驗組合下非穩定非飽和滲流數值模擬，輸出每組試驗下的測點模擬滲壓監測序列。以最低庫水位為界，將心墻分為飽和區和非飽和區。根據上游水位變化，選取代表性水位對應時間節點作為滲壓控制節點。以待反演參數為自變量，各組有限元模型控制節點下計算滲壓值為因變量擬合響應面方程。基于最小二乘法原理，利用心墻實際測點滲壓監測序列和擬合的響應面方程構建目標函數。將構造的多目標函數輸入到NSGA-II算法中進行優化，輸出Pareto最優解集，最終確定反演參數。反演參數代入有限元模型，計算得到心墻擬合滲壓值，結合實測滲壓值，采用AIC、BIC準則對每個模型進行評價，最終篩選出最優模型。模型識別方法流程如圖1所示。

圖1 模型識別方法

4 工程實例

4.1 有限元模型及參數

PBG大壩壩體斷面尺寸和材料分區如圖2所示，選取大壩滲流監測河床0+240斷面，建立壩-地基非飽和滲流二維有限元模型，如圖3所示。壩體最大壩高186m，壩頂寬度14m，上游壩坡比為1∶2和1∶2.25，下游壩坡比為1∶1.8。心墻頂寬4m，底寬96m，上、下游坡比均為1∶0.25；在心墻上、下游側設置雙層的反濾層，反濾層與堆石間設過渡層。壩基下有3層覆蓋層，最大厚度75.4m，具有中等-強透水性。壩基設主、副兩道各厚1.2m的混凝土防滲墻，墻中心間距14m，墻底嵌入基巖1.5m。主防滲墻底帷幕深度為50m；副墻在主墻上游側，墻頂插入心墻內部10m，墻底帷幕深度為10m。其地基模擬范圍為：上下游方向自坡腳延伸2倍壩高，深度方向延伸2倍壩高。離散后的壩-地基體系有限元模型含11555個單元，其中壩體5912個單元，防滲墻及墻底帷幕102個單元，壩基覆蓋層①816個單元，覆蓋層②1311個單元，基巖3414個單元。壩體壩基各類材料的滲透系數見表1。根據PBG工程心墻滲壓計的布置和運行情況，選取0+240斷面測點分析，測點布置如圖4所示，其中正常在測測點有P12、P13、P15、P16、P19、P20、P21。

4.2 非飽和滲流參數反演

通過敏感性分析，結合上述提到的3種非飽和滲流模型，確定需要反演的參數有心墻上區滲透系數、心墻下區滲透系數、弱風化基巖滲透系數以及非飽和特性參數。輸入2014年PBG工程運行時庫水位歷時過程線，每2個月的一天作為時間控制點，共6個控制節點。3種模型的參數分別進行正交試驗設計，輸入有限元中計算得到各組正交試驗的滲壓監測序列。采用不含交叉項的3次多項式，各時間控制節點試驗滲壓進行響應面方程擬合。各測點基于滲壓實測值和響應面擬合值的誤差平方和構建目標函數，采用多目標優化對多個目標函數求解，篩選得到合理反演參數。將反演參數代回有限元模型進行非飽和非穩定滲流計算，同時不考慮非飽和特性進行非穩定滲流計算。反分析誤差分析結果見表2，擬合效果選取心墻上區測點P12和心墻下區測點P15、P16進行展示，如圖5—7所示。

圖2 PBG工程河床0+240斷面圖

圖3 PBG0+240斷面各分區有限元網格

表1 壩體及壩基各區材料滲透系數 單位：cm/s

圖4 PBG堆石壩心墻特征點示意圖

對比非飽和非穩定滲流與飽和非穩定滲流，非飽和特性對心墻下區滲壓影響較小，2種滲流的心墻下區測點的平均絕對誤差分別為1.18m(VG模型)和1.97m，但是對心墻非飽和區滲壓影響較大，2種滲流的平均絕對誤差分別為3.53m(VG-G模型)和8.75m，飽和非穩定滲流的最大絕對誤差達到35.50m。因此非飽和滲流有限元模型不僅能反映飽和區滲流變化，也能反映非飽和區滲流變化，使心墻滲流的模擬更符合實際情況。同時在3種模型的對比中，對于心墻上區測點，VG模型、VG-G模型和BC模型的平均絕對誤差分別為2.83、3.53、2.53m，對于心墻下區測點，VG模型、VG-G模型和BC模型的平均絕對誤差分別為1.18、0.98、0.90m，因此BC模型對心墻飽和區和非飽和區的滲壓變化的擬合情況更加準確。

4.3 模型識別

基于3種非飽和滲流模型計算的擬合滲壓值和心墻實際滲壓值，分別采用AIC和BIC準則對VG模型、BC模型、VG-G模型3種非飽和滲流模型進行識別，AIC和BIC準則計算結果見表3，如圖8所示。VG模型、VG-G模型、BC模型的AIC均值分別為604.48、522.96、374.06，BIC均值分別為619.40、543.84、388.97，2種準則結果一致，均表明3種模型中BC模型最優，VG-G模型稍差，VG模型最差。

表2 反演結果誤差分析 單位：m

圖5 P12 VG、VG-G和BC模型反分析及非穩定滲壓對比

圖6 P15 VG、VG-G和BC模型反分析及非穩定滲壓對比

單從非飽和區測點P12分析，VG-G模型的AIC值為689.09，大于BC模型的298.98，小于VG模型的787.32，表明對于非飽和區，BC模型最優，VG-G模型稍差，VG模型最差。飽和區測點模型識別結果表明，BC模型和VG-G模型的AIC和BIC值均接近，平均相差28.89，VG模型更劣于VG-G模型，平均相差73.17。總而言之，3種模型對心墻飽和區的擬合差異遠小于心墻非飽和區，BC模型對于心墻滲壓變化的擬合情況更加準確，更能整體反映大壩心墻滲壓場的變化情況。

圖7 P16 VG、VG-G和BC模型反分析及非穩定滲壓對比

圖8 非飽和滲流模型識別結果圖

表3 非飽和滲流模型識別結果

5 結語

本文以PBG礫石土心墻堆石壩0+240斷面為基礎，以SWCC和滲透系數函數為基礎，構建了3種非飽和非穩定滲流模型，結合實測滲壓值和反演參數代入有限元計算得到模型擬合滲壓值使用AIC、BIC準則對非飽和滲流模型進行最優模型識別。結果表明，3種非飽和滲流模型的擬合滲壓值與實測值誤差較小，可以反映心墻滲壓變化；同時在這3種模型中，反演結果誤差分析與AIC、BIC準則計算結果一致，均表明BC模型最為優越，對于心墻滲壓變化擬合效果更加準確，更能整體反映堆石壩心墻滲壓場的實際變化情況。