基于馬爾科夫鏈種群競爭的貝葉斯有限元模型修正

關(guān)鍵詞：模型修正；貝葉斯估計；馬爾科夫鏈蒙特卡洛；種群競爭

0 引言

在結(jié)構(gòu)健康監(jiān)測領(lǐng)域中，隨著建立精確有限元模型評估結(jié)構(gòu)服役狀態(tài)的方法不斷完善，有限元法逐漸成為最廣泛、最常用、最基本的數(shù)值分析手段。但由于實際工程結(jié)構(gòu)存在材料的參數(shù)誤差、邊界條件的不明確、外部荷載的變化、環(huán)境因素以及測量噪聲等影響，造成實際工程問題具有來源未知、難以測量的不確定性［1-2］。因此，已有的一些確定性模型修正方法往往處于理論效果好，但實際應(yīng)用效果不佳的尷尬位置。為了解決不確定性造成的理論與實際之間的差距，采用概率或統(tǒng)計手段量化與結(jié)構(gòu)相關(guān)的不確定性，已逐漸成為建立有限元模型關(guān)注的重點。其中，數(shù)學(xué)統(tǒng)計分析中貝葉斯學(xué)派以可充分考慮不確定性而聞名。

近年來，大量國內(nèi)外學(xué)者開始研究基于貝葉斯理論的概率模型修正方法，原因在于貝葉斯推理過程具有可同時考慮結(jié)構(gòu)的先驗信息和當(dāng)前測量值的優(yōu)勢［3］。因此，其逐步成為模型修正領(lǐng)域的重點研究方向，實現(xiàn)了多種工程結(jié)構(gòu)的模型修正和健康監(jiān)測。其中BECK等［4］首次將貝葉斯理論應(yīng)用于工程結(jié)構(gòu)模型修正，并建立了修正的基本流程體系，通過Laplace漸進法近似得到后驗概率分布對兩自由度的框架修正；同時，KATAFYGIOTIS 等［5］ 提出了一種基于Metropolis-Hastings（MH）隨機游走算法的馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo， MCMC）模擬方法，通過模擬采樣近似得到后驗概率分布，并以兩自由度模型驗證了方法的正確性；易偉建等［6］采用MH法對混凝土框架結(jié)構(gòu)進行修正，編制了MCMC方法的損傷識別程序；WAN 等［7］提出一種延緩拒絕自適應(yīng)Metropolis（Delayed Rejection Adaptive Metropolis， DRAM）算法，并對一座復(fù)雜的斜拉人行橋修正，結(jié)果表明DRAM算法擁有較好的樣本遍歷性；GREEN［8］針對非線性問題，提出了一種模擬退火思路的新式MCMC模擬方法，有效地提高了模型修正的效率；CHING等［9］從模型類的選擇出發(fā)，采取避開復(fù)雜概率密度函數(shù)的采樣方法，提出了一種過渡馬爾科夫鏈蒙特卡洛（TransitionalMarkov Chain Monte Carlo，TMCMC）算法，通過3 個例子說明了TMCMC方法在貝葉斯模型更新、模型類選擇和模型平均方面的有效性；ZENG等［10］利用振動響應(yīng)或模態(tài)參數(shù)來估計結(jié)構(gòu)參數(shù)和估計本身相關(guān)的不確定性，采用差分進化自適應(yīng)Metropolis （DifferentialEvolutionary Adaptive Metropolis， DREAM）算法通過并行運行多個Markov鏈抽取樣本提高接受率，并通過十層剪力建筑和三層框架結(jié)構(gòu)的數(shù)值示例驗證了算法的有效性；WU 等［11］為了解決Expectation?Maximization（EM）算法因時間復(fù)雜導(dǎo)致的計算困難，以MCMC采樣用潛在變量的模擬樣本逼近EM算法，提出了MCMCEM算法，并驗證了其比傳統(tǒng)EM 算法的先進性；HUANG等［12］將貝葉斯概率法和微擾模型相結(jié)合的方法，該方法可僅利用少量有噪聲干擾的節(jié)點振動響應(yīng)，對結(jié)構(gòu)進行健康監(jiān)測；劉綱等［13］針對MCMC方法面對高維待修正參數(shù)收斂難效率低問題，引進相關(guān)向量機作為回歸模型代替有限元計算提高效率；SHERRI等［14］針對不確定性量化問題提出一種改進的進化MCMC方法，并使用斯諾克更新器改進和擴展更新過程，通過一個F型框架結(jié)構(gòu)驗證了方法的有效性；彭珍瑞等［15］通過新鳥巢更新思想改進標(biāo)準(zhǔn)MH算法，同時使用支持向量機（Support Vector Machines，SVM）建立代理模型，有效地解決了MH算法不收斂、拒絕率高的問題，并在一個三自由度平面桁架上驗證了方法的可行性。

由上述研究可以看出，貝葉斯理論的核心內(nèi)容主要是通過查閱資料或根據(jù)已有先前經(jīng)驗得到先驗分布，采用當(dāng)前實測信息來修正先驗分布，使其不斷接近后驗概率分布，從而有效處理參數(shù)不確定性。但是后驗分布的求解難度較高，原因在于其分母涉及復(fù)雜且難以求解的高維積分，通過常規(guī)的數(shù)值積分手段效果甚微，引進MCMC采樣算法，以近似解等效替代精確解，因此采樣算法的精確度將直接影響修正結(jié)果。目前已有的標(biāo)準(zhǔn)MH采樣算法，對單個參數(shù)的模型采樣效果較好，在多參數(shù)的采樣過程中容易出現(xiàn)大面積滯留現(xiàn)象，同時樣本采樣拒絕率過高會導(dǎo)致收斂困難，甚至無法收斂。針對以上缺陷，本文在標(biāo)準(zhǔn)MH算法基礎(chǔ)上引入種群競爭算法和差進化算法，可有效地改善高維參數(shù)模型修正的不足，并通過一個桁架結(jié)構(gòu)數(shù)值算例驗證了本文算法的高效性。

1 貝葉斯有限元模型修正理論

采用貝葉斯方法描述有限元模型修正過程中，不確定性參數(shù)分布通過后驗概率密度函數(shù)的數(shù)值求解得到，其根據(jù)貝葉斯方法具體表述為［16］［17］115-132

π （θ | D， M ） = [ P （D| θ，M ） π （θ，M ） ] / [ P （D，M ） ]（1）式中，M 為系統(tǒng)的結(jié)構(gòu)模型類別（每個結(jié)構(gòu)類別的定義的修正參數(shù)不同）；θ 為不確定參數(shù)；D 為結(jié)構(gòu)體系的實測數(shù)據(jù)樣本（本文采用固有頻率和振型）；π（θ，M）為先驗概率密度分布，表示不確定參數(shù)θ 在給定的模型M 中且不含D 時的初始認(rèn)知；P（D|θ，M）為似然函數(shù)，表示實測數(shù)據(jù)與有限元解之間的差異；π（θ|D，M）為給定模型M 和實測數(shù)據(jù)D 的修正參數(shù)的后驗概率分布函數(shù)；由于本文只研究單個模型，后續(xù)公式推導(dǎo)中將M 省略，且P（D，M）與修正參數(shù)θ 不相關(guān)，可用常數(shù)c 表示。因此，π（θ|D）與P（D|θ）π（θ）存在正比例關(guān)系為

式中，N 為用于修正的結(jié)構(gòu)響應(yīng)量階數(shù)；θ 為修正參數(shù)組集的向量（θ∈RN）；Z*為實測獲得的結(jié)構(gòu)響應(yīng)量組集的向量（Z*∈RN）；Z（θ）為有限元模型計算值組集的向量［Z（θ）∈RN］；ε為測試波動、環(huán)境變化等不確定性引起的差值組集的向量（ε∈RN）；其中，RN為對應(yīng)階數(shù)的實數(shù)集范圍，R 為實數(shù)集，上標(biāo)N 為用于修正的結(jié)構(gòu)響應(yīng)量階數(shù)。故后驗概率分布函數(shù)表示為

貝葉斯修正過程知之非難，行之不易，要想取得良好效果需解決以下兩個難題：①結(jié)構(gòu)的復(fù)雜性越大，后驗概率密度函數(shù)解析解的獲取越困難；②未知參數(shù)的維數(shù)以及樣本空間的復(fù)雜性會顯著影響修正過程，且不確定參數(shù)的全局解不容易獲得。因此，貝葉斯修正過程往往與MCMC方法相結(jié)合，通過MCMC方法以采樣代替復(fù)雜高維積分計算，近似獲取π（θ|D），從而得到近似的后驗分布函數(shù)。

2 基于種群的多鏈差分競爭算法

傳統(tǒng)MCMC方法只能得到單條Markov鏈，如應(yīng)用廣泛的標(biāo)準(zhǔn)MH 算法通過隨機游走采樣生成一條Markov鏈，依據(jù)Markov鏈中樣本分布近似得到后驗分布，而其采樣機制的受限于固定的建議分布方差，若方差較大則易導(dǎo)致建議點時常被拒絕沿用上次迭代值，遍歷性差出現(xiàn)滯留現(xiàn)象；方差較小則遍歷步距小、收斂速度慢，甚至出現(xiàn)迭代結(jié)束而樣本未更新，導(dǎo)致該方法在面對復(fù)雜結(jié)構(gòu)或后驗概率密度復(fù)雜的情況時，存在采樣效率低且精度不佳的缺陷。本節(jié)通過建立Markov鏈種群，以多鏈角度引入差分算法以克服滯留和難收斂缺陷。在此基礎(chǔ)上，加入競爭算法，通過種群中多鏈差分競爭的方式選擇更優(yōu)修正方案，提出一種基于種群內(nèi)部競爭的算法以提高采樣效率及修正精度，從而高效逼近目標(biāo)分布。

2. 1 Markov鏈種群

種群競爭需要在多條Markov鏈中進行，本文采用標(biāo)準(zhǔn)MH算法，對同一目標(biāo)函數(shù)設(shè)定不同建議分布方差，采樣得到多條Markov鏈，建立Markov鏈種群［Y1，Y2，…，YN］。在理論層面Markov鏈中N 越大種群攜帶測試信息越豐富效果越好，但綜合考慮種群競爭計算性價比本文設(shè)定N=7，其原因在于剛好能滿足完成一次完整的種群競爭，而交叉、變異采用的Markov鏈不重復(fù)使用，建議分布方差在范圍［0. 1，1. 0］依據(jù)Markov鏈數(shù)均勻取值。

同時，本方法建立Markov鏈種群的方法不局限于標(biāo)準(zhǔn)MH算法，其他能生成有效Markov鏈的采樣算法同樣適用。建立Markov鏈種群的具體步驟如下：

2. 2 多鏈差分算法

多鏈差分算法類似遺傳進化均包含變異、交叉和選擇3個操作運算，但與遺傳進化不同的是變異這一操作運算通過利用種群中個體與個體之間的差量實現(xiàn)條件分布的更新［18］。本文中個體為已構(gòu)建的種群中的Markov鏈，由于構(gòu)建每條鏈的目標(biāo)函數(shù)相同，但每條單鏈的初始建議分布，鏈內(nèi)接受概率不同，以及多鏈差分過程中鏈間變異、交叉概率的不同，導(dǎo)致部分交換有助于進一步改善Markov鏈的質(zhì)量，且多鏈相對單鏈包含更多模態(tài)空間信息，通過選擇操作運算接受可提升部分，可有效改善傳統(tǒng)單鏈面對高維模態(tài)空間能力受限問題。其中多鏈差分算法變異具體操作方式如下：

因此，以Markov 鏈為個體的差分算法基于進化MCMC方法，利用多鏈之間的相互作用關(guān)系得到新的優(yōu)化建議，可解決高維參數(shù)模型修正過程中的滯留缺陷，提高修正精度和效率。

2. 3 種群差分競爭算法

種群競爭算法是一種激勵算法，在多鏈差分算法基礎(chǔ)上，采用不斷競爭與學(xué)習(xí)機制提高差分進化效率，充分利用種群中Markov 鏈攜帶信息，采用較少的Markov鏈獲得較高的精度，是一種基于競爭尋找解決目標(biāo)函數(shù)最優(yōu)方案的新型改進差分算法。

競爭的具體思路是：Markov鏈在多鏈差分算法變異，交叉操作的基礎(chǔ)上融入競爭，在選擇操作中定義處于拒絕狀態(tài)的為失敗者，反之為勝利者；通過競爭算法提升多鏈差分的質(zhì)量，即失敗者通過選取勝利者Markov鏈的部分，試圖學(xué)習(xí)勝利者的優(yōu)勢進而縮小與目標(biāo)函數(shù)的差距得到最優(yōu)方案；同時這又導(dǎo)致勝利者優(yōu)勢只是暫時的領(lǐng)先，多鏈差分算法中的變異、交叉、選擇3項操作也反過來輔助競爭算法，使勝利者繼續(xù)突破尋找更佳的解，失敗者在快速學(xué)習(xí)勝利者的基礎(chǔ)上仍有超過勝利者的可能，勝利者與失敗者的定位在激烈的競爭中隨時可能交換。

由于種群中Markov鏈的前部分處于燃燒階段，如一開始便加入競爭會影響算法的收斂速度，因此競爭算法是人為假設(shè)在m 次迭代后出現(xiàn)競爭直到迭代結(jié)束的整個過程，因此迭代階段按該算法思路，可根據(jù)是否含有競爭劃分為競爭期和正常期（競爭期多鏈差分算法和競爭算法均存在，正常期只含有多鏈差分算法）。

其中，正常期算法具體思路與操作步驟如下：步驟1：開始t=1～m 時，迭代處于正常期，待修正參數(shù)更新方法是在Markov鏈種群中隨機一個鏈記為Yts，根據(jù)式（11）突變得到Y(jié)tsp。

此外由于Markov鏈燃燒階段的樣本是未收斂的，若以此為修正模型參數(shù)會降低精度，根據(jù)經(jīng)驗公式［17］115-132可知，只要迭代步數(shù)滿足T ≥ 300d，得到的樣本個數(shù)對于修正后參數(shù)值的精確幾乎不造成影響，因此，本文迭代步數(shù)設(shè)定為T = 1 × 104 gt; 300d，迭代步數(shù)充分滿足修正要求，其中，正常期與競爭期的m 雖是人為假設(shè)，但其選值需大于Markov鏈?zhǔn)諗糠纸琰c10%的迭代步數(shù)，若取值太小，在Markov鏈?zhǔn)諗糠纸琰c之內(nèi)，將會影響算法的收斂；若取值太大，超過收斂分界點太多，將會消除競爭作用，延長正常期，影響修正精度和效率。所提方法的具體步驟流程如圖1所示。

3 數(shù)值算例

3. 1 平面桁架

某平面桁架結(jié)構(gòu)如圖2所示，以此算例檢驗所提算法的高效性。該桁架為對稱斜三角桁架，由21個桿件鉸接組成，共有12個節(jié)點、21個單元和21個自由度，如圖3所示。按桿件單元位置將其劃分為5個小組，其中編號1、2、3、10、11號單元為第1組；編號1、4、12、13、18號單元為第2組；編號13、14、15、19、21號單元為第3組；編號5、6、16、19、20號單元為第4組；編號6、7、8、9、17號單元為第5組。選取第1、3、5組的彈性模量和密度為模型的待修正參數(shù)，一共6個待修正參數(shù)，分別為θ1 =E1 /E0，θ2 =ρ2 /ρ0，θ3 =E3 /E0，θ4=ρ4 /ρ0，θ5 =E1 /E0，θ6=ρ6 /ρ0，其中桁架基數(shù)分別為E0 =210 GPa，ρ0 =7 800 kg/m3。桁架組號劃分內(nèi)容及結(jié)構(gòu)參數(shù)如表1所示。

本文假定6 個試驗?zāi)Ｐ蛥?shù)服從正態(tài)分布N｛［0. 8，0. 12］，［1. 0，0. 12］，［1. 2，0. 12］，［1. 4，0. 12］，［1. 6，0. 12］，［1. 8，0. 12］｝，且參數(shù)之間相互獨立。通過概率統(tǒng)計手段，隨機抽取一組待修正參數(shù)代入桁架有限元模型，計算得到前6階模態(tài)，為了充分考慮實際工程測試中的不確定性因素，得到更接近實測情況的測試數(shù)據(jù)，本文共模擬計算35次，以得到的35次數(shù)據(jù)均值為結(jié)構(gòu)實測響應(yīng)。

對樣本進行獨立檢驗，由圖4可知，6個參數(shù)35次模擬得到的樣本點絕大部分處于正態(tài)分布的假設(shè)概率密度山形圖的山峰表面，表明隨機采樣得到樣本點滿足試驗獨立假設(shè)要求，所得到的樣本點具有足夠的代表性。

將不確定參數(shù)θ1～θ6均除以其對應(yīng)的初始值θ0k （k=1、2、3、4、5、6），得到歸一化數(shù)據(jù)后繪制置信橢圓圖（誤差橢圓圖）如圖5所示。圖中的點為迭代1×104次后得到的不確定性參數(shù)樣本，橢圓區(qū)域表示95%置信度的樣本空間，可以看出不確定性參數(shù)基本落位在置信橢圓之內(nèi)。

圖6給出了標(biāo)準(zhǔn)MH算法和本文所提算法的收斂曲線對比圖。其中，誤差定義為6個待修正參數(shù)的絕對誤差之和，其計算式為

如圖6所示，在收斂速度上，可明顯觀察到本文所提算法在第200 次迭代步后基本收斂，而MH算法在400次迭代步才逐漸收斂；在收斂精度上，也可看出本文所提算法的誤差更小；在收斂路徑上，MH算法呈鋸齒狀，表明采樣空間未充分探索，而本文所提算法小范圍波動，具備良好的遍歷性。

圖7給出了MH算法和本文所提算法的修正對比圖。由圖7分析可知，當(dāng)完成1×104迭代步時兩種算法均已收斂，但MH算法得到的Markov鏈6個參數(shù)均明顯出現(xiàn)大塊區(qū)域的滯留現(xiàn)象，在較高維模型修正中樣本的采樣效率極低；而本文所提算法融入差分思想有效改善高維參數(shù)模型修正滯留缺點，同時收斂速度較快。再根據(jù)圖7（a）、圖7（b）對參數(shù)θ4和θ5局部對比放大可知（圖7右側(cè)圖），本文所提算法迭代得到的參數(shù)樣本滯留區(qū)域較小、接受率和采樣效率較高、采樣具有良好的遍歷性，且融入的種群競爭和差分步驟使修正在穩(wěn)定后仍在收斂值附近小幅度起伏波動，進而提高模型修正精度，相比MH算法可更好地模擬出接近真實情況的修正參數(shù)后驗分布的樣本空間。

提取前800次迭代步的修正過程如圖8所示。處于正常期時，以Markov鏈為個體的差分算法有效地改善了傳統(tǒng)MCMC采樣效率低且難收斂的問題；200次迭代步加入競爭算法進入競爭期后，差分算法與競爭算法相互提供激勵條件，可觀察到6個修正參數(shù)收斂后的波動幅度更小，進一步提升了模型修正精度。

理論上對比修正后模型得到的響應(yīng)與實測響應(yīng)是否一致可判斷是否達到修正要求，但在研究中僅以此來判斷模型修正效果是片面的，還需結(jié)合修正參數(shù)值。因此，本文給出了不同算法修正前后模態(tài)響應(yīng)和參數(shù)的對比數(shù)據(jù)表（表2、表3）。從表中可以看出：從頻率響應(yīng)角度分析，首先，初始模型響應(yīng)與真實響應(yīng)的前6階頻率均存有5% 左右的誤差，修正前第2 階誤差最大，達到10. 58%。標(biāo)準(zhǔn)MH算法修正后前6階頻率誤差均值為1%左右，修正后最大誤差為1. 77%，而本文所提算法修正后頻率誤差均值在0. 21%以內(nèi)，且修正后最大誤差僅為0. 21%；從修正前后參數(shù)對比分析，MH算法表現(xiàn)遠不如本文所提算法，其中θ3和θ5修正后相對誤差分別為5. 84%和5. 1%，而本文所提算法的各參數(shù)相對誤差均在1%以內(nèi)。對比分析表明本文所提算法識別結(jié)構(gòu)參數(shù)精度較高，對高維參數(shù)模型修正效果遠勝于標(biāo)準(zhǔn)MH算法。

3. 2 加噪試驗

為了使桁架模型修正更加合乎實際，以及進一步驗證所提方法的高效性，本文對試驗得到的測試數(shù)據(jù)加入10%的噪聲，用來模擬實際工程問題中測試遇到的不確定性因素（如實測人員操作規(guī)范問題、設(shè)備裝置精度問題等）。為對比不同算法在同一噪聲條件下的修正效果，根據(jù)式（26）隨機加入10%的測試噪聲，且假定模態(tài)噪聲在不同階和不同參數(shù)之間相互獨立。

圖9（b）為Markov鏈和每個修正參數(shù)的概率直方圖。首先對采樣過程圖分析可知，加噪10%對本文所提算法的收斂性影響不大，Markov鏈均可收斂，且每個參數(shù)的概率直方圖表明，加噪10%后的參數(shù)也大致在真實值±0. 04E0/ρ0左右，雖略有差距但在可允許接受范圍之內(nèi)，表明了本文所提算法具有良好的魯棒性。對加噪后得到的修正參數(shù)實施正態(tài)概率檢驗，得到結(jié)果如圖10所示。從圖10可知，加噪后計算得到的樣本點基本位于假設(shè)參數(shù)分布線上，表明采樣得到的樣本點滿足正態(tài)分布假設(shè)，結(jié)果可充分接受（對縱坐標(biāo)刻度進行了調(diào)整，以便觀察正態(tài)率，不影響檢驗數(shù)據(jù)是否近似服從正態(tài)分布）。

因此，結(jié)合上述分析可知，本文所提算法可以在考慮各種不確定性因素下，可精準(zhǔn)識別結(jié)構(gòu)模型參數(shù)變化，高效地完成有限元模型修正，且算法具備較高的修正精度和良好的魯棒性。

4 結(jié)論

針對不確定性有限元模型修正問題，以貝葉斯統(tǒng)計理論為核心，在標(biāo)準(zhǔn)MH算法基礎(chǔ)上引進了競爭算法和差分算法，提出了一種基于Markov鏈種群競爭差分進化采樣的改進算法，通過一個六參數(shù)桁架結(jié)構(gòu)模型修正數(shù)值算例驗證了所提算法的精度和效率，且面對10%的測試噪聲干擾下仍保持一定精度，并與傳統(tǒng)MH算法計算結(jié)果對比，得出以下結(jié)論：

1）本文所提算法面對高維參數(shù)模型修正采樣接受率、遍歷性更好，通過提高采樣效率的方式極大地提升了模型修正速率，實現(xiàn)了高維參數(shù)的不確定性有限元模型修正。

2）本文所提算法通過種群競爭和差分手段，以Markov鏈之間攜帶信息差異為建議優(yōu)化方向，進一步提高了算法對高維參數(shù)模型修正精度。

3）本文所提算法在10%的測試噪聲干擾下，仍能保證較高精度，得到可靠的識別結(jié)果，表明所提算法具備良好的魯棒性。

本研究目前使用數(shù)值算例驗證了所提算法的高效性，在保證采樣效率的同時具有較高的精度和穩(wěn)定性。后續(xù)將陸續(xù)開展結(jié)構(gòu)模態(tài)試驗研究，為實際結(jié)構(gòu)中考慮不確定性的有限元模型修正提供一種有效穩(wěn)定的手段。