基于Gibbs抽樣的馬爾科夫蒙特卡羅方法在結構物理參數識別及損傷定位中的研究

劉書奎，吳子燕，張玉兵

（西北工業大學 力學與土木建筑學院，西安 710072）

近些年來結構物理參數的識別引起了人們的廣泛關注［1－4］。但由于結構材料的不均勻性，施工質量的易變性，作用荷載的隨機性等不確定性因素的影響，觀測數據和結構模型均具有強烈的本質不確定性，從而導致結構物理參數識別問題成為不確定性問題［5］。因此，必須在確定性物理參數識別研究的基礎上，發展能夠合理反映不確定性特性的物理參數識別概率方法。

貝葉斯方法能夠綜合先驗和條件信息，并且給出明確的后驗信息，因此被廣泛地應用于不確定性問題。在土木工程領域，Beck等人［6］首先提出了系統辨識的貝葉斯統計模型，但這是一種漸漸逼近的方法，要保證其精確性就必須得到足夠多的模態數據。為解決該問題，Beck和 Au等人［7］又提出了基于 Metropolis-Hastings抽樣的馬爾科夫蒙特卡羅方法，該方法的一個局限在于僅對低維問題有效。Ching和Cheng［8］提出了變異的蒙特卡羅馬爾科夫方法（TMCMC），避免了直接從后驗分布中抽樣，但構造的馬爾科夫鏈的“burn-in”階段非常長，基于Gibbs抽樣的馬爾科夫蒙特卡羅方法可以很好的解決高維參數識別問題［9］。文獻［5］中利用了Gibbs抽樣來進行物理參數概率識別，但所利用的初始條件是確定性的結構模態參數。由于環境影響、測試誤差以及分析過程中的簡化，引入了不確定因素，造成實測模態參數與真值之間不可避免地存在偏差［10］。

本文視結構模態參數為隨機變量，根據結構動力特征方程，將模態參數與結構物理參數通過線性結構識別模型聯系起來，采用基于Gibbs抽樣的馬爾科夫蒙特卡羅方法來對結構物理參數進行后驗抽樣，進而根據后驗樣本的統計特性進行結構物理參數識別及損傷定位。

1 線性回歸模型

線性回歸模型通常表示如下［11］：

其中，Y∈Rn為觀測模型輸出，X∈Rn×m為觀測回歸量矩陣，θ∈Rm為模型參數向量，Xθ為理論模型輸出，e∈Rn為預測誤差，服從均值為零，協方差矩陣為的n維正態分布。

假設模型參數θ的先驗分布為不正常的均勻分布，即為：

預測誤差方差的先驗分布為倒伽馬分布IG（α，β），即為：

當α=β=0時，該先驗分布就轉變為通常的無信息先驗分布，即：

由貝葉斯估計理論，模型參數θ和預測誤差方差的后驗聯合分布為：

若使用 Gibbs抽樣，還需知道θ、的條件后驗分布［11］：

2 線性結構識別模型

在結構健康監測領域，通常用線性結構模型來進行模型更新。一般振動測試數據都是在很低幅值的激勵下得到的，因此很多結構包括損傷結構都近似地表現出線性行為［9］。

結構第i階模態的理論頻率和振型向量φi應滿足如下特征方程：

其中φi=［φi1，φi2，…，φim］T，m為結構自由度數目。當采用集中質量矩陣時，該方程可以表示為：

其中，θ1，θ2，…，θn為無量綱的歸一化參數（以下稱其為結構剛度參數），且恒大于0，表示某一構件對總體剛度矩陣的貢獻，當θi小于1時，即可判斷該構件發生損傷，由θi大小，可判斷損傷的程度。由于結構質量對損傷敏感程度較低，在上述方程中視結構質量為常量，當需要對其進行識別時，運用同樣的處理方法對式（10）做適當變換即可。

對式（10）做變換，并加入預測誤差項，可得線性結構識別模型：

對應于線性回歸模型Y=Xθ+e。

3 Gibbs抽樣過程

本文通過Matlab編程，按照以下步驟實現Gibbs抽樣：

（2）a.根據式（11）并結合的條件后驗分布即式（7），抽取

4 損傷概率

剛度參數向量θ中的元素與結構構件存在一一對應關系，本文利用θ中某一元素出現降低來標識損傷的位置，以其降低的比例d來刻畫損傷的程度。

5 算例分析

算例為某三層剪切型結構（圖1）。本文模擬了三種損傷類型：類型1（DP1）：結構第一層的剛度降低30%；類型2（DP2）：結構第二層的剛度降低30%，第三層的剛度降低50%；類型3（DP3）：結構第一層剛度降低50%、第二層剛度降低40%、第三層剛度降低30%。在進行Gibbs抽樣之前，需要先識別出結構的模態參數（本文中利用的是結構的圓頻率），識別結構模態參數的方法有時間序列法、隨機減量法、NExT、隨機子空間法以及最近由任宜春等提出的基于改進L－P小波的時變模態參數識別方法［12］等。為考慮結構自振圓頻率的隨機性，本文假設其服從以其理論值為均值的正態分布，其分布規律如表1所示。在Gibbs抽樣過程中結構自振圓頻率均根據該統計規律隨機產生。

圖1 三層剪切型結構Fig.1 Three-story shear-building

表1 結構一階自振圓頻率統計特性Tab.1 Statistical properties of the first-order natural frequency

方便起見，對以上四種類型，均取剛度參數向量初值θ=［0.8 0.8 0.8］T，θ1、θ3、θ3的 Gibbs抽樣曲線如下圖所示。

由圖2可以看出，不同損傷類型，當結構某處發生剛度降低時，對應的結構剛度參數θi的抽樣曲線也會發生明顯變化，從而可以清楚地判斷出結構的損傷位置及程度。

圖3顯示的是不同參數空間下的抽樣結果，參數之間表現出近似的線性關系，這表明了識別結果的可靠性。圖3同樣可以判斷損傷的位置及程度，如在圖3.1中，以“+”號表示的圖像（DP1）相比以“o”號表示的圖像（UD）在θ1方向上發生了明顯的左移，而在θ2方向上并未發生明顯的偏移，這說明在DP1下結構的第一層柱子k1發生了剛度降低，而結構的第二層柱子k2并未明顯地出現剛度變化。值得注意的是，在圖3.3中以“+”號表示的圖像（DP1）和“o”號表示的圖像（UD）發生了重合，這是由于結構在損傷類型DP1下僅第一層柱設定了剛度降低，第二層和第三層柱未設定損傷，圖3.3 是在（θ2，θ3）空間的抽樣結果，這又一次驗證了識別結果的可靠性。

由式（12），可以得到結構某一層柱子在不同損傷類型下發生不同剛度降低的超越概率（圖4），亦即其降低的超過某一比例d的概率，當d取值不同時，超越概率也就不同，這樣取一系列d值，最終得到圖4。該圖除了能夠判斷結構的損傷位置及程度外，還可以給出結構在某處損傷程度的概率度量，這也正是貝葉斯方法在處理不確定性問題時的優越之處。

圖4 結構在不同損傷類型下對應不同剛度降低的超越概率Fig.4 Estimated damage probability cures for all damaged patterns

表2 結構剛度參數后驗樣本的統計特性Tab.2 Posterior statistical properties of the stiffness parameters

表2為結構剛度參數后驗樣本的統計特性。由于在結構響應中加入了噪聲，考慮了結構自振圓頻率的隨機性，且僅采用了結構一階頻率等因素的影響，當損傷位置及程度增加時，識別的誤差有所增大，但能滿足工程基本要求。

6 結論

本文首先通過對結構的動力特征方程進行一系列的變化，得到線性結構識別模型，進一步由貝葉斯更新理論得到模型中參數的后驗分布形式。利用結構的模態參數，并考慮其隨機性，應用基于Gibbs抽樣的馬爾科夫蒙特卡羅方法對線性結構識別模型的條件后驗分布進行了抽樣，避免了對后驗分布的直接抽樣，并克服了Metropolis-Hastings抽樣僅對低維問題有效的缺陷，成功地實現了結構物理參數識別及損傷定位。數值算例表明：Gibbs抽樣結果除了可以以不同的方式標識結構的損傷程度及位置且識別的誤差較小外，還能給出識別值的概率度量。

［1］Moaveni B，He X，Conte J P.Damage identification of a composite beam using finite element model updating［J］.Computer-Aided Civil and Infrastructure Engineering，2008（23）：339－359.

［2］張 坤，羅紹湘，段忠東.基于動態響應靈敏度同步反演結構物理參數與輸入的算法［J］.振動與沖擊，2009，28（9）：143 －148，167.

［3］Muto M M.Application of stochastic simulation methods to system identification［D］.California Institute of Technology，Pasadena，California，2007.

［4］Cheung S H，Beck J L.Bayesian Model Updating Using Hybrid Monte Carlo Simulation with Application to Structural Dynamic Models with Many Uncertain Parameters ［J］.Journal ofEngineering Mechanics， 2009， 135（4）：243－255.

［5］李小華.結構參數識別的貝葉斯方法及應用研究［D］.中國地震局工程力學研究所，2009.

［6］ Beck J L，Katafygiotis L S.Updating models and their uncertaintiesI：Bayesian statistical［J］. Journalof Engineering Mechanics，1998，124（4）：455 －461.

［7］Beck J L，Au S K.Bayesian updating of structural models and reliability using Markov chain Monte Carlo simulation［J］. JournalofEngineering Mechanics，2002，128：380－391.

［8］ Ching J，Chen Y J.Transitional markov chain monte carlo method for bayesian model updating，model class selection，and model averaging［J］.Journal of Engineering Mechanics，2007，133（7）：816 －832.

［9］Ching J，Muto M，Beck J L.Structural model updating and health monitoring with incomplete modal data using Gibbs sampler［J］.Computer － Aided Civil and Infrastructure Engineering，2006，21：242 －257.

［10］易偉建，吳高烈，徐 麗.模態參數不確定性分析的貝葉斯方法研究［J］.計算力學學報，2006，26（6）：700 －705.

［11］ Scott M L.Introduction to applied bayesian statistics and estimation for social scientists［M］.Springer，2007.

［12］任宜春，張杰峰，易偉建.基于改進L-P小波的時變模態參數識別方法［J］.振動與沖擊，2009，28（3）：144 －148.