基于MCMC法的混凝土壩壩體壩基變形模量隨機反演

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(5) 計算接受率h：

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

抽樣過程中，抽樣終止以馬爾可夫鏈收斂為前提，因此，馬爾可夫鏈收斂與否對模型參數(shù)后驗估計具有重要影響。BGR診斷法[13]是馬爾可夫鏈收斂性診斷常用方法之一，該方法基于區(qū)間統(tǒng)計特性及長度進行診斷，定義診斷指標為：

式中：L 為 馬爾可夫鏈總的序列區(qū)間長度或統(tǒng)計特性；l為馬爾可夫鏈單鏈區(qū)間長度均值或統(tǒng)計特性均值。若馬爾可夫鏈收斂，值應接近于1。

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

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

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

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

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

1.3 壩體壩基變形模量隨機反演

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

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

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

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

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

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

2 工程算例

2.1 工程概況

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

圖 1 11 號壩段有限元模型及垂線測點Fig. 1 FEM mesh and vertical monitoring points of 11# dam section

參數(shù)反演時，選取正垂線測點 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日的實測順河向位移進行分析，依據(jù)式（4）分離壩體位移水壓、溫度及時變等主要分量，各垂線的復相關(guān)系數(shù)R分別為0.964，0.951和0.977，反演分析中輸入不同水位下實測位移統(tǒng)計模型水壓分量D，即對模型依次施加水位為335，340，345，355和360 m的上游靜水壓力，結(jié)果如表1所示。

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

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

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

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

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

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

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

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

圖 5 馬氏鏈收斂性診斷Fig. 5 Convergence diagnosis of Markov chains

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

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

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

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

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

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

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

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

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

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

圖 9 觀測樣本標準差與壩體彈性模量、壩基變形模量變異系數(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é) 語

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

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