貝葉斯理論在既有砌體損傷定位中的應用研究

夏恒，崔鴻知

（長沙理工大學 土木工程學院，湖南長沙 410114）

0 引言

目前，部分既有砌體[1]屬于歷史文化保護遺產，其歷史人文價值與日俱增，但材料老化、荷載作用、結構疲勞、超期服役等因素會導致砌體結構出現退化和損壞的現象，由于結構復雜、損傷位置隱蔽等情況，很難對結構進行準確的損傷定位，假如對其損傷放任不管，隨著損傷的累積，將可能造成重大安全事故和經濟損失，故對這些砌體結構進行損傷識別研究具有重要的意義。

在20 世紀50 年代，我國的傳統檢測方法（筒壓法、扁頂法、原位單磚雙剪法、原位軸壓法、點荷法）屬于定期的人工檢測，在進行檢測時還需預知結構損傷的大致位置，對于一些難以測量的部位，其人工檢測的危險程度也會增加。到了20 世紀80 年代，隨著科學技術的發展，基于動力學的損傷識別[2]方法逐漸成為主流，其原理是結構發生損傷時，相關參數會發生改變，通過在結構上布置相應的傳感器可對結構的健康狀況進行實時監測。

在基于動力學的損傷定位方法上， Pandey[3]提出絕對變化量指標，根據損傷前后應變模態差的絕對值進行損傷定位。董聰[4]根據結構應變的一階變分關系，推導出了應變模態相較于位移振型和自振頻率更能有效、準確地進行損傷定位分析，并且在進行試驗后數據表明：損傷位置應變模態的變化十分明顯，且變化值的大小與損傷程度成正比關系。Vanik[5]提出基于參數與貝葉斯理論相結合的在線健康監測方法。顧培因[6]提出一種在無法獲取損傷前的模態數據下的一種直接指標法，該方法可以直接利用損傷后的應變模態參數進行損傷定位。殷志翔[7]利用應變模態差對索網結構實現了損傷識別。鄒云峰[8]提出了基于應變模態響應重構的損傷識別方法。吳韜[8-10]等人將動力特性用于砌體結構的損傷識別。

研究發現，目前砌體結構的損傷定位研究大部分屬于確定性分析，即損傷前的參數都有準確的數據。但在現實中，大部分的砌體結構因為各種原因，很難獲取損傷前的參數。本文以湖南省長沙市湖南大學財經學院的省級保護性建筑砌體柱為例，考慮測量誤差和模型誤差的影響，以改良的應變模態差值公式嵌入到貝葉斯模型修正的框架中。該方法可以有效地處理損傷定位中的不確定性因素，能準確地識別結構損傷位置，并且還能將此次獲得的后驗信息作為下次更新的先驗信息，從而實現長期的動態觀測。

1 貝葉斯基本理論及MCMC 方法

1.1 貝葉斯基本理論

貝葉斯統計推斷[11]方法是引入概率學的知識來做出判斷，它是一種承認生產實踐的繼續性方法，將相關先驗信息（如歷史可靠數據）參與決策。推斷基礎來自于貝葉斯公式，其公式如下：

式（1）中P（θi）是事件θi 發生的先驗分布，與事件D 無關，先驗信息主要來源于可靠的歷史數據和專家經驗，假設先驗分布為正態分布N（μ，λ2），可以利用先驗信息確定μ 和λ2的值；P(D|θi)是指在事件θi 發生的基礎上，事件D 的條件分布；是一個全概率數，起一個正則化因子的作用；P(θi|D)是在事件D 發生的基礎上，事件θi 的分布，也被稱為后驗分布。一般情況下后驗概率是難以表達的，對后驗概率的計算可以采用MCMC-馬爾科夫鏈蒙特卡羅法。

1.2 MCMC-馬爾科夫鏈蒙特卡羅法

在馬爾科夫鏈蒙特卡羅法中進行抽樣是至關重要的，可以降低其樣本的相關性。不同的抽樣方法，將產生不同的馬爾科夫鏈蒙特卡羅算法。1987 年，Durane[13]將馬爾可夫鏈蒙特卡羅法與分子動力學相結合，提出了混合蒙特卡羅算法。混合蒙特卡羅算法結合了梅特羅波利斯-黑斯挺算法、吉布斯算法以及確定性分子動力學方法，可以用于高維復雜分布的抽樣模擬，改進了建議概率分布的產生機制，避免陷入局部解。

混合蒙特卡羅法的算法流程[14]如下：

（1）鏈指標k=0，初始化位置向量x0和動力向量p0。

（2）鏈指標k=k+1，從正態分布N（0，M）中隨機產生一個新的動力向量pt。

（3）賦予蛙跳算法初值[x(0),p(0)=(xk-1,pt)]并允許L 時間步以獲得一個建議狀態（x*，p*）=[x(t+Lε)，pt+Lε)]。

（4）計算接受概率α=min[1,exp(-△H)],其中△H=H（ x*，p*）-H（xk-1,pt）。如果接受建議狀態（xk，pk）=（x*，p*），否則（xk，pk）= （xk-1，pt）。

（5）重復步驟（2）—（4），直至到達預設的鏈長為止。

再根據后驗樣本可計算后驗分布的各階矩，進行相應的統計推斷。

2 有限元模型修正方法

部分既有砌體由于年代久遠，無法準確獲取結構損傷前的參數，只能根據現場測量數據和相關歷史資料初步建立無損的有限元模型，但所建立的模型與實際情況存在差異。為了修正有限元模型，本文把結構自身物理參數，如材料密度、彈性模量、以及泊松比等作為修正參數，以實際結構的測量值與有限元模型的動力響應特征值的相關信息來構建目標函數[15]，以相關信息差值作為修正依據。

式（2）中：y（r）表示為現實結構與有限元模型動力響應特征值的差值，可以是應力、應變、頻率、振型、阻尼、應變模態等；ME表示為現實結構測量的數據 ；MA（r）表示為有限元模型的特征值，其中r 表示為修正的參數向量。有限元模型修正流程如圖1 所示。

圖1 有限元模型修正流程圖

3 基于貝葉斯理論的損傷定位方法

3.1 基于應變模態的損傷識別原理

在基于動力學的結構損傷定位中，用于損傷定位的參數需要滿足兩個條件：一是對局部損傷敏感；二是位置坐標的單調函數[4]。

在文獻[4]中證實了應變模態相較于其他參數更能有效、準確解決損傷定位問題。

由動力學方程：

文獻[4]推導出結構應變的一階變分函數為：

式（4）中{Δε}為應變的變化；[Δψ]為應變模態的變化；[Уr]結構自振頻率矩陣；[φ]為位移模態矩陣；[ψ]為應變模態矩陣；[ΔУr]結構自振頻率變化；[Δφ]位移模態的變化；其中{Δε}和[Δψ]的變化在位置上存在明確的對應關系，是具有進行損傷識別的能力。

損傷定位的公式以應變模態差的絕對值為依據，其公式為：

式子（5）中[ψI]和[ψD]表示為無損結構和有損結構的應變模態矩陣；i 表示坐標位置；j 表示模態階數。在損傷位置處，應變模態差Δψ會發生突變，Δψ的大小表示損傷的嚴重程度，Δψ 恒大于0。

本文綜合考慮不確定性因素帶來的影響，實際的應變模態差的絕對值調整為：

式（6）中,f1、f2分別是測量誤差和模型誤差，假設f1、f2都服從正態分布，則其均值分別為μ1、μ2，方差分別為σ21、σ22。

3.2 似然函數與后驗分布

由改良的定位公式（6），可推出似然函數表示為：

確定似然函數和先驗分布后，后驗分布[16]則為：

式（9）中，P(D)為邊緣似然函數，其值為固定常數。其中參數θ 中的每一個θj，其邊緣概率密可以通過積分獲得為：

3.3 損傷識別流程

本文將砌體結構識別流程分為三個主要部分：（1）先基于歷史資料初步建立有限元模型，再激勵砌體結構獲取相關參數數據，進行動力特性分析，最后以相關信息差值對有限元模型進行修正；（2）將現場收集的數據作為先驗信息，再通過馬爾科夫鏈產生候選樣本數據，依據混合蒙特卡羅法計算接受概率，確定樣本集；（3）通過貝葉斯方法，計算應變模態的后驗分布，判斷是否收斂。具體流程圖如圖2 所示。

圖2 損傷識別流程圖

4 砌體結構損傷定位研究

4.1 有限元模型材料參數和邊界條件

燒結磚砌體柱來源于湖南省長沙市湖南大學財經學院的省級保護性建筑的第一層柱體，結構尺寸為720 mm×470 mm×25 440 mm，如圖3 所示。實體結構多處帶有明顯裂縫，且砌塊間的砂漿風化程度嚴重，幾乎散失了強度。損傷前的有限元模型依據現場測量數據、相關資料和GB 50003—2011《砌體結構設計規范》確定燒結磚為MU15、砂漿為M7.5，砌塊密度設為1 800 kg/m3、彈性模量設為1 600（ff 為砌體抗壓強度，單位為MPa），計算為3 312 MPa、體積剪切模量為1 577 MPa、剪切模量1 440 MPa、線性膨脹系數設為5×10-6/℃、收縮率設為-0.1 mm/m、摩擦系數設為0.70、熱導率設為1.4 W/(m·K)、泊松比設為0.15。

圖3 燒結磚砌體柱

有限元模型采用ANSYS Workbench 建模，單塊磚的模型為240 mm×115 mm×53 mm，邊界條件設置為：上部承受荷載為均布荷載，方向垂直向下無偏心，根據GB 50009—2012《建筑結構荷載規范》和設計資料，其均布荷載值設為20 kN/m2，底部采用 Fix supporied（固定約束），加速度設為9.8 m/s2，方向豎直向下。由于模型整齊，單元類型采用六面體SOLID186 單元，單元尺寸為50 mm，網格劃分采用程序自動控制。有限元模型如圖4 所示。

圖4 有限元模型及檢測點布置圖

4.2 試驗方案和數據結果

在實體結構上取五點分別為a、b、c、d、e，這五點處于結構正面，如圖4 所示。其中a、b、c 三處是帶有明顯裂縫的位置，d、e 兩處無裂縫，通過比較現實結構和有限元模型的應變模態差值判斷結構損傷位置，差值越大說明損傷越嚴重。對實體結構a、b、c、d、e 測量點分別測量3 次，以其平均值作為先驗樣本（表1），再經過混合蒙特卡羅法計算，設置樣本接受的概率，經過多次循環后，使樣本值向目標值收斂。由于前面的樣本值高度依賴于初始值，故舍棄前200 個樣本。本文分析前四階模態下應變模態差值的后驗分布。

表1 實際結構檢測點應變模態

有限元模型在修正參數后得到無損傷前應變模態云圖（圖5），檢測點a、b、c、d、e 的應變模態如表2 所示，實際結構和有限元模型檢測點的應變模態如圖6 所示。

表2 有限元模型檢測點應變模態

圖5 前四階應變模態云圖

圖6 實際結構與有限元模型檢測點應變模態

從圖6 中可以明顯發現：測量點處應變模態值的大小與結構損傷程度成正相關。測量點（a、b、c）在一階模態和四階模態下的應變模態值明顯高于二階模態和三階模態下的應變模態值，其應變模態差值在一階模態和四階模態下有著更為明顯的峰值，故在損傷定位時，針對此結構采用一階模態和四階模態能更準確地識別損傷位置。在檢測點（d、e）處，各階模態上的應變模態差值相較于測量點（a、b、c），應變模態差值保持在穩定低值范圍內。將五點的應變模態差值作為先驗信息帶入計算模型中，一階模態下測量點處的應變模態差值統計結果如表3 所示。

表3 測量點處一階應變模態差值

考慮測量誤差和模型誤差的影響（式6），其均值和標準差在不同的采樣點取值不同，將修正的應變模態差值的均值帶入馬爾科夫鏈中，生成樣本集，通過貝葉斯算法確定后驗分布圖（圖7），具體數值如表4 所示。

表4 測量點處一階應變模態差值的修正值

圖7 各測量點應變模態差值的后驗分布圖

從圖7 和表4 中可以發現經過貝葉斯推斷后，其應變模態差值的均值和標準差有相對明顯的修正，在帶有裂縫的a、b、c點應變模態差值與沒帶裂縫d、e 點應變模態差值存在明顯不同，在帶有裂縫處的檢測點其應變模態差值遠遠高于沒有裂縫處的檢測點，如圖8 所示。通過應變模態差值的大小可以判斷出a、b、c 點相較于d、e 點是存在損傷，并且其損傷的程度與應變模態差值成正相關。

圖8 應變模態差值折線圖

5 結論

貝葉斯理論屬于一種概率統計理論，該方法可借助先驗信息（包括可靠的歷史數據以及專家經驗），再基于總體信息和樣本信息得出合理推斷。本文將該方法與有限元方法共同運用于砌體柱的損傷定位研究，結論如下：

（1）砌體柱的損傷前有限元模型可以通過現場測量數據和相關資料以及經驗進行初步的建立，利用實際結構與有限元模型的動力特征值的差值作為修正依據，以貝葉斯算法作為優化方法，做到信息差最小，完成有限元模型的建立。

（2）貝葉斯推斷和有限元相結合的方法在結構損傷定位研究中具有可行性，貝葉斯理論可以有效地處理損傷定位中不確定性因素的影響，從而得出更為合理的推斷。