王一凡，徐毓威，曾 俊，于 浩，潘中祺

（1.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢430072；2.武漢大學(xué)水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，武漢430072；3.長江勘測規(guī)劃設(shè)計(jì)研究有限責(zé)任公司，武漢430010）

0 引 言

裂隙巖體的滲透系數(shù)是研究地下水運(yùn)動(dòng)［1］、污染物遷移［2］的重要水文地質(zhì)參數(shù)，也是滲流分析與計(jì)算中的關(guān)鍵參數(shù)，其取值對(duì)于水利水電工程的設(shè)計(jì)與安全評(píng)價(jià)具有重要影響［3，4］。在工程實(shí)際中，巖體的滲透特性參數(shù)通常采用現(xiàn)場水文地質(zhì)試驗(yàn)方法來獲取，如壓水試驗(yàn)、抽水試驗(yàn)、微水試驗(yàn)等，但由于技術(shù)和經(jīng)濟(jì)條件的制約，通常僅在某些重點(diǎn)部位開展，因此研究滲流過程中巖體滲透系數(shù)的確定方法具有重要意義。

在工程施工與運(yùn)行過程中，通常會(huì)積累較為豐富的滲流與滲壓監(jiān)測資料，基于滲流、滲壓資料開展?jié)B流場的反饋分析研究已較為廣泛［5］，例如李守巨等［6］在滲壓和滲流量觀測值基礎(chǔ)上，采用Levenberg-Marquardt 優(yōu)化BP 神經(jīng)網(wǎng)絡(luò)，反演得到了工程樞紐區(qū)巖體及混凝土帷幕的滲透系數(shù)取值；劉武等［7］基于正交設(shè)計(jì)、有限元正分析、BP神經(jīng)網(wǎng)絡(luò)、遺傳算法建立起了多目標(biāo)非穩(wěn)定滲流反饋分析方法，解決了長河壩水電站基坑涌水問題。

Pest（Parameter Estimation）是由澳大利亞Water Numerical Computing咨詢公司John Doherty開發(fā)的獨(dú)立于模型的非線性參數(shù)估計(jì)與不確定性分析程序［8，9］，近年來在地下水與土壤水運(yùn)動(dòng)模型中應(yīng)用較為廣泛。Elcl 等［10］使用PEST 程序確定了地下水模型（MODFLOW）的水力學(xué)參數(shù)，研究了土耳其伊茲密爾地區(qū)地下水污染物的擴(kuò)散規(guī)律；Fang［11］等在Root Zone Water Quality Model（RZWQM）模型中估算了土壤水力學(xué)參數(shù)，發(fā)現(xiàn)該程序在參數(shù)尋優(yōu)方面具有高的效率；胡丹等［12］在MODFLOW-HY‐DRUS 耦合模型中對(duì)水文地質(zhì)參數(shù)和土壤水力學(xué)參數(shù)進(jìn)行優(yōu)化，顯著提高了模擬精度，但是該程序在水電工程巖體滲流分析中的應(yīng)用較少。

本文根據(jù)在工程實(shí)際中兩類常用的滲流場監(jiān)測數(shù)據(jù)：滲流量與滲壓構(gòu)建目標(biāo)函數(shù)，采用PEST 程序結(jié)合含自適應(yīng)罰函數(shù)的Signorini 變分不等式滲流分析方法［13］，建立起了一種全新的滲流場反饋分析方法，并通過一個(gè)簡單算例，驗(yàn)證了該方法的可靠性，并成功地將該方法應(yīng)用于我國西南某水電站巖體穩(wěn)態(tài)滲流場分析，取得了良好的擬合效果。

1 穩(wěn)定滲流理論及其反問題

1.1 穩(wěn)定滲流分析基本方程

在水利水電滲流計(jì)算中，由于巖體裂隙網(wǎng)絡(luò)的復(fù)雜性，通常采用連續(xù)介質(zhì)力學(xué)方法。在不考慮巖土體變形以及流體的可壓縮性情況下，流體的連續(xù)性方程可表示為：

考慮巖土體中流體運(yùn)動(dòng)服從達(dá)西定律，結(jié)合式（1）可得到巖土體滲流控制方程：

式中：h=h（x，y，z）為代求水頭函數(shù)；Kx、Ky、Kz分別為x、y、z方向的滲透系數(shù)。

基于式（2）的控制方程求解實(shí)際滲流問題時(shí)，通常會(huì)涉及以下4類邊界條件。

（1）水頭邊界條件。

（2）流量邊界條件。

（3）出滲面Signorini型互補(bǔ)邊界條件。

式中：Γs為潛在出滲邊界。

（4）自由面邊界條件。

式中：Γf={(x，y，z)|?=z}為自由面，即濕區(qū)與干區(qū)的分界面。

1.2 穩(wěn)定滲流反問題表述

假設(shè)在穩(wěn)定滲流計(jì)算中需要確定的巖體滲透系數(shù)有X組，而通過監(jiān)測儀器獲得的滲壓和滲流量數(shù)據(jù)分別為M、N組，那么穩(wěn)定滲流的反問題可以表述為：在待確定的巖體滲透系數(shù)的范圍內(nèi)，尋找一組最優(yōu)的組合，使得由滲流量和滲壓值構(gòu)成的如下目標(biāo)函數(shù)最?。?/p>

式中：Φ為目標(biāo)函數(shù)；Hi和分別為第i支滲壓計(jì)測點(diǎn)的計(jì)算水頭與監(jiān)測水頭；Qj和分布為第j條排水廊道測點(diǎn)的計(jì)算流量與監(jiān)測流量；wH與wQ分別表示滲壓權(quán)重系數(shù)與流量權(quán)重系數(shù)。

2 PEST參數(shù)優(yōu)化算法

PEST 程序是介于高斯牛頓法與梯度下降法之間的一種非線性優(yōu)化方法，同時(shí)具備高斯牛頓法的快速收斂性和梯度下降法的全局搜索性，能夠在多維的參數(shù)空間內(nèi)優(yōu)化模型輸入?yún)?shù)。該程序主要基于高斯-馬奎特-列文伯格（Gauss-Marquart-Levenberg）算法，對(duì)模型的輸入?yún)?shù)進(jìn)行優(yōu)化調(diào)整，進(jìn)而求取目標(biāo)函數(shù)（模型計(jì)算值與實(shí)際觀測值的差異函數(shù)）的最小值。具體過程如下：

假定n個(gè)模型輸入?yún)?shù)存儲(chǔ)于向量x中，m個(gè)模擬值存儲(chǔ)于向量y中，則向量x與向量y的關(guān)系為：

式中：M為模型輸入?yún)?shù)與模擬值相關(guān)聯(lián)的非線性函數(shù)。

將y=M（x）在參數(shù)初值及其計(jì)算值（x0，y0）處通過一階泰勒公式展開得到：

帶入雅可比矩陣表示為：

式中：矩陣J是y=M（x）在x0處的雅可比矩陣，大小為m行n列，記錄了n個(gè)模擬值分別對(duì)m個(gè)模型輸入?yún)?shù)的偏導(dǎo)數(shù)。

定義參數(shù)反演的非線性目標(biāo)函數(shù)為：

式中：y1是模擬值y對(duì)應(yīng)的真實(shí)值；w為相應(yīng)的權(quán)重向量。

由于對(duì)（8）式以一階泰勒公式展開，要求x-x0足夠小，所以參數(shù)優(yōu)化過程需要不斷更新參數(shù)向量x和雅可比矩陣J，使得目標(biāo)函數(shù)最小，定義參數(shù)更新向量：u=x-x0，結(jié)合（11）式可表示為：

3 算例驗(yàn)證

為了驗(yàn)證PEST 反演算法的有效性與可行性，本文參考武曉煒等［14］使用GA-BP算法開展的研究論文，建立了同樣的土石壩二維模型［14］，模型共有單元5 768個(gè)，節(jié)點(diǎn)9 459 個(gè)，上游水位40 m，如圖1所示。其中，I 為壩體、Ⅱ?yàn)榉罎B帷幕、Ⅲ為第一層基巖、Ⅳ為第二層基巖，依據(jù)工程類比設(shè)置4類材料的滲透系數(shù)理論值分別為2.0×10-3m/s、6.0×10-7m/s、6.5×10-4m/s、1.0×10-6m/s，參數(shù)反演變化范圍以理論值為中心，上下跨越兩個(gè)數(shù)量級(jí)。模型內(nèi)總共布置了3個(gè)滲壓計(jì)P1、P2、P3和一個(gè)量水堰Q1。

為了確定4類材料的滲透系數(shù)，本文分別采用了PEST程序和GA-BP 算法進(jìn)行反演分析。其中，在反演過程中，采用GABP 算法需要對(duì)每個(gè)參數(shù)設(shè)計(jì)7 個(gè)水平的49 組正交試驗(yàn)和10 組均勻試驗(yàn)，并基于參數(shù)輸入、滲流量與滲壓輸出數(shù)據(jù)，通過反復(fù)試算確定了4-9-34-4 的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，最終通過遺傳算法進(jìn)行優(yōu)化；而采用PEST 程序，只需直接調(diào)用模型計(jì)算63 次，就能完成參數(shù)優(yōu)化?；趦煞N算法獲得的參數(shù)優(yōu)化結(jié)果如表1所示，由表可以看出，在相同的參數(shù)變化范圍內(nèi)，采用PEST 法獲得的4 種材料滲透系數(shù)值同理論值的相對(duì)誤差分別為0.61%、0.78%、0.69%、0.69%，均不超過1%，并且小于GA-BP 法的相對(duì)誤差（2.00%、1.50%、6.15%、9.00%）。所以，同GA-BP 法相比，PEST 反演算法不僅在過程上直接簡便，而且具有更高的準(zhǔn)確度。

表1 待反演參數(shù)設(shè)置與優(yōu)化結(jié)果

4 工程實(shí)例

為了研究PEST 參數(shù)反演算法在大型水利水電工程滲流分析中的適用性，本文以我國西南某水電站為背景，對(duì)其2014-2015年的蓄水周期內(nèi)滲流計(jì)算過程中涉及的巖體滲透系數(shù)開展反演分析。

4.1 工程概況

該水電站位于金沙江上游高山峽谷段，谷底寬闊平緩，兩岸山體陡峭雄厚，河谷呈典型的對(duì)稱“U”型結(jié)構(gòu)。工程樞紐區(qū)由混凝土雙曲拱壩，泄洪消能設(shè)施及兩岸引水發(fā)電建筑物等組成，最大壩高285.5 m，正常蓄水位600 m，裝機(jī)容量12 600 MW，壩段控制流域面積占金沙江的96%，約為45.44 萬km2。

壩區(qū)河床基巖及兩岸邊坡巖性主要為二疊系上統(tǒng)峨眉山玄武巖，下部埋藏二疊系下統(tǒng)茅口組石灰?guī)r，玄武巖與灰?guī)r間夾有2～3 m厚的泥頁巖夾層。樞紐區(qū)主要地質(zhì)構(gòu)造為層間和層內(nèi)錯(cuò)動(dòng)帶，整體上巖體性質(zhì)較好，左右岸各發(fā)育有12 條層間錯(cuò)動(dòng)帶，其平均厚度約為0.1～0.6 m。依據(jù)鉆孔壓水?dāng)?shù)據(jù)，并考慮巖體的風(fēng)化程度隨埋深增加而減弱，可人為地將玄武巖劃分為中等透水，弱透水上段、弱透水下段、微微透水、新鮮巖體等5個(gè)滲透分區(qū)，如圖2所示。

4.2 有限元模型

根據(jù)該水電站樞紐區(qū)的工程地質(zhì)條件、水文地質(zhì)條件和樞紐布置情況，考慮到工程左右岸屬于獨(dú)立的地下水流動(dòng)系統(tǒng)，本文建立了該水電站右岸三維有限元模型如圖3所示，模型共劃分實(shí)體單元433 萬個(gè)，節(jié)點(diǎn)134 萬個(gè)，上下游相距3 220 m，右邊界距河床中心線約1 775 m，底部高程-180 m。同時(shí)，采用排水子結(jié)構(gòu)技術(shù)［15］對(duì)廠壩區(qū)防滲排水系統(tǒng)進(jìn)行了精細(xì)模擬，如圖4、圖5所示。為了開展廠壩區(qū)三維滲流場計(jì)算，設(shè)置模型的邊界條件具體為：上下游庫水和河水淹沒范圍內(nèi)的節(jié)點(diǎn)取定水頭邊界；右岸山體側(cè)邊界與下游側(cè)邊界取隔水邊界；將與底層廊道相連向上匯水的溢流型排水孔設(shè)置為定水頭邊界；其他排水廊道內(nèi)排水孔設(shè)置為Signorini型潛在溢出邊界。

4.3 反演計(jì)算

考慮到樞紐區(qū)巖體的風(fēng)化程度差異，結(jié)合地質(zhì)勘測資料與工程施工情況，本研究設(shè)置了強(qiáng)風(fēng)化，弱風(fēng)化上段、弱風(fēng)化下段、微風(fēng)化、新鮮巖體、灰?guī)r、泥頁巖夾層、混凝土材料、防滲帷幕、廠房開挖擾動(dòng)區(qū)等共計(jì)13 個(gè)待反演滲透系數(shù)，并依據(jù)有關(guān)鉆孔壓水試驗(yàn)結(jié)果與工程經(jīng)驗(yàn)合理設(shè)置待反演滲透系數(shù)上下限、各向異性比值，同時(shí)選取廠壩區(qū)具有代表性的排水廊道6條（廠區(qū)5條、壩區(qū)1條）、滲壓計(jì)14支（廠區(qū)6支，壩區(qū)8支）分別作為滲流量監(jiān)測點(diǎn)和滲壓監(jiān)測點(diǎn)，其空間布置如圖4所示，并由這些監(jiān)測點(diǎn)數(shù)據(jù)共同構(gòu)成參數(shù)反演目標(biāo)函數(shù)。

由于水庫的蓄水過程會(huì)導(dǎo)致區(qū)域水文地質(zhì)條件的顯著改變，而巖體的滲透系數(shù)可能會(huì)受到水力耦合效應(yīng)的影響而產(chǎn)生一定程度的變化，因此本文僅針對(duì)該水電站的2014-2015年的蓄水周期內(nèi)滲流場特征，以上游庫水位在600 m 正常蓄水位時(shí)作為計(jì)算工況，利用PEST 程序?qū)υO(shè)定的13 個(gè)待反演參數(shù)進(jìn)行優(yōu)化，參數(shù)設(shè)置與優(yōu)化結(jié)果如表2所示。

表2 待反演參數(shù)設(shè)置與優(yōu)化結(jié)果

反演過程中目標(biāo)函數(shù)變化曲線如圖6所示，該圖表明在總計(jì)7 次迭代，調(diào)用模型114 次過程中，目標(biāo)函數(shù)呈現(xiàn)逐次遞減的趨勢，并且具有快速收斂性。

依據(jù)PEST 程序參數(shù)反演獲取的滲透系數(shù)，本文分別對(duì)庫水位為550、560、570、580、590、600 m 的滲流場特征進(jìn)行了模擬。其中排水廊道內(nèi)流量測點(diǎn)計(jì)算值與實(shí)測值對(duì)比如圖7所示，結(jié)果表明計(jì)算值與監(jiān)測值總體吻合較好，其中廠區(qū)滲流量明顯高于壩區(qū)，并且同上游庫水位具有一定的相關(guān)性，但是由于部分量水堰的埋藏深度與位置的差異會(huì)使得流量監(jiān)測數(shù)據(jù)具有較大的離散性，導(dǎo)致計(jì)算值與監(jiān)測值的差別相對(duì)較大，如Q1、Q4等。

廠壩區(qū)代表性滲壓測點(diǎn)的計(jì)算值與實(shí)測值對(duì)比如圖8所示，二者同樣吻合較好，其中帷幕前滲壓計(jì)（P1）與上游水位呈明顯的相關(guān)性，且變化幅度較大，而帷幕后滲壓計(jì)的水頭由于受防滲帷幕和排水孔幕的控制，波動(dòng)較小，較為平穩(wěn)。

為了評(píng)價(jià)反演結(jié)果的準(zhǔn)確性與可靠性，本文選取了3 類廣泛用于時(shí)間序列預(yù)測模型結(jié)果檢驗(yàn)的誤差評(píng)估指標(biāo)：平均絕對(duì)誤差MAE、均方根誤差RMSE，平均絕對(duì)百分比誤差MAPE，其表達(dá)式如下所示：

式中：Pi是第i個(gè)工況下的模型預(yù)測值；Hi是第i個(gè)工況下的模型實(shí)測值，N計(jì)算工況總數(shù)。

利用上述的3 類誤差評(píng)價(jià)指標(biāo)，對(duì)研究工況下的滲流量測點(diǎn)與滲壓測點(diǎn)的計(jì)算結(jié)果展開了誤差分析，如圖9、圖10所示。結(jié)果表明，各工況下的滲壓測點(diǎn)的MAE、RMSE、MAPE三類指標(biāo)均不超過6 m、6 m、2%；流量測點(diǎn)的MAE、RMSE、MAPE三類指標(biāo)基本不超過2.5 L/s、2.5 L/s、25%，依據(jù)相關(guān)工程經(jīng)驗(yàn)，能夠認(rèn)為反演計(jì)算得到的參數(shù)是準(zhǔn)確和可靠的。

5 結(jié) 論

本文基于PEST 參數(shù)率定程序，結(jié)合含自適應(yīng)罰函數(shù)的Signorini 變分不等式滲流分析方法，依據(jù)特定時(shí)間序列內(nèi)的滲流量與滲壓監(jiān)測數(shù)據(jù)構(gòu)建的目標(biāo)函數(shù)，建立了一種對(duì)巖體滲流參數(shù)進(jìn)行反演分析的方法，并成功應(yīng)用于實(shí)際工程。并通過采用平均絕對(duì)誤差MAE、均方根誤差RMSE，平均絕對(duì)百分比誤差MAPE等三類誤差評(píng)價(jià)指標(biāo)對(duì)根據(jù)反演參數(shù)獲得的滲流量與滲壓計(jì)算值與監(jiān)測值偏差進(jìn)行分析，結(jié)果表明，流量測點(diǎn)的MAE、RMSE、MAPE三類指標(biāo)基本不超過2.5 L/s、2.5 L/s、25%；該工程滲壓測點(diǎn)的MAE、RMSE、MAPE三類指標(biāo)均不超過6 m、6 m、2%，進(jìn)而證明了本文提出的反演方法的可靠性與準(zhǔn)確性。 □