程 井，李培聰，李同春，袁 平

(1. 河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098；2. 中冶長(zhǎng)天國(guó)際工程有限責(zé)任公司，湖南 長(zhǎng)沙 410007)

壩體、壩基物理力學(xué)參數(shù)的選定是壩工設(shè)計(jì)和分析壩體應(yīng)力、變形以及裂縫形成機(jī)理的基礎(chǔ)。傳統(tǒng)反演方法主要建立在確定性分析基礎(chǔ)上。然而，由于混凝土壩工作環(huán)境的復(fù)雜多變性，壩體變形是一個(gè)隨機(jī)變量，變形過程是一個(gè)隨機(jī)過程，觀測(cè)值只是一個(gè)樣本序列或樣本值的實(shí)現(xiàn)，即觀測(cè)的變形具有不確定性，只能從概率意義來研究變形的均值與方差[1]。因此，考慮不確定性的隨機(jī)反演方法更符合實(shí)際情況。

目前，混凝土重力壩主要通過對(duì)原型觀測(cè)資料分析，建立確定性數(shù)學(xué)模型，反求大壩的材料參數(shù)以及某些結(jié)構(gòu)特性[2]。近年來，部分學(xué)者開始考慮水工結(jié)構(gòu)反演中的不確定性問題。楊杰等[3]提出基于最大熵原理的貝葉斯不確定性反分析法；蘇懷智等[4]在傳統(tǒng)的確定性模型中引入?yún)^(qū)間數(shù)因子，考慮了不確定性因素對(duì)反演結(jié)果的影響；杜永峰等[5]基于模糊理論提出考慮測(cè)量數(shù)據(jù)不確定性的結(jié)構(gòu)物理參數(shù)識(shí)別方法；雷鵬等[6-7]進(jìn)一步研究了參數(shù)不確定性區(qū)間反演。貝葉斯方法是一種考慮不確定性因素的有效反演方法，馬爾可夫鏈蒙特卡羅方法(MCMC)是目前貝葉斯方法的標(biāo)準(zhǔn)抽樣方法，在巖土工程、水資源優(yōu)化以及土壤學(xué)等學(xué)科已有學(xué)者進(jìn)行相關(guān)研究[8-10]。然而，對(duì)于許多復(fù)雜問題，似然函數(shù)難以獲得，因此，Marjoram等[11]首次提出了無似然函數(shù)的馬爾可夫鏈蒙特卡羅方法，在此基礎(chǔ)上，Daniel等[12]通過調(diào)節(jié)馬爾可夫鏈的容差，提高該方法的計(jì)算效率，簡(jiǎn)化了許多復(fù)雜的工程問題。

本文基于Bayesian理論，考慮參數(shù)先驗(yàn)信息分布及觀測(cè)值的不確定性，采用無似然函數(shù)的馬爾可夫鏈蒙特卡羅（MCMC without likelihoods）方法進(jìn)行彈性模量的隨機(jī)反演，研究參數(shù)后驗(yàn)分布的統(tǒng)計(jì)特性，并分別對(duì)觀測(cè)值出現(xiàn)波動(dòng)情況（人為觀測(cè)誤差及擬合中出現(xiàn)誤差）與觀測(cè)值進(jìn)行比較，評(píng)價(jià)模型誤差和預(yù)測(cè)精度。

1 基于現(xiàn)代貝葉斯理論的參數(shù)隨機(jī)反演

貝葉斯原理中，令X為一連續(xù)型分布的隨機(jī)變量，其先驗(yàn)概率密度函數(shù)為π(X)，則后驗(yàn)概率密度函數(shù)f(X∣D)為：

式中：X為隨機(jī)變量（待反演參數(shù)）；D為觀測(cè)樣本；f (X∣D) 為隨機(jī)變量的后驗(yàn)分布；f (D∣X) 為似然函數(shù)，數(shù)值上等于D在X上的條件分布；π(X) 為隨機(jī)變量（參數(shù)X）的先驗(yàn)分布。

1.1 馬爾可夫鏈蒙特卡羅方法（MCMC）

對(duì)于復(fù)雜問題，若先驗(yàn)分布和后驗(yàn)分布不滿足共軛分布條件，常規(guī)貝葉斯方法難以得到后驗(yàn)分布估計(jì)值（均值、標(biāo)準(zhǔn)差）的解析式，因此，需要借助于數(shù)值方法或近似方法進(jìn)行模擬[8]。為解決該問題，采用馬爾可夫鏈蒙特卡羅方法計(jì)算后驗(yàn)分布。而MCMC方法的基本思想則是通過隨機(jī)抽樣建立一個(gè)最終平穩(wěn)分布為所求后驗(yàn)分布的馬爾可夫鏈，通過馬爾可夫鏈得到后驗(yàn)分布的樣本，進(jìn)而得到對(duì)應(yīng)的期望值和標(biāo)準(zhǔn)差。常用的MCMC方法有Gibbs抽樣、Metropolis抽樣等基本算法。采用的無似然函數(shù)MCMC算法具體步驟[11]如下：

(1) 初始化馬爾可夫鏈初始狀態(tài)X0=x0；

(2) 在i時(shí)刻（第i次循環(huán)），馬爾可夫鏈狀態(tài)Xi=xi，由轉(zhuǎn)移核q(Xi→Xi+1)采樣Xi+1；

(3)根據(jù)Xi+1及已知先驗(yàn)計(jì)算模型M，產(chǎn)生隨機(jī)樣本D′；

(4) 如果D′=D（已知抽樣樣本），進(jìn)入下一步，否則從第（2）步重新開始；

(5) 計(jì)算接受率h：

(6) 若滿足接收條件，則接收Xi+1，否則從第(2)步重新開始。

抽樣過程中，抽樣終止以馬爾可夫鏈?zhǔn)諗繛榍疤幔虼?，馬爾可夫鏈?zhǔn)諗颗c否對(duì)模型參數(shù)后驗(yàn)估計(jì)具有重要影響。BGR診斷法[13]是馬爾可夫鏈?zhǔn)諗啃栽\斷常用方法之一，該方法基于區(qū)間統(tǒng)計(jì)特性及長(zhǎng)度進(jìn)行診斷，定義診斷指標(biāo)為：

式中：L 為 馬爾可夫鏈總的序列區(qū)間長(zhǎng)度或統(tǒng)計(jì)特性；l為馬爾可夫鏈單鏈區(qū)間長(zhǎng)度均值或統(tǒng)計(jì)特性均值。若馬爾可夫鏈?zhǔn)諗浚祽?yīng)接近于1。

1.2 觀測(cè)數(shù)據(jù)（樣本）不確定性分析

壩體位移是一個(gè)隨機(jī)變量，變形過程是一個(gè)隨機(jī)過程，觀測(cè)值只是一個(gè)樣本序列或樣本值的實(shí)現(xiàn)，因此具有不確定性。根據(jù)文獻(xiàn)[2]大壩位移統(tǒng)計(jì)模型表達(dá)式為：

式中：δH，δT，δθ分別為水壓分量、溫度分量、時(shí)效分量?？紤]到混凝土壩工作環(huán)境的復(fù)雜多變性以及各種不確定因素（如測(cè)量誤差、儀器精度、擬合誤差等），假設(shè)對(duì)于某個(gè)區(qū)間內(nèi)的每個(gè)值有：

式（6）表明，大壩位移由兩部分組成，一部分是由水壓分量δH、溫度分量δT和時(shí)效分量δθ等組成，另一部分ε ～N(0, σ2)為不確定性因素造成的隨機(jī)誤差。則對(duì)于水壓分量有：

式中：H為上游水位測(cè)值與壩底高程之差；ai為水壓因子回歸系數(shù)；εH為水壓分量中不確定性因素造成的隨機(jī)誤差。

1.3 壩體壩基變形模量隨機(jī)反演

采用MCMC without likelihoods算法進(jìn)行參數(shù)隨機(jī)反演，具體步驟：

（1）根據(jù)觀測(cè)壩體位移變形數(shù)據(jù)統(tǒng)計(jì)回歸分析[2]，得到不同水位下實(shí)測(cè)位移統(tǒng)計(jì)模型水壓分量D以及不確定性因素造成的隨機(jī)誤差 εH～N(0,σH2)；

（2）確定壩體、壩基變形模量X的先驗(yàn)分布π(X)；

（3）由于觀測(cè)的變形具有不確定性，由統(tǒng)計(jì)模型分析得到的壩體位變形與實(shí)際觀測(cè)值必然存在差距，因此，抽樣過程中必須考慮不確定性因素造成的隨機(jī)誤差εH～N(0,σH2)；

（4）以后驗(yàn)概率密度函數(shù)f (X∣D)為目標(biāo)函數(shù)，采用MCMC without likelihoods算法產(chǎn)生壩體、壩基變形模量隨機(jī)樣本，其中通過隨機(jī)誤差εH～N(0,σH2)篩選相應(yīng)的有效樣本，2.4節(jié)通過人為調(diào)整隨機(jī)誤差，進(jìn)一步研究了隨機(jī)誤差對(duì)壩體、壩基變形模量分布的影響；

（5）將有效的馬爾可夫鏈樣本作為后驗(yàn)概率密度函數(shù)f (X∣D)的樣本，計(jì)算后驗(yàn)分布f (X∣D)的估計(jì)值（均值，標(biāo)準(zhǔn)差）。

2 工程算例

2.1 工程概況

龍灘水電站前期工程于2006年10月開始下閘蓄水，2008年12月初步達(dá)到正常蓄水位。大壩壩頂高程382 m，建基面高程216.43 m，壩高165.57 m，壩頂寬度14 m。以11號(hào)壩段為典型壩段，根據(jù)實(shí)測(cè)壩體水平位移反演壩體、壩基變形模量。如圖1所示，二維有限元模型范圍為：向上、下游各延伸250 m、壩基建基面以下延伸235 m。有限元節(jié)點(diǎn)布置時(shí)，盡可能考慮將位移測(cè)點(diǎn)安排在單元節(jié)點(diǎn)上，采用四節(jié)點(diǎn)單元，共2 383個(gè)節(jié)點(diǎn)，2 244個(gè)單元。

圖 1 11 號(hào)壩段有限元模型及垂線測(cè)點(diǎn)Fig. 1 FEM mesh and vertical monitoring points of 11# dam section

參數(shù)反演時(shí)，選取正垂線測(cè)點(diǎn) PL11-1（高程 342.00～379.00 m），PL11-2-2（高程 270.00～342.00 m）與PL11-3（高程222.75～270.00 m）2010年5月28日至2011年8月15日的實(shí)測(cè)順河向位移進(jìn)行分析，依據(jù)式（4）分離壩體位移水壓、溫度及時(shí)變等主要分量，各垂線的復(fù)相關(guān)系數(shù)R分別為0.964，0.951和0.977，反演分析中輸入不同水位下實(shí)測(cè)位移統(tǒng)計(jì)模型水壓分量D，即對(duì)模型依次施加水位為335，340，345，355和360 m的上游靜水壓力，結(jié)果如表1所示。

表 1 不同水位下實(shí)測(cè)位移統(tǒng)計(jì)模型水壓分量DTab. 1 Water pressure component D given by statistical model at different water levels

觀測(cè)值只是一個(gè)樣本序列或樣本值的實(shí)現(xiàn)，也就是說觀測(cè)的變形具有不確定性。因此，觀測(cè)值必須考慮各種不確定因素。由圖2(a)所示，擬合值與實(shí)測(cè)值存在誤差，誤差來源于上述各種不確定因素。根據(jù)1.2節(jié)所述引入水壓分量的隨機(jī)誤差 εH～N(0,σH2)，圖 2(b)，(c)和 (d)分別為 222.75，270.00和 342.00 m 高程擬合總體位移及其分量。

圖 2 龍灘觀測(cè)位移擬合值與實(shí)測(cè)值Fig. 2 Fitted values and measured values of observation displacements of Longtan Dam

2.2 單參數(shù)反演分析

假定壩基變形模量為確定量，壩體彈性模量為待反演隨機(jī)變量。采用無似然函數(shù)的MCMC算法來反演單個(gè)參數(shù)，初步認(rèn)為壩體彈性模量服從區(qū)間為[10，60]的均勻分布，模擬時(shí)馬氏鏈樣本數(shù)取20 000個(gè)。

圖3為無似然函數(shù)MCMC抽樣樣本。圖4為正態(tài)分布Q-Q檢驗(yàn)圖，其中橫坐標(biāo)為分位數(shù)，縱坐標(biāo)為樣本值；圖中點(diǎn)近似為一條直線，說明樣本基本服從正態(tài)分布。圖5為BGR診斷法馬氏鏈?zhǔn)諗啃栽\斷圖，通過比較可以看出樣本數(shù)量達(dá)到7 500時(shí)，診斷指標(biāo)穩(wěn)定在1附近，所以在以下分析中有效鏈長(zhǎng)度取為10 000～20 000。圖6比較了參數(shù)的先驗(yàn)分布、后驗(yàn)分布概率密度曲線及其正態(tài)擬合曲線，結(jié)果表明：（1）壩體彈性模量均值由先驗(yàn)分布的45 GPa變?yōu)楹篁?yàn)分布均值32.09 GPa，基本符合龍灘重力壩實(shí)際壩體彈性模量；（2）彈性模量后驗(yàn)分布區(qū)間比先驗(yàn)分布明顯收窄，壩體彈性模量反演結(jié)果的標(biāo)準(zhǔn)差為3.81 GPa，變異系數(shù)為0.119。

圖 3 MCMC抽樣樣本Fig. 3 Samples in MCMC

圖 4 正態(tài)分布Q-Q檢驗(yàn)圖Fig. 4 Q-Q diagram of normal distribution

圖 5 馬氏鏈?zhǔn)諗啃栽\斷Fig. 5 Convergence diagnosis of Markov chains

圖 6 彈性模量后驗(yàn)分布概率密度分布曲線Fig. 6 Probability density curves of posterior distribution

2.3 多參數(shù)反演分析

假定壩體彈性模量及壩基變形模量均為待反演隨機(jī)變量，采用無似然函數(shù)的MCMC算法反演多個(gè)參數(shù)。反演時(shí)，初步認(rèn)為兩參數(shù)均服從區(qū)間為[10，60]的均勻分布，模擬時(shí)馬氏鏈樣本數(shù)取20 000個(gè)。

圖7為壩體彈性模量、壩基變形模量MCMC抽樣樣本，其Q-Q檢驗(yàn)圖表明兩個(gè)變量的馬氏鏈樣本均服從正態(tài)分布。圖8比較了兩個(gè)參數(shù)的先驗(yàn)分布、后驗(yàn)分布概率密度曲線及其正態(tài)擬合曲線，結(jié)果表明：（1）壩體彈性模量均值由先驗(yàn)分布的45 GPa變?yōu)楹篁?yàn)分布均值33.49 GPa，壩基變形模量的均值由先驗(yàn)分布35 GPa變?yōu)楹篁?yàn)分布均值29.53 GPa，基本符合龍灘重力壩實(shí)際情況；（2）壩體彈性模量、壩基變形模量標(biāo)準(zhǔn)差分別為3.78 和1.90 GPa，變異系數(shù)分別為0.11和0.06，標(biāo)準(zhǔn)差與變異系數(shù)均明顯縮小，說明參數(shù)的不確定性顯著減小。

圖 7 MCMC抽樣樣本Fig. 7 Samples in MCMC

圖 8 參數(shù)后驗(yàn)分布概率密度分布曲線Fig. 8 Probability density curves of posterior distribution

采用經(jīng)典的最小二乘法反演壩體、壩基變形模量，與本文提出的反演方法進(jìn)行對(duì)比，在同一壩段，同樣觀測(cè)數(shù)據(jù)情況下，兩種反演方法的結(jié)果比較如表2所示。結(jié)果表明，兩種算法都能較好地給出參數(shù)的反演結(jié)果，最小二乘法盡管誤差平方和更小，但缺少對(duì)不確定因素的考慮，只給出了一個(gè)確定值的結(jié)果；而MCMC算法充分考慮了不確定因素的影響，反演結(jié)果能更好地反映壩體、壩基變形模量的分布規(guī)律。

表 2 兩種反演方法結(jié)果Tab. 2 Comparison between two inversion methods

2.4 壩體壩基變形模量敏感性分析

位移的不確定性決定了壩體彈性模量及壩基變形模量反演結(jié)果的不確定性。通過模擬不同的σH進(jìn)行壩體彈性模量及壩基變形模量的敏感性分析。正垂線PL11-1（高程342.00 m）與PL11-3（高程222.75 m）測(cè)點(diǎn)實(shí)測(cè)樣本標(biāo)準(zhǔn)差分別為σH= 0.12 mm和σH= 0.24 mm。圖9為不同的觀測(cè)樣本標(biāo)準(zhǔn)差（隨機(jī)誤差）與壩體彈性模量、壩基變形模量標(biāo)準(zhǔn)差的關(guān)系，易發(fā)現(xiàn)隨著觀測(cè)樣本隨機(jī)性（波動(dòng)）的增大，反演的壩體、壩基變形模量變異系數(shù)相應(yīng)增大。

圖 9 觀測(cè)樣本標(biāo)準(zhǔn)差與壩體彈性模量、壩基變形模量變異系數(shù)關(guān)系Fig. 9 Relationships between standard deviations of observation samples and variation coefficients of elastic modulus of dam body and deformation modulus of dam foundation

3 結(jié) 語

針對(duì)混凝土壩材料力學(xué)參數(shù)反演中存在大量不確定性的問題，提出了高混凝土重力壩壩體彈性模量與壩基變形模量的MCMC隨機(jī)反演方法，給出了無似然函數(shù)的馬爾可夫鏈蒙特卡羅法在工程中的實(shí)際應(yīng)用，并研究了壩體變形不確定性對(duì)反演結(jié)果的影響。研究結(jié)果表明：（1）該隨機(jī)反演方法可以得出所需反演參數(shù)（壩體彈性模量、壩基變形模量）的分布；（2）反演參數(shù)的后驗(yàn)分布變異性較先驗(yàn)分布有明顯降低；（3）在考慮變形觀測(cè)值離散性變化情況下，后驗(yàn)分布變異性與觀測(cè)值離散性呈正相關(guān)關(guān)系。

MCMC隨機(jī)反演方法計(jì)算簡(jiǎn)便、高效，具有很強(qiáng)的工程實(shí)用性，為混凝土壩參數(shù)隨機(jī)反演提供了重要計(jì)算方法。但有時(shí)為了得到較為精確的參數(shù)尾部分布，需要的抽樣數(shù)較大，這就導(dǎo)致有限元計(jì)算的總耗時(shí)巨大，因此需要結(jié)合有限元原理對(duì)該算法進(jìn)行優(yōu)化。