基于Gibbs抽樣的馬爾科夫蒙特卡羅方法在結(jié)構(gòu)物理參數(shù)識別及損傷定位中的研究

劉書奎，吳子燕，張玉兵

（西北工業(yè)大學(xué) 力學(xué)與土木建筑學(xué)院，西安 710072）

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

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

本文視結(jié)構(gòu)模態(tài)參數(shù)為隨機變量，根據(jù)結(jié)構(gòu)動力特征方程，將模態(tài)參數(shù)與結(jié)構(gòu)物理參數(shù)通過線性結(jié)構(gòu)識別模型聯(lián)系起來，采用基于Gibbs抽樣的馬爾科夫蒙特卡羅方法來對結(jié)構(gòu)物理參數(shù)進行后驗抽樣，進而根據(jù)后驗樣本的統(tǒng)計特性進行結(jié)構(gòu)物理參數(shù)識別及損傷定位。

1 線性回歸模型

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

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

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

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

當(dāng)α=β=0時，該先驗分布就轉(zhuǎn)變?yōu)橥ǔ５臒o信息先驗分布，即：

由貝葉斯估計理論，模型參數(shù)θ和預(yù)測誤差方差的后驗聯(lián)合分布為：

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

2 線性結(jié)構(gòu)識別模型

在結(jié)構(gòu)健康監(jiān)測領(lǐng)域，通常用線性結(jié)構(gòu)模型來進行模型更新。一般振動測試數(shù)據(jù)都是在很低幅值的激勵下得到的，因此很多結(jié)構(gòu)包括損傷結(jié)構(gòu)都近似地表現(xiàn)出線性行為［9］。

結(jié)構(gòu)第i階模態(tài)的理論頻率和振型向量φi應(yīng)滿足如下特征方程：

其中φi=［φi1，φi2，…，φim］T，m為結(jié)構(gòu)自由度數(shù)目。當(dāng)采用集中質(zhì)量矩陣時，該方程可以表示為：

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

對式（10）做變換，并加入預(yù)測誤差項，可得線性結(jié)構(gòu)識別模型：

對應(yīng)于線性回歸模型Y=Xθ+e。

3 Gibbs抽樣過程

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

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

4 損傷概率

剛度參數(shù)向量θ中的元素與結(jié)構(gòu)構(gòu)件存在一一對應(yīng)關(guān)系，本文利用θ中某一元素出現(xiàn)降低來標(biāo)識損傷的位置，以其降低的比例d來刻畫損傷的程度。

5 算例分析

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

圖1 三層剪切型結(jié)構(gòu)Fig.1 Three-story shear-building

表1 結(jié)構(gòu)一階自振圓頻率統(tǒng)計特性Tab.1 Statistical properties of the first-order natural frequency

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

由圖2可以看出，不同損傷類型，當(dāng)結(jié)構(gòu)某處發(fā)生剛度降低時，對應(yīng)的結(jié)構(gòu)剛度參數(shù)θi的抽樣曲線也會發(fā)生明顯變化，從而可以清楚地判斷出結(jié)構(gòu)的損傷位置及程度。

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

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

圖4 結(jié)構(gòu)在不同損傷類型下對應(yīng)不同剛度降低的超越概率Fig.4 Estimated damage probability cures for all damaged patterns

表2 結(jié)構(gòu)剛度參數(shù)后驗樣本的統(tǒng)計特性Tab.2 Posterior statistical properties of the stiffness parameters

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

6 結(jié)論

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

［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］張 坤，羅紹湘，段忠東.基于動態(tài)響應(yīng)靈敏度同步反演結(jié)構(gòu)物理參數(shù)與輸入的算法［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］李小華.結(jié)構(gòu)參數(shù)識別的貝葉斯方法及應(yīng)用研究［D］.中國地震局工程力學(xué)研究所，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］易偉建，吳高烈，徐 麗.模態(tài)參數(shù)不確定性分析的貝葉斯方法研究［J］.計算力學(xué)學(xué)報，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小波的時變模態(tài)參數(shù)識別方法［J］.振動與沖擊，2009，28（3）：144 －148.