基于貝葉斯的鋼筋混凝土梁模型修正方法

宋彥朋，陳 輝,，黃 斌

(1.武漢工程大學 郵電與信息工程學院，湖北 武漢 430073；2.武漢理工大學 土木工程與建筑學院，湖北 武漢 430070)

在土木工程結(jié)構(gòu)中，計算模型通常采用有限元格式。土木工程結(jié)構(gòu)尤其是城市道路橋梁結(jié)構(gòu)，主要是由鋼筋混凝土構(gòu)成。混凝土或鋼筋混凝土結(jié)構(gòu)相對于鋼結(jié)構(gòu)而言，材料內(nèi)部的非均勻性、各向異性、樣本參數(shù)的個體間差異顯著。因此有限元仿真模型與實際結(jié)構(gòu)存在比鋼結(jié)構(gòu)更為顯著的差異。為了減小差異并建立足夠精確的有限元模型，有必要根據(jù)測量數(shù)據(jù)更新有限元模型的參數(shù)[1～3]。

使用測量數(shù)據(jù)來修改結(jié)構(gòu)的計算模型是一個非?；钴S的研究方向。傳統(tǒng)的基于動力的模型修正方法主要是基于結(jié)構(gòu)固有頻率、振型、模態(tài)曲率、柔度矩陣等[4]。這些方法基本都是基于特征值方程的動力學反問題。在動力試驗中往往只能得到前幾階模態(tài)數(shù)據(jù)，因此在基于有限元方法的模型修正中往往會求解不適定方程，造成病態(tài)解[5]。很多學者采用不同的方法來避免病態(tài)解。例如，Ren等[6]采用模態(tài)截斷奇異值分解算法來修正自由模態(tài)的鋼筋混凝土實驗梁，并將修正后的模型作為基準模型用于該梁的損傷識別，取得了很好的效果；Hua等[7]將Tikhonov正則化用于平面數(shù)值桁架的有限元模型修正來避免修正方程求解的病態(tài)問題；李英超等[8]將奇異值截斷正則化算法(Truncated Singular Value Decomposition，TSVD)用于三維導管架平臺的縮尺簡化模型的模型修正，很好地避免了多解問題。從以上方法中可以看出，正則化技術是處理方程病態(tài)問題的一種非常有效的手段，但正則化參數(shù)或者奇異值截斷數(shù)以及測量誤差對求解結(jié)果也有較大影響，另外，這些文章多半都是基于數(shù)值算例或者單一材料的實驗算例，對鋼筋混凝土結(jié)構(gòu)的研究較少。另一方面，基于Beck等[9,10]提出的基于貝葉斯框架的模型修正或結(jié)構(gòu)參數(shù)識別的方法也取得了巨大的發(fā)展，很多學者也將其應用到混凝土結(jié)構(gòu)的有限元模型修正中。該方法最大的特點是通過隨機抽取參數(shù)樣本，正向計算結(jié)構(gòu)響應并和測量結(jié)果比較來選擇合理的待修正參數(shù)，從根本上避免逆向求解問題。易偉建等[11]運用貝葉斯方法，通過測量頻率和振型MAC值構(gòu)建目標函數(shù)，對一個四層混凝土框架結(jié)構(gòu)進行了模型修正，并將修正后的模型作為基準模型用于該結(jié)構(gòu)的損傷識別。房長宇等[12]用貝葉斯方法基于測量頻率對預應力混凝土梁進行了模型修正。Wan等[13]基于測量頻率，將貝葉斯方法用于實際混凝土人行天橋的模型修正中。從上述基于貝葉斯框架的模型修正方法中可以看出，該方法最重要的兩個問題是目標函數(shù)的構(gòu)建和抽樣方法的選取。

考慮到混凝土材料的特殊性以及施工及環(huán)境的影響，混凝土結(jié)構(gòu)各區(qū)域材料參數(shù)差異較其他材料顯著，因此需要盡可能地對其作全局修正，但這會導致修正參數(shù)較多，傳統(tǒng)的吉布斯(Gibbs)抽樣方法或MH(Metropolis-Hasyings)抽樣方法限制了該方法的使用；另外，要進行全局修正，通常需要對大小發(fā)生變化的修正參數(shù)進行定位，但和振型相關MAC值只能體現(xiàn)整體模態(tài)差異，無法體現(xiàn)局部差異，而振型是位置坐標的函數(shù)，因此在目標函數(shù)的構(gòu)建上應該充分使用該測量信息?；谪惾~斯方法的最新進展，本文考慮采用貝葉斯模型修正方法對鋼筋混凝土實驗梁進行模型修正，以避免模型修正中的病態(tài)問題。同時基于測量頻率和模態(tài)來構(gòu)建目標函數(shù)，采用更為高效的DRAM(Delayed Rejection Adaptive Metropolis)抽樣方法來提高抽樣的遍歷性。為了確認結(jié)構(gòu)真實模型，文章通過對比實際物理特征對修正結(jié)果進行了確認。

1 貝葉斯模型修正的基本原理

基于貝葉斯方法的隨機模型修正是結(jié)合先驗知識(主觀信息)和測試數(shù)據(jù)(客觀信息)，采用MCMC(Markov Chain Monte Carlo)方法推斷修正參數(shù)的后驗概率分布。其表達式如下：

=cp(x/θ)π(θ)∝p(x/θ)π(θ)

(1)

式中：x為觀測信息；π(θ)為待修參數(shù)θ的先驗分布；p(x/θ)為在θ給定下的條件分布，通常稱為似然函數(shù)；c是一個不依賴于θ的常數(shù)因子。

在土木工程或者機械結(jié)構(gòu)中，一般選用結(jié)構(gòu)的彈性模量、剛度、密度、質(zhì)量或某些幾何尺寸作為修正參數(shù)，標記為θ，觀測值記為DN，N為觀測次數(shù)。將上述貝葉斯原理引入到有限元模型修正中，可將修正參數(shù)的后驗概率分布寫為：

p(θ/DN)=cp(DN/θ)p(θ)

(2)

式中：p(θ/DN)為結(jié)構(gòu)后驗概率；p(DN/θ)為似然函數(shù)；p(θ)為待修正參數(shù)θ的先驗分布；μ0為先驗參數(shù)均值，covσ0為先驗參數(shù)的協(xié)方差矩陣；y為測量響應；y(θ)為經(jīng)模型仿真計算(如有限元計算)得到的響應；covy為測量信息協(xié)方差矩陣；c為與θ無關的常數(shù)。根據(jù)貝葉斯假設采用無偏廣義先驗，p(θ)值應為1，待修正參數(shù)的后驗概率可以寫為：

(3)

式中：J(θ)為目標函數(shù)，可表示為：

(4)

2 目標函數(shù)的建立

結(jié)構(gòu)的響應通常為頻率和振型，在基于貝葉斯方法的模型修正方法中，常用的是頻率、振型和MAC值或者它們的組合來構(gòu)造目標函數(shù)。實際模態(tài)測試中，頻率的測試較為容易和精確，而且不會受測點缺失的影響；但模態(tài)振型的測試中僅前幾階低階模態(tài)能較為準確測量到。同時在測試過程中也經(jīng)常會出現(xiàn)部分傳感器信號缺失或者故障的情況發(fā)生，這會造成模態(tài)信息不完整。因此，本文從實際模態(tài)測試角度構(gòu)造目標函數(shù)以修正結(jié)構(gòu)參數(shù)。

考慮每次測量的頻率與振型和抽樣計算結(jié)果的誤差，可將目標函數(shù)(4)進一步定義為：

(5)

將式(5)表示的目標函數(shù)J(θ)代入式(3)中可以得到待修正參數(shù)的后驗概率密度的表達式。為了找到修正參數(shù)的后驗概率密度函數(shù)的最大值，可用MCMC方法近似計算后驗概率密度函數(shù)的分布，即采用MCMC抽樣使目標函數(shù)J(θ)最小，最終可得到最優(yōu)修正參數(shù)。

3 計算后驗概率密度的DRAM抽樣方法

在采用MCMC抽樣計算待修正參數(shù)后驗概率密度的過程中，產(chǎn)生Markov鏈的常用方法是Gibbs 抽樣[14]和MH抽樣[15]。這兩種方法的不足是當參數(shù)較多時，采樣過程平滑段出現(xiàn)逐漸增多，且參數(shù)之間易產(chǎn)生相互影響最終難以收斂。近年來，發(fā)展了一種新的抽樣方法即DRAM抽樣，它是一種更高效的抽樣方法[16]，本文采用了這種抽樣方法。該方法的抽樣步驟如下：

(1)假設在t時刻，Markov鏈當前值為θt，然后根據(jù)建議分布q(θt,C0)隨機抽取候選樣本θt+1，其中C0為初始抽樣方差。

(2)根據(jù)式(3)，計算接受概率p。

(3)從[0，1]均勻分布中隨機產(chǎn)生一個變量u，當p>u時，接受候選樣本θt+1，即θt+1=θt；當p≤u時，拒絕候選樣本θt+1。

(4)根據(jù)二次建議分布q2(θt,C0)隨機抽取候選樣本θt+1，根據(jù)式(3)，計算接受概率p2，然后執(zhí)行步驟(3)。

(5)當執(zhí)行到第t+1步，且t+1≥N0(非適應抽樣次數(shù))時，更改建議分布的協(xié)方差矩陣為Ct+1，從(1)開始執(zhí)行。其中任意抽樣時刻t建議分布的協(xié)方差矩陣Ct滿足：

(6)重復步驟(1)～(5)，直至Markov鏈趨于平穩(wěn)后迭代終止。然后計算從收斂開始后修正參數(shù)樣本的均值和方差等統(tǒng)計特征。

4 基于貝葉斯的鋼筋混凝土梁有限元模型修正

4.1 修正參數(shù)的選擇

(6)

式中：Δ[K]i為結(jié)構(gòu)第i個單元的剛度變化量；αi為結(jié)構(gòu)第i個單元的修正系數(shù)。為確保本文提出方法的有效性，在數(shù)值算例中驗證了幾種實際測量中常見工況下的模型修正效果。

4.2 基于貝葉斯的鋼筋混凝土梁有限元模型修正

對鋼筋混凝土梁的初始有限元模型進行計算，得到其自振頻率和位移振型。同時對實際模型進行模態(tài)測試，分析并得到有效的測量信息。將初始模型自由度和測量自由度進行匹配，得到匹配后的初始有限元模型對應的頻率和位移振型。將匹配后的所有振型進行質(zhì)量歸一化。最后將測量頻率和計算頻率以及歸一化以后的位移模態(tài)代入式(3)中，用DRAM方法進行抽樣并計算得到實際結(jié)構(gòu)修正參數(shù)的后驗概率，最終可以得到修正參數(shù)的優(yōu)值。

4.3 修正結(jié)果驗證

模型修正的目的是使得修正后結(jié)構(gòu)的計算響應與測量結(jié)果一致，同時修正后的參數(shù)必須能和真實結(jié)構(gòu)的物理特征對應。因此，在用本文提出的方法得到鋼筋混凝土梁的修正參數(shù)后，首先應將修正參數(shù)代入修正后的有限元模型中計算響應，并和測量結(jié)果進行對比以評估修正結(jié)果的有效性；同時，為了確保修正結(jié)果具有真實意義，還應該對鋼筋混凝土梁構(gòu)造明顯的物理特征(如可見的細小裂縫)，來進一步驗證修正結(jié)果是否具有真實意義。

5 數(shù)值算例

取一懸臂梁，如圖1所示，其密度ρ=7.8×103 kg/m3，初始彈性模量E=2.10×l011Pa。將梁劃分為七個單元，假定懸臂梁彈性模量發(fā)生了變化，變化后各單元彈性模量真實值如表1所示。假定待修正參數(shù)為各單元彈性模量，相對于初始彈性模量E的變化量用參數(shù)αi(i=1,2,…,7)表示。實際模態(tài)測試中，測量信息往往有限，同時高階模態(tài)測量誤差較大。能有效使用的往往是少數(shù)低階模態(tài)。因此，需要就測量模態(tài)階數(shù)對參數(shù)識別結(jié)果的影響進行分析。為此，分四種工況討論，每種工況的測量信息為：(1)結(jié)構(gòu)前二階頻率和前二階位移振型；(2)結(jié)構(gòu)前三階頻率和前三階位移振型；(3)結(jié)構(gòu)前四階頻率和前四階位移振型；(4)結(jié)構(gòu)前四階頻率。

圖1 懸臂梁有限元模型/mm

表1 各單元彈性模量真值

分別將以上四種工況下的測量頻率和位移模態(tài)信息代入目標函數(shù)進行參數(shù)識別。仿真測量誤差假定為變異系數(shù)為0.01的正態(tài)分布隨機變量。四種工況下各識別參數(shù)樣本的Markov鏈如圖2所示。在工況1～3下，懸臂梁中待識別參數(shù)的均值如圖3所示。

圖2 四種工況下采用DRAM抽樣的各單元彈模變化參數(shù)樣本的Markov鏈

圖3 懸臂梁中待識別參數(shù)的均值

從圖2a～2c中可以看出，不同測量信息下，三種基于頻率與振型識別結(jié)果中，所有單元的修正參數(shù)均在真值附近波動，并呈現(xiàn)穩(wěn)定的收斂現(xiàn)象，說明本方法識別結(jié)果具有很強的穩(wěn)定性；另外，還發(fā)現(xiàn)可靠的測量信息越豐富，Markov鏈的穩(wěn)定性越高。另外，圖2d的結(jié)果說明，僅基于前四階測量頻率的Markov鏈并不收斂，說明僅依賴于前幾階低階頻率的貝葉斯方法無法有效識別較多結(jié)構(gòu)參數(shù)的變化，可能需要更多高階頻率參與修正，或采用多層MCMC方法。綜合整個仿真算例來看，本文提出的方法在測量的少數(shù)低階頻率與振型信息下，能有效修正結(jié)構(gòu)模型。

6 混凝土梁的模型修正試驗

通過對一實際鋼筋混凝土梁進行模型修正，以驗證本文提出方法的有效性。試驗梁幾何尺寸和配筋等情況如圖4a所示。實際模型如圖4b所示。混凝土的力學參數(shù)為：密度ρ= 2400 kg/m3，彈性模量E=26 GPa，剪切彈性模量G= 0.4E，泊松比μ=0.33。充分考慮梁的剪切效應，選用Timoshenko梁單元類型梁單元，梁被分成110個單元用于有限元計算以保證計算結(jié)果收斂。修正過程中，考慮到修正參數(shù)過多容易使抽樣結(jié)果發(fā)散，因此將修正模型作如下處理：將110個單元中的1～10單元作為一個修正參數(shù)(單元)，以此類推。因此總共修正參數(shù)(單元)為11個。

圖4 鋼筋混凝土試驗梁尺寸和動力測試

為了準確地驗證修正結(jié)果的物理意義，將試驗梁緩慢加載至輕微裂縫出現(xiàn)后停止加載，使之產(chǎn)生明顯的裂縫這一物理特征，但同時要保證試驗梁處于正常使用極限范圍之內(nèi)。在裂縫處做好紅色標記以便于和修正結(jié)果比較。最終在單元4，6，7，8處有不同程度的裂縫產(chǎn)生，如圖4b所示，為非對稱裂縫。其中4，6單元的裂縫深度在梁深度的一半附近，7單元裂縫達到裂縫深度40%附近，基本在跨中，8單元裂縫深度接近梁深度的20%，各裂縫寬度均不超過0.5 mm；將試驗梁放置在軟橡膠墊塊上進行自由模態(tài)測試。模態(tài)試驗中，采用六個加速度傳感器分三批對12個測點進行測量，傳感器布置如圖4b所示，加速度響應圖和頻譜圖如圖4c，4d所示。采用工作模態(tài)分析(Operational Modal Analysis，OMA)識別結(jié)構(gòu)頻率和模態(tài)，所得結(jié)果如表2和圖5所示。各階頻率的測量誤差變異系數(shù)為0.02，在1 h內(nèi)總共進行了12次模態(tài)測試。將12次測量模態(tài)各測點變異系數(shù)中的最大值作為振型的變異系數(shù)，為0.02。測量誤差均假定為零均值的正態(tài)分布。其中振型已進行質(zhì)量歸一化。為了提高修正精度和計算效率，認為梁兩端的單元剛度基本一致，跨中出現(xiàn)裂縫的單元剛度變化較大。因此，考慮到試驗梁質(zhì)量保持不變，將1～3單元和9～11單元的剛度相對于初始有限元剛度的變化量分別作為兩個待修正參數(shù)，中間的4～8單元的剛度相對于初始有限元剛度的變化量作為另外五個待修正參數(shù)，總共7個待修正參數(shù)。取測量結(jié)果中的前四階頻率和位移振型代入目標函數(shù)對待修正參數(shù)進行抽樣。值得注意的是在式(4)中調(diào)整目標函數(shù)中頻率與位移振型的權重比例為10000以使兩者數(shù)量級一致。

表2 測量頻率 Hz

圖5 混凝土梁初始有限元模型和實際測量的振型

圖6為本文方法得到的Markov鏈，從圖中可以看出，抽樣超過3000次后各單元待修正參數(shù)開始收斂，因此取燃燒期為3000。根據(jù)收斂后的Markov鏈進行最大后驗估計，得到修正參數(shù)均值和標準差，如圖7所示。對比圖7中的剛度變化和圖4b中靜力加載后記錄的裂縫位置和深度，可以發(fā)現(xiàn)：剛度下降的位置和裂縫位置對應，裂縫深度占梁高度的比例和剛度下降量比較接近，單元4中裂縫深度為梁高的50%，識別結(jié)果顯示剛度下降約為40%；單元5中裂縫深度為梁高的50%，但裂縫寬度很小，識別結(jié)果顯示剛度下降接近10%，單元6中裂縫深度為梁高的60%，識別結(jié)果顯示剛度下降量接近57%，單元7中裂縫深度為梁高的40%，識別剛度下降量接近38%；單元8中裂縫深度為梁高的15%，識別剛度下降量接近14%；同時注意到梁首端三個單元剛度略有下降，末端三個單元的剛度均有所增大。整個修正結(jié)果和試驗梁的實際剛度變化位置和大小情況基本吻合。

圖6 試驗混凝土梁的DRAM抽樣結(jié)果

圖7 混凝土梁的貝葉斯模型修正結(jié)果

為了進一步對修正結(jié)果進行驗證，將修正后的參數(shù)代入新的有限元模型中進行計算，得到修正后的結(jié)構(gòu)響應。圖8比較了修正前后試驗梁的前六階計算頻率相對于測量頻率的誤差。從圖8中可以看出，修正前結(jié)構(gòu)頻率和測量頻率之間的最大誤差接近26%，最小頻率誤差也接近10%；修正后試驗梁前四階計算頻率和測量頻率之間的相對誤差很小，未參與修正過程的第五、六兩階頻率相對于測量結(jié)果也有大幅下降。這進一步說明修正結(jié)果是有效的。

圖8 修正前后頻率誤差

圖9比較了修正前后的前四階計算位移振型和初始有限元振型。從圖9中可以看出，測量的振型和初始有限元振型有一定的誤差，并且隨著振型階次的提高，測量振型和初始值之間的誤差越來越大。這再次顯示了仿真和試驗的區(qū)別；圖10是修正前后前四階振型MAC值，可以發(fā)現(xiàn)，經(jīng)過修正后的計算振型和測量結(jié)果更加接近，前四階振型MAC值均更接近1。另外值得注意的是，從圖9e可以看出，第五階測量振型明顯偏離實際情況的允許范圍，這時因為高階的振型由于儀器精度等原因的限制，會產(chǎn)生脫離實際的誤差，因此若作為測量信息使用則會產(chǎn)生錯誤的修正結(jié)果。

圖9 初始有限元和修正前后的前五階位移振型

圖10 修正前后模態(tài)MAC值

7 結(jié) 論

本文針對鋼筋混凝土梁模型修正中的方程病態(tài)問題及模態(tài)測試中測點信息不完整和模態(tài)信息不完備的實際情況，采用貝葉斯方法進行模型修正，目標函數(shù)的構(gòu)建充分考慮了模態(tài)測試的特點，抽樣方法采用高效的DRAM抽樣方法。數(shù)值算例結(jié)果表明，基于DRAM抽樣和模態(tài)測量數(shù)據(jù)的貝葉斯模型修正方法能有效修正結(jié)構(gòu)有限元模型參數(shù)，即使在僅用前兩階模態(tài)的情況下也能進行準確修正，說明該方法具有很好的魯棒性和靈活性。模態(tài)試驗結(jié)果顯示：該方法得到的修正結(jié)果能使試驗梁結(jié)構(gòu)修正后的前六階計算頻率與測量值吻合；同時使前四階振型的MAC值也進一步提高并非常接近1；且靜力加載導致的裂縫位置和程度與修正的位置和大小有很好的對應關系，說明了該方法能有效應用于實際鋼筋混凝土梁結(jié)構(gòu)的模型修正。