基于響應(yīng)面法的堆石體流變模型參數(shù)反演

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

（1.河南省河口村水庫運(yùn)行中心，河南 濟(jì)源 454650；2.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098；3.南京水利科學(xué)研究院 水災(zāi)害防御全國重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029）

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

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

1 堆石體流變模型及響應(yīng)面函數(shù)構(gòu)建

當(dāng)前，指數(shù)模型、九參數(shù)模型、Kelvin 模型、Maxwell 模型、Burgers 模型等流變模型常被用來描述堆石體材料的流變本構(gòu)關(guān)系［12］。本文采用Burgers 模型，Burgers 模型由Kelvin 模型和Maxwell 模型組成，其形式為

響應(yīng)面法處理復(fù)雜系統(tǒng)的輸入與輸出的轉(zhuǎn)換關(guān)系，實(shí)際上是構(gòu)造一個(gè)包含未知參數(shù)的多項(xiàng)式函數(shù)，以此替代隱式函數(shù)，借助統(tǒng)計(jì)方法求解最優(yōu)的響應(yīng)面函數(shù)［13－14］。因此，應(yīng)用響應(yīng)面法最關(guān)鍵的步驟是構(gòu)造高精度、高效率的響應(yīng)面函數(shù)。

綜合考慮計(jì)算精度和計(jì)算量等，對(duì)于n個(gè)隨機(jī)變量x1，x2，…，xn，通常采用不含交叉項(xiàng)的二次多項(xiàng)式來描述系統(tǒng)輸入與輸出間的響應(yīng)關(guān)系，公式為

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

對(duì)于堆石體流變參數(shù)反演分析而言，主要是用響應(yīng)面函數(shù)建立流變參數(shù)與大壩流變變形間的非線性函數(shù)關(guān)系。參考堆石體Burgers 流變模型和響應(yīng)面法的基本原理，構(gòu)造含有指數(shù)項(xiàng)的堆石體流變變形響應(yīng)面函數(shù)，表達(dá)式為

式中：sk（X）為測(cè)點(diǎn)k豎向位移的流變分量，xi為待反演參數(shù)。

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

參數(shù)無量綱化處理公式為

綜上，堆石體流變變形響應(yīng)面函數(shù)形式如下：

2 基于響應(yīng)面法的Burgers 模型參數(shù)反演分析

2.1 響應(yīng)面函數(shù)求解

以測(cè)點(diǎn)P的流變變形為例，詳細(xì)闡述響應(yīng)方程待定系數(shù)的求解步驟。

1）對(duì)模型待反演參數(shù)進(jìn)行組合，基于參數(shù)組合對(duì)各參數(shù)加、減擾動(dòng)量，以此增加組合數(shù)量，共構(gòu)造13 個(gè)參數(shù)組合。參數(shù)組合公式為

2）利用有限元法計(jì)算各組合中壩體的流變變形，將流變變形結(jié)果作為式（9）等號(hào)的左邊項(xiàng)，形成測(cè)點(diǎn)P的響應(yīng)面方程組。

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

2.2 堆石體流變參數(shù)反演分析

確定各測(cè)點(diǎn)的響應(yīng)面函數(shù)后，采用最優(yōu)化方法確定響應(yīng)面函數(shù)系數(shù)，目標(biāo)尋優(yōu)函數(shù)為流變分量計(jì)算值與實(shí)測(cè)值的均方根誤差最小：

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

將反演得到的參數(shù)作為有限元模型參數(shù)，再次計(jì)算流變變形，比較流變變形計(jì)算值與實(shí)測(cè)值，若兩者相差較大，則以反演參數(shù)為基準(zhǔn)值重新進(jìn)行模型參數(shù)組合、響應(yīng)面方程求解和參數(shù)反演等，直到目標(biāo)函數(shù)滿足精度要求。

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

3 工程實(shí)例

3.1 工程概況

某混凝土面板堆石壩由面板、墊層、過渡層、主堆石區(qū)、次堆石區(qū)、下游塊石護(hù)面及壩前輔助防滲料區(qū)等組成。最大壩高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，設(shè)計(jì)洪水位756.2 m，正常蓄水位755.0 m，總庫容2.65 億m3，為多年調(diào)節(jié)水庫。堆石壩于1998 年1 月13 日開始填筑，至2000 年9 月填筑至760.0 m，澆筑基本完成。水庫于1999 年10 月11 日下閘蓄水，至2000 年9 月24 日庫水位為745.08 m。

3.2 有限元模型建立

依據(jù)工程資料建立堆石壩三維有限元模型（見圖1），模型包括面板、墊層、過渡層、主堆石區(qū)和次堆石區(qū)，共有17 854 個(gè)單元、19 545 個(gè)節(jié)點(diǎn)。

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

3.3 特征點(diǎn)選取

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

圖2 壩體監(jiān)測(cè)儀器布置

3.4 流變模型參數(shù)組合

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

表1 流變模型參數(shù)組合

3.5 有限元正演計(jì)算

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

3.6 各測(cè)點(diǎn)響應(yīng)面函數(shù)求解

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

由式（12）得到W2 測(cè)點(diǎn)的響應(yīng)面函數(shù)系數(shù)，其響應(yīng)面函數(shù)表達(dá)式見式（13）。同理，可得其他測(cè)點(diǎn)的響應(yīng)面函數(shù)。

3.7 響應(yīng)面函數(shù)測(cè)試

為了檢驗(yàn)響應(yīng)面函數(shù)是否能準(zhǔn)確擬合Burgers 模型參數(shù)與堆石體流變變形間的非線性關(guān)系并且達(dá)到代替有限元數(shù)值計(jì)算的精度，構(gòu)造測(cè)試樣本參數(shù)組，利用構(gòu)造好的響應(yīng)面函數(shù)計(jì)算流變位移，同時(shí)采用有限元法計(jì)算流變位移，將響應(yīng)面法和有限元法計(jì)算的流變位移進(jìn)行對(duì)比分析，據(jù)此判斷響應(yīng)面法的精度。

采用正交試驗(yàn)設(shè)計(jì)方法［18］構(gòu)造測(cè)試樣本，以6 個(gè)典型測(cè)點(diǎn)的流變位移為正交試驗(yàn)指標(biāo)，主堆石體和次堆石體的6 個(gè)流變參數(shù)為試驗(yàn)因素，各參數(shù)均采用2水平進(jìn)行試驗(yàn)，依據(jù)表1 中的堆石體參數(shù)初始值設(shè)計(jì)因素水平，見表2。依據(jù)設(shè)計(jì)要求選用正交設(shè)計(jì)表L8（26）構(gòu)造8 個(gè)測(cè)試樣本參數(shù)組，見表3。

表2 試驗(yàn)因素水平

表3 測(cè)試樣本參數(shù)

采用表3 中的參數(shù)組合和有限元法計(jì)算壩體的流變位移，提取典型測(cè)點(diǎn)的流變位移值。利用構(gòu)建的響應(yīng)面函數(shù)計(jì)算壩體的流變位移，對(duì)表3 中的參數(shù)進(jìn)行無量綱化處理，初始值仍取表1 中的參數(shù)組合1，代入構(gòu)建好的響應(yīng)面函數(shù)，得到典型測(cè)點(diǎn)的流變位移值，見表4。

表4 響應(yīng)面函數(shù)流變位移計(jì)算結(jié)果 cm

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

3.8 流變參數(shù)反演

將得到的響應(yīng)面函數(shù)和面板堆石壩6 個(gè)測(cè)點(diǎn)的實(shí)測(cè)流變位移代入目標(biāo)函數(shù)，求解目標(biāo)函數(shù)最小值。使用L－M 算法程序進(jìn)行尋優(yōu)，得到模型參數(shù)的反演值，見表5。

表5 流變參數(shù)反演結(jié)果

3.9 算法評(píng)價(jià)

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

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

表6 流變位移計(jì)算值與實(shí)測(cè)值比較

由表6 可以看出，采用反演的參數(shù)、有限元法計(jì)算得到的特征點(diǎn)流變位移與實(shí)測(cè)值接近，誤差較小，最大相對(duì)誤差為5.64%，說明響應(yīng)面法的反演參數(shù)精度較高。

4 結(jié)束語

本文采用響應(yīng)面法和L－M 算法進(jìn)行Burgers 模型參數(shù)反演分析，通過工程算例進(jìn)行驗(yàn)證，得出如下結(jié)論：

1）采用響應(yīng)面法建立堆石體Burgers 流變模型參數(shù)與面板堆石壩運(yùn)行期流變位移間的響應(yīng)關(guān)系，求解得到響應(yīng)面函數(shù)，通過構(gòu)造測(cè)試樣本參數(shù)組檢驗(yàn)可得，響應(yīng)面函數(shù)的計(jì)算精度滿足要求，在反演分析中使用響應(yīng)面函數(shù)代替有限元法計(jì)算是可行的。

2）使用L－M 算法搜索目標(biāo)函數(shù)最小值，進(jìn)而求解最優(yōu)參數(shù)組合，通過算例分析表明響應(yīng)面法的反演參數(shù)精度較高，說明響應(yīng)面法應(yīng)用于面板堆石壩的堆石體流變參數(shù)反演分析是可行的。