基于響應面法的堆石體流變模型參數反演

解楓贊，包騰飛，李子陽，張玉霞，任 杰

（1.河南省河口村水庫運行中心，河南 濟源 454650；2.河海大學 水利水電學院，江蘇 南京 210098；3.南京水利科學研究院 水災害防御全國重點實驗室，江蘇 南京 210029）

堆石體變形除受材料應力狀態的影響外，還受時間的影響。大量工程監測資料分析表明，混凝土面板堆石壩在蓄水運行后一定時期內會發生較大變形，說明堆石體具有流變特性［1］。堆石體的流變變形對面板及趾板等防滲體結構安全造成很大的影響［2］。例如：羅馬尼亞里蘇混凝土面板堆石壩（壩高為60 m）的流變變形導致面板產生大量裂縫［3］；天生橋混凝土面板堆石壩的流變變形造成面板脫空、產生裂縫［4］。因此，有必要研究堆石體的流變特性并預測其時空變化規律。堆石體流變變形的預測精度取決于流變參數的準確性，流變參數主要通過工程類比或室內試驗獲得［5］，但工程類比法精度低，室內試驗只能采用縮尺三軸試驗，無法準確獲得原級配堆石體的流變參數。

采用實測變形資料反演堆石體流變參數是獲取真實材料力學參數的主要方法之一。反演參數方法有單純形法［6］、復合形法［7］、遺傳算法［8］、蟻群算法［9］等。若待反演參數較多，采用上述反演參數方法時需要進行大量的有限元分析計算，費時費力。響應面法采用響應面模型代替有限元模型，可以大大提高計算效率。響應面法是一種數學統計方法，能夠處理復雜系統的輸入（基本變量）和輸出（系統響應）的轉換關系。當系統輸入與輸出之間為隱式函數關系時，可以通過響應面法建立兩者間的顯式函數關系［10－11］。本文利用響應面法建立流變參數與面板堆石壩運行期流變位移之間的響應關系，利用有限元正分析計算結果求解響應面函數的系數，建立流變模型參數與流變位移間的顯式函數關系，基于面板堆石壩變形監測資料，建立目標尋優函數并求解，從而反演流變參數。

1 堆石體流變模型及響應面函數構建

當前，指數模型、九參數模型、Kelvin 模型、Maxwell 模型、Burgers 模型等流變模型常被用來描述堆石體材料的流變本構關系［12］。本文采用Burgers 模型，Burgers 模型由Kelvin 模型和Maxwell 模型組成，其形式為

響應面法處理復雜系統的輸入與輸出的轉換關系，實際上是構造一個包含未知參數的多項式函數，以此替代隱式函數，借助統計方法求解最優的響應面函數［13－14］。因此，應用響應面法最關鍵的步驟是構造高精度、高效率的響應面函數。

綜合考慮計算精度和計算量等，對于n個隨機變量x1，x2，…，xn，通常采用不含交叉項的二次多項式來描述系統輸入與輸出間的響應關系，公式為

式中：s（X）為系統輸出；X為系統輸入，X＝{x1，x2，…，xn} ；a、bi、ci為待定系數。

對于堆石體流變參數反演分析而言，主要是用響應面函數建立流變參數與大壩流變變形間的非線性函數關系。參考堆石體Burgers 流變模型和響應面法的基本原理，構造含有指數項的堆石體流變變形響應面函數，表達式為

式中：sk（X）為測點k豎向位移的流變分量，xi為待反演參數。

然而，上述各反演參數的數量級相差較大，主堆石體材料與次堆石體材料的參數差異也很顯著，如果直接將上述參數值代入響應面函數，必定會造成函數計算誤差較大。因此，有必要對參數進行無量綱化處理，即將每個模型參數都與其初始值相除。需要指出的是，用于有限元計算的參數是未經無量綱化處理的，參數無量綱化只在求解響應面函數時進行，對最優反演結果進行無量綱化逆操作，即可得到真正的反演值［15］。

參數無量綱化處理公式為

綜上，堆石體流變變形響應面函數形式如下：

2 基于響應面法的Burgers 模型參數反演分析

2.1 響應面函數求解

以測點P的流變變形為例，詳細闡述響應方程待定系數的求解步驟。

1）對模型待反演參數進行組合，基于參數組合對各參數加、減擾動量，以此增加組合數量，共構造13 個參數組合。參數組合公式為

2）利用有限元法計算各組合中壩體的流變變形，將流變變形結果作為式（9）等號的左邊項，形成測點P的響應面方程組。

3）若待定系數的個數為Q，那么其與待反演參數個數n的關系為Q＝2n＋1。采用上述參數組合方法恰好可形成2n＋1 個方程，則待定系數個數和線性方程個數相等，基于線性方程組唯一解的原則，求解得到測點P的唯一響應面函數。同理，可得到其他測點的響應面函數。

2.2 堆石體流變參數反演分析

確定各測點的響應面函數后，采用最優化方法確定響應面函數系數，目標尋優函數為流變分量計算值與實測值的均方根誤差最小：

采用Levenberg－Marquart（L－M）算法求解目標函數的最小值，進而反演得出堆石體流變模型參數。L－M算法是介于梯度下降法與高斯－牛頓法之間的非線性優化算法，既具有梯度下降法的全局尋優特性又具有高斯－牛頓法的局部收斂特性，能夠處理冗余參數問題［16］，對過參數化問題敏感性低，迭代收斂速度快，不易陷入局部最優。計算公式如下：

將反演得到的參數作為有限元模型參數，再次計算流變變形，比較流變變形計算值與實測值，若兩者相差較大，則以反演參數為基準值重新進行模型參數組合、響應面方程求解和參數反演等，直到目標函數滿足精度要求。

總結反演分析基本流程如下：1）采用有限元法計算面板堆石壩流變變形［17］，即求解流變位移的正分析過程；2）分析面板堆石壩運行期的流變位移數據，選擇用于反分析的測點和時間段；3）基于Burgers 流變模型的數學形式，構造響應面函數；4）設計模型參數組合并進行有限元正分析，采用組合參數和計算的流變位移求解響應面函數系數；5）構造測試樣本參數組，檢驗響應面函數的計算精度；6）定義參數反演的目標函數，采用L－M 算法搜索最優解，確定反演參數；7）由反演參數進行有限元分析得到特征點流變位移，與實測值比較，若兩者相差較大，則以反演參數為基準值重復進行步驟3）～步驟6），直到目標函數滿足精度要求。

3 工程實例

3.1 工程概況

某混凝土面板堆石壩由面板、墊層、過渡層、主堆石區、次堆石區、下游塊石護面及壩前輔助防滲料區等組成。最大壩高120.0 m，壩頂高程760.0 m，壩頂長259.8 m、寬9.0 m，上下游壩坡分別為1 ∶1.40、1 ∶1.35。面板頂部厚0.3 m、底部厚0.7 m。水庫校核洪水位759.1 m，設計洪水位756.2 m，正常蓄水位755.0 m，總庫容2.65 億m3，為多年調節水庫。堆石壩于1998 年1 月13 日開始填筑，至2000 年9 月填筑至760.0 m，澆筑基本完成。水庫于1999 年10 月11 日下閘蓄水，至2000 年9 月24 日庫水位為745.08 m。

3.2 有限元模型建立

依據工程資料建立堆石壩三維有限元模型（見圖1），模型包括面板、墊層、過渡層、主堆石區和次堆石區，共有17 854 個單元、19 545 個節點。

圖1 某混凝土面板堆石壩三維有限元模型

3.3 特征點選取

在該混凝土面板堆石壩布置水平和垂直位移監測儀器，利用垂直位移進行反演分析。壩體垂直位移監測儀器為安裝在最大斷面的NSC－1 型水管式沉降儀，測點布置見圖2。選取W2、W3、W6、W7、W10、W11 測點作為反演分析的特征點。這6 個測點位于壩高1／3～2／3 處，其沉降位移較大，能很好地反映堆石壩的沉降變形特征。

圖2 壩體監測儀器布置

3.4 流變模型參數組合

流變模型參數組合見表1。組合1 為堆石體參數初始值，采用式（9）參數組合方法生成另外12 個參數組合。

表1 流變模型參數組合

3.5 有限元正演計算

該混凝土面板堆石壩于2000 年9 月24 日蓄水至正常水位745.08 m，在此之后水位維持在正常水位附近。本文選取2000 年9 月24 日至2004 年9 月17 日作為流變參數反演時段，該時段內起止2 d 的水位相同，2 d 的壩體豎向位移增量可認為是面板堆石壩4 a間豎向位移的流變變形。W2、W3、W6、W7、W10、W11 測點的實測流變分量分別為5.89、6.13、5.48、5.73、3.77、3.95 cm。利用表1 的13 個參數組合進行有限元分析，得到6 個測點在13 個參數組合的流變位移。

3.6 各測點響應面函數求解

按照式（6）和式（7）對模型參數進行無量綱化處理，將無量綱化后的模型參數和其對應的流變位移代入已定義好的響應面方程，構建方程組并求解待定系數。以測點W2 為例，W2 測點的響應面函數系數a、bi、ci組成待求矩陣A1×13，無量綱化后的模型參數構成矩陣M13×13，不同參數組合下的W2 測點流變位移計算值構成矩陣S1×13，則有：

由式（12）得到W2 測點的響應面函數系數，其響應面函數表達式見式（13）。同理，可得其他測點的響應面函數。

3.7 響應面函數測試

為了檢驗響應面函數是否能準確擬合Burgers 模型參數與堆石體流變變形間的非線性關系并且達到代替有限元數值計算的精度，構造測試樣本參數組，利用構造好的響應面函數計算流變位移，同時采用有限元法計算流變位移，將響應面法和有限元法計算的流變位移進行對比分析，據此判斷響應面法的精度。

采用正交試驗設計方法［18］構造測試樣本，以6 個典型測點的流變位移為正交試驗指標，主堆石體和次堆石體的6 個流變參數為試驗因素，各參數均采用2水平進行試驗，依據表1 中的堆石體參數初始值設計因素水平，見表2。依據設計要求選用正交設計表L8（26）構造8 個測試樣本參數組，見表3。

表2 試驗因素水平

表3 測試樣本參數

采用表3 中的參數組合和有限元法計算壩體的流變位移，提取典型測點的流變位移值。利用構建的響應面函數計算壩體的流變位移，對表3 中的參數進行無量綱化處理，初始值仍取表1 中的參數組合1，代入構建好的響應面函數，得到典型測點的流變位移值，見表4。

表4 響應面函數流變位移計算結果 cm

比較有限元法計算位移和響應面函數解析的流變位移值可知，兩者的相對誤差均小于1.5%，表明利用響應面函數模擬面板堆石壩有限元法計算過程的精度較高，在反演分析中使用響應面函數代替有限元法計算是可行的，同時證明本文在構建響應面函數、設置流變模型參數組合、選取測點等方面是合理的。

3.8 流變參數反演

將得到的響應面函數和面板堆石壩6 個測點的實測流變位移代入目標函數，求解目標函數最小值。使用L－M 算法程序進行尋優，得到模型參數的反演值，見表5。

表5 流變參數反演結果

3.9 算法評價

在響應面法反演模型參數過程中，為反演6 個未知參數，需要進行13 次有限元正分析，而采用常規的間接法反演參數一般需要進行70 余次有限元計算［19］。顯然，使用響應面法進行參數反演分析具有更高的計算效率。

根據反演的流變參數進行有限元正分析，選取2001—2005 年每年9 月正常高水位的時間點作為有限元流變計算的時間點，計算各測點在典型時間節點的流變位移，其中以各典型時間節點相對于2000 年9 月24 日的豎向位移值為流變位移，結果見表6。

表6 流變位移計算值與實測值比較

由表6 可以看出，采用反演的參數、有限元法計算得到的特征點流變位移與實測值接近，誤差較小，最大相對誤差為5.64%，說明響應面法的反演參數精度較高。

4 結束語

本文采用響應面法和L－M 算法進行Burgers 模型參數反演分析，通過工程算例進行驗證，得出如下結論：

1）采用響應面法建立堆石體Burgers 流變模型參數與面板堆石壩運行期流變位移間的響應關系，求解得到響應面函數，通過構造測試樣本參數組檢驗可得，響應面函數的計算精度滿足要求，在反演分析中使用響應面函數代替有限元法計算是可行的。

2）使用L－M 算法搜索目標函數最小值，進而求解最優參數組合，通過算例分析表明響應面法的反演參數精度較高，說明響應面法應用于面板堆石壩的堆石體流變參數反演分析是可行的。