許云鵬，吳震宇，尹 川

(四川大學(xué)水利水電學(xué)院 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065)

1 概述

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

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

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

2 非飽和滲流模型

達(dá)西定律被證明既適合飽和土體，也適合非飽和土體。與飽和土體相比，非飽和土體滲透系數(shù)不是一個(gè)常數(shù)，它隨著土體飽和度或基質(zhì)吸力的變化而變化。根據(jù)達(dá)西定律及流體運(yùn)動(dòng)連續(xù)性方程可導(dǎo)出二維非穩(wěn)定飽和—非飽和滲流控制方程為：

(1)

(2)

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

對于穩(wěn)定流條件，基質(zhì)吸力剖面的時(shí)間依賴性消失，等式(1)的右側(cè)變?yōu)榱?。以下介紹了本工作中用于描述飽和和非飽和流動(dòng)的主要本構(gòu)模型和參數(shù)。Brooks和Corey基于土體含水量試驗(yàn)提出了有效飽和度和基質(zhì)吸力之間關(guān)系的土水特征曲線模型：

(3)

式中，h—為吸力水頭；hb—進(jìn)氣壓力水頭。

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

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

(4)

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

滲透系數(shù)隨基質(zhì)吸力變化而變化。BC模型和VG模型分別給出了滲透系數(shù)函數(shù)關(guān)系方程如下：

(5)

(6)

式中，Ks—飽和土滲透系數(shù)。

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

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

(7)

式中，α—滲透系數(shù)開始下降的起點(diǎn)；β—滲透系數(shù)隨基質(zhì)吸力下降的斜率。

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

3 模型識別方法

3.1 模型識別準(zhǔn)則

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

AIC=2k-2lnL

(8)

式中，k—參數(shù)個(gè)數(shù)，L—極大似然函數(shù)。

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

BIC=klnN-2lnL

(9)

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

BIC準(zhǔn)則考慮了樣本數(shù)量，避免精度過高選擇復(fù)雜模型。使用AIC、BIC準(zhǔn)則進(jìn)行模型識別時(shí)，值越小模型越優(yōu)。

3.2 模型識別方法流程

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

圖1 模型識別方法

4 工程實(shí)例

4.1 有限元模型及參數(shù)

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

4.2 非飽和滲流參數(shù)反演

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

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

圖3 PBG0+240斷面各分區(qū)有限元網(wǎng)格

表1 壩體及壩基各區(qū)材料滲透系數(shù) 單位：cm/s

圖4 PBG堆石壩心墻特征點(diǎn)示意圖

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

4.3 模型識別

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

表2 反演結(jié)果誤差分析 單位：m

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

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

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

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

圖8 非飽和滲流模型識別結(jié)果圖

表3 非飽和滲流模型識別結(jié)果

5 結(jié)語

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