基于PEST的高壩樞紐區巖體穩態滲流場反演分析

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

（1.武漢大學水資源與水電工程科學國家重點實驗室，武漢430072；2.武漢大學水工巖石力學教育部重點實驗室，武漢430072；3.長江勘測規劃設計研究有限責任公司，武漢430010）

0 引 言

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

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

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

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

1 穩定滲流理論及其反問題

1.1 穩定滲流分析基本方程

在水利水電滲流計算中，由于巖體裂隙網絡的復雜性，通常采用連續介質力學方法。在不考慮巖土體變形以及流體的可壓縮性情況下，流體的連續性方程可表示為：

考慮巖土體中流體運動服從達西定律，結合式（1）可得到巖土體滲流控制方程：

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

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

（1）水頭邊界條件。

（2）流量邊界條件。

（3）出滲面Signorini型互補邊界條件。

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

（4）自由面邊界條件。

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

1.2 穩定滲流反問題表述

假設在穩定滲流計算中需要確定的巖體滲透系數有X組，而通過監測儀器獲得的滲壓和滲流量數據分別為M、N組，那么穩定滲流的反問題可以表述為：在待確定的巖體滲透系數的范圍內，尋找一組最優的組合，使得由滲流量和滲壓值構成的如下目標函數最小：

式中：Φ為目標函數；Hi和分別為第i支滲壓計測點的計算水頭與監測水頭；Qj和分布為第j條排水廊道測點的計算流量與監測流量；wH與wQ分別表示滲壓權重系數與流量權重系數。

2 PEST參數優化算法

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

假定n個模型輸入參數存儲于向量x中，m個模擬值存儲于向量y中，則向量x與向量y的關系為：

式中：M為模型輸入參數與模擬值相關聯的非線性函數。

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

帶入雅可比矩陣表示為：

式中：矩陣J是y=M（x）在x0處的雅可比矩陣，大小為m行n列，記錄了n個模擬值分別對m個模型輸入參數的偏導數。

定義參數反演的非線性目標函數為：

式中：y1是模擬值y對應的真實值；w為相應的權重向量。

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

3 算例驗證

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

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

表1 待反演參數設置與優化結果

4 工程實例

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

4.1 工程概況

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

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

4.2 有限元模型

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

4.3 反演計算

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

由于水庫的蓄水過程會導致區域水文地質條件的顯著改變，而巖體的滲透系數可能會受到水力耦合效應的影響而產生一定程度的變化，因此本文僅針對該水電站的2014-2015年的蓄水周期內滲流場特征，以上游庫水位在600 m 正常蓄水位時作為計算工況，利用PEST 程序對設定的13 個待反演參數進行優化，參數設置與優化結果如表2所示。

表2 待反演參數設置與優化結果

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

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

廠壩區代表性滲壓測點的計算值與實測值對比如圖8所示，二者同樣吻合較好，其中帷幕前滲壓計（P1）與上游水位呈明顯的相關性，且變化幅度較大，而帷幕后滲壓計的水頭由于受防滲帷幕和排水孔幕的控制，波動較小，較為平穩。

為了評價反演結果的準確性與可靠性，本文選取了3 類廣泛用于時間序列預測模型結果檢驗的誤差評估指標：平均絕對誤差MAE、均方根誤差RMSE，平均絕對百分比誤差MAPE，其表達式如下所示：

式中：Pi是第i個工況下的模型預測值；Hi是第i個工況下的模型實測值，N計算工況總數。

利用上述的3 類誤差評價指標，對研究工況下的滲流量測點與滲壓測點的計算結果展開了誤差分析，如圖9、圖10所示。結果表明，各工況下的滲壓測點的MAE、RMSE、MAPE三類指標均不超過6 m、6 m、2%；流量測點的MAE、RMSE、MAPE三類指標基本不超過2.5 L/s、2.5 L/s、25%，依據相關工程經驗，能夠認為反演計算得到的參數是準確和可靠的。

5 結 論

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