橋梁結構多測組振型融合的兩階段快速貝葉斯方法研究

劉玲玲，顏王吉，李 丹，任偉新，秦 超

(1.合肥工業大學土木與水利工程學院，合肥230009；2.澳門大學智慧城市物聯網國家重點實驗室，澳門999078；3.澳門大學土木與環境工程系，澳門999078 4.深圳大學土木與交通工程學院，廣東深圳518060)

模態參數識別是結構狀態評估、有限元模型修正和損傷識別的基礎[1]。環境激勵下的模態參數識別方法具有無需貴重的激勵設備、不影響結構的正常使用等特點，被廣泛應用于工程實踐中，主要分為時域法、頻域法和時頻域法等。其中，時域法主要有隨機子空間法（SSI）、隨機減量技術與特征系統實現算法的聯合方法[2]和極點對稱模態分解方法[3]等，頻域法有頻域分解法和振動響應傳遞比方法[4]等。

自然環境條件下橋梁結構的動力響應測試數據具有幅值小和隨機性強的特點。與測試噪聲及計算模型等因素有關的不確定性影響廣泛存在于參數識別問題中[5-7]。因此，有必要在模態參數識別過程中考慮不確定性的影響。近年來，基于概率統計分析原理的貝葉斯理論在模態參數識別領域發揮了重要作用[8]。例如，Yuen和Katafygiotis[9]提出了基于貝葉斯模態參數識別的FFT 驅動方法，將模態域的結構模型嵌入到貝葉斯系統識別框架中，可識別出頻率、阻尼比、振型等參數的最優值并量化其不確定性。貝葉斯FFT方法具有嚴格的理論推導，然而，其需要識別的參數過多，計算量巨大，不能直接應用于工程實際。為解決這一問題，Au[10]提出了快速貝葉斯FFT 方法，利用矩陣論和函數優化理論相關方法極大地降低了貝葉斯FFT方法對于優化和計算協方差矩陣的難度；Li 等[11]在貝葉斯FFT 方法中引入期望最大化算法，提高了其計算效率和魯棒性。

振型可以提供橋梁結構變化的空間信息，在橋梁結構損傷識別和模型修正中具有重要價值。然而，在橋梁結構全尺寸環境振動實驗中，因傳感器數量不足以覆蓋整個橋梁結構等原因，常需劃分多個測組進行多次測試，如圖1 所示。無法直接獲得并處理全橋實測數據從而識別橋梁的整體振型。如何從各測組數據中準確地識別頻率、阻尼比、局部振型，并利用局部振型準確地獲得橋梁整體振型至關重要?！熬植孔钚《朔ā迸c“整體最小二乘法”[12]等方法可快速獲得整體振型，但都無法量化整體振型的不確定性。隨后，Au等[13]提出了一種快速貝葉斯模態參數識別方法，通過融合各測組的FFT 信息求解整體振型。在Au 等[13]的基礎上，Zhu 等14]引入期望最大化算法提高了振型最優值的計算效率，可適用于多模態甚至密集模態的情況。文獻[15]在Au等[13]的基礎上推導了Hessian矩陣的解析表達式，可進一步量化整體振型的后驗不確定性。

針對橋梁多測組環境振動實驗，本文提出了多次測試條件下橋梁振型融合的兩階段快速貝葉斯方法。首先，采用快速貝葉斯FFT 方法[10]求解各測組頻率、阻尼比和局部振型的最優值及協方差。然后，基于貝葉斯原理[16-17]融合局部振型統計信息，通過解析迭代優化算法快速求解整體振型的最佳估計，并可計算負對數似然函數對整體振型的Hessian矩陣。該方法可利用局部振型的不確定性考慮各測組的數據質量的影響，并量化整體振型的不確定性。為驗證該方法在實際橋梁結構模態分析中的有效性，對一座斜拉橋的實測數據進行了分析處理。同時，論文引入模態置信準則（MAC）概率指標[18]，評估了整體振型的后驗不確定性，并研究了各測組數據質量及頻帶寬度等因素對整體振型不確定性的影響規律。研究表明，采用該方法處理實橋數據效率高，且能合理地量化識別結果的不確定性。

1 多測組工況下的橋梁模態參數識別及振型融合問題

假設在橋梁結構環境振動實驗中，共有nt個測組，測組i包含ni個測點，實驗共nl個測點，每個測組至少有一個參考點與其他任意測組共享，則三者關系滿足如圖1所示，采集橋梁環境振動實驗數據。橋梁模態參數識別及振型融合不確定性量化可分兩個階段實現：

圖1 多測組環境振動實驗測點布置及兩階段快速貝葉斯方法模態分析步驟示意圖

(1)第1 階段：運用快速貝葉斯FFT 方法識別各測組模態參數最優值及其協方差。識別參數包括：模態頻率f、阻尼比ζ、局部振型ψ、模態激勵的功率譜S及預測誤差的功率譜σ2。此外，令∈Rni和Cψr,i分別表示測組i的第r階模態局部振型最優值及其協方差，上標⌒表示參數最優值，并令m表示識別模態的階數。

(2)第2階段：運用振型融合的快速貝葉斯方法求解整體振型最優值及其協方差。令φr∈Rnl表示包含所有測點的第r階模態整體振型，并令局部振型φr,i∈Rni表示φr的組成部分，對應測組i。

2 第1階段：快速貝葉斯FFT方法

2.1 理論基礎

假定一離散為nd個自由度的線性系統在環境振動試驗中共有nl( ≤nd)個測點，分為nt個測組進行測量，每組共有ni個測點。對第i個測組加速度響應∈Rni建?？傻茫?/p>

其中：xj是第i組模型加速度響應，ej∈Rnl(j=1,2,…Ns)為預測誤差，Ns為每個通道的采樣數。待識別的模態參數為λθ={f,ζ,ψ,S,σ2}。xj的傅里葉變換為[10]

其中：i2=-1；Δt為采樣間隔；Fk=Re ?k，Gk=Im ?k分別為的傅里葉變換?k的實部和虛部；k=2,3,…,Nq，對應的頻率為fk=(k-1)/NsΔt，其中Nq=int[Ns/2]+1。

對于分離模態，假定在所選的頻帶內僅包含一階模態。根據貝葉斯定理：

其中：c0是歸一化常數為先驗概率密度函數（PDF），在先驗信息未知時，通常假定其為均勻分布[10,23]。由于該先驗分布在后面推導的負對數似然函數式（8）中體現為一個常數，因而該假定并不影響整體振型最優值和協方差的求解。此時，參數后驗PDF為[10]

其中：Ck是Zk的協方差矩陣，其表達式為

式中：ψ∈Rni×m為模態振型矩陣；I2ni∈R2ni×2ni為單位矩陣；Hk為模態響應的譜密度矩陣，其元素表達式為

式中：βik=f(i)/fk為頻率比，f(i)、ζi分別為第i階模態的頻率和阻尼比，Sij為第i個和第j個模態激勵間的互功率譜密度。

由式（4）可進一步推得后驗PDF 與負對數似然函數成正比，即：

其中：L(λθ)為目標函數[10]：

式中：Nf為k的總個數其中：βk=f/fk；

2.2 參數的最優值及其后驗不確定性

對式（9）中矩陣A進行特征值分解，其最大特征值所對應的特征向量即為振型ψ的最優值，進而求出歸一化的振型最優值/‖ψ‖。

將式（8）關于ψ最小化[10]：

后驗協方差矩陣等價于負對數似然函數L(λθ)的Hessian矩陣的逆矩陣[10]。局部振型的協方差矩陣Cψr,i為中對應于振型ψ的ni×ni階方陣[19]。

通過上述快速貝葉斯FFT方法可求出各測組頻率、阻尼比、局部振型的最優值及其協方差矩陣。假定任意測組局部振型的PDF 近似服從高斯分布[16,20]，則對于局部振型φr,i∈Rni，其概率分布公式為

其中：N為正態分布符號。

3 第2 階段：快速貝葉斯振型融合方法

3.1 局部振型與整體振型的關聯

其中：Li∈Rni×nl，若第i個測組的第p個測點對φr的第q個自由度有貢獻，則元素Li(p,q)=1，否則為0。

識別的局部振型已對其歐幾里得范數歸一化，同樣地，將φr,i對其歐幾里得范數歸一化：

在2.2 節中已說明，測組i的局部振型可通過高斯分布很好地近似，則以φr表示并反映(,Cψr,i)貢獻的似然函數為[16,20]

在式（14）中已省略與φr獨立的歸一化常數。

3.2 整體振型后驗概率密度函數

根據貝葉斯定理，整體振型φr的后驗PDF如下：

其中：?={,Cψr,i:i=1,…,nt}，c1是歸一化常數；p(φr)為先驗PDF。假定從不同測組中識別的局部振型在統計上是獨立的，因而反映?貢獻的似然函數如下[16]：

假定p(φr)為均勻分布，φr的后驗PDF 滿足p(φr| ?)∝exp(-Las(φr))，Las為目標函數[16]：

式（17）中的整體振型φr受到約束條件φrTφr=1 的限制。

當局部振型協方差矩陣為單位矩陣I∈Rni×ni時，每個測組分配的權重相同，此時，目標函數簡化為

“整體最小二乘法”[12]是通過已識別的局部振型和待識別整體振型之間的平方差構造目標函數，為所有測組分配相同的權重，其目標函數與式（18）一致。當假定局部振型協方差矩陣為I∈Rni×ni時，即不考慮局部振型的協方差時，該表達式退化為“整體最小二乘法”。因此，快速貝葉斯振型融合方法是對整體最小二乘法的延伸，它可利用局部振型協方差考慮不同測組數據的不確定性信息。

3.3 最優值與協方差矩陣的求解

最小化式（17）即可求解φr的最優值。引入輔助變量其中i=1,2,…,nt，在的約束條件下，利用拉格朗日乘數法，式（17）可變換為[16]

其中：γr、βr,i(i=1,2,…,nt)均為拉格朗日乘子，分別使得約束條件‖φr‖2=1和成立。

3.3.1 整體振型最優值

式（19）中待識別的參數為λas={φr,γr,βr,i,χr,i:i=1,2,…,nt}。直接求解參數最優值的計算過程具有高維及非線性特征，因而采用迭代求解策略計算。分別求出χr,i、βr,i、φr和γr最優值的解析表達式，然后通過線性迭代優化算法確定其最優值[16]。

(1)χr,i和βr,i最優值：

其中：sgn( ? )表示符號函數。

(2)φr和γr的最優值：

其中：

3.3.2 整體振型后驗不確定性

整體振型φr的不確定性近似服從以參數最優值為中心的高斯分布。將目標函數式（17）中的整體振型歸一化，即可得：

令表示式（24）在最優值處求得的Hessian 矩陣，其表達式為[17]:

為整體評估最優整體振型的后驗不確定性，本文引入MAC概率指標[18]，用來衡量振型變量φr與振型最優值之間的相似性。當φr的不確定性越小時，φr越接近于，計算的MAC 值越接近于1，為計算方便，選擇計算φr的期望及之間的MAC 值。該指標的近似解析表達式為[17–18]

4 計算流程

針對橋梁結構分組多次測試的環境振動試驗，兩階段快速貝葉斯方法的計算過程如圖2所示。

圖2 2階段快速貝葉斯方法的計算流程

5 實橋分析

常德白馬湖公園虹橋是一座4跨曲線鋼箱梁斜拉人行橋，全長120 m，跨徑布置為27 m+33 m+33 m+27 m，如圖3 所示。該橋的環境振動實驗分為8個測組，每組共有12個測點。其中豎向測點有6個，布置在橋梁兩側；橫向水平測點有3 個；全橋共設置3 個參考點，布置在中跨跨中附近偏向中墩側的斷面處。全橋加速度傳感器布置情況如圖4 所示。分組情況見表1。實驗采樣頻率為500 Hz，采樣時長為15 min。實測的詳細方案參見文獻[21–22]。

表1 測點分組布置

圖3 虹橋實景圖

圖4 橋面測點布置圖[22]

采用快速貝葉斯FFT方法識別模態參數時需要先選取頻帶。假定在所選頻帶中僅包含1 階分離模態，并以頻率初始值為中心選取頻帶，范圍為f0(1±κζ0)[10,23]，此處f0和ζ0分別為頻率和阻尼比的初始值，κ為頻帶寬度[23]。如圖5所示。

圖5 第2測組功率譜圖

本文選取各階模態頻率初始值f0分別為3.784 Hz、4.761 Hz、5.981 Hz、6.470 Hz，設各階模態頻帶寬度κ為2，并取各階模態初始阻尼比ζ0為0.01。首先求解振型ψ的最優值，并通過解決一個數值優化問題來確定f、ζ、S和σ2的最優值，然后計算參數的后驗協方差矩陣，并從中提取局部振型協方差矩陣Cψr,i。利用SSI 和快速貝葉斯FFT 方法識別的前4階豎向模態頻率及阻尼比值如表2所示。2種方法識別的頻率值基本一致，阻尼比值也很接近，由此可見，采用快速貝葉斯FFT 方法識別頻率和阻尼比的精度較高。第5、7列為采用快速貝葉斯FFT方法識別的頻率及阻尼比的變異系數，前4 階豎向模態頻率變異系數最大為0.12%，不確定性很小，而阻尼比變異系數最大達到15%，不確定性較大。在識別得到局部振型最優值及其協方差矩陣后，在參數最優值解析表達式的基礎上，通過線性優化迭代算法快速計算整體振型的最佳估計，并計算負對數似然函數對整體振型的Hessian 矩陣，以量化整體振型的不確定性。整個計算過程包含2 個階段，每個階段的計算時間均不超過2 min，其中，從單個測組數據中識別每階模態的運算可在幾十秒內完成。

表2 模態識別結果對比

圖6 為采用貝葉斯振型融合方法識別的前4 階豎向模態整體振型圖。以1-E(Mr)評估局部振型及整體振型的不確定性，不同測組局部振型不確定性及整體振型不確定性的評估結果如表3所示。將不同測組的局部振型不確定性變化繪于圖7，其中實線代表不同測組局部振型不確定性，虛線代表整體振型不確定性。局部振型不確定性可表征各測組數據的質量，在振型融合過程中，各測組間相互影響，數據質量差的測組將會干擾數據質量好的測組，因此應被分配較小的權重。由圖7 可見，個別測組的局部振型不確定性較大，而對應的整體振型不確定性有所下降，表明數據質量差的測組被賦予的權重較小，從而有效地抑制了它們帶來的影響。由圖7（c）和圖7（d）可知，相比第4階模態的測組2，快速貝葉斯振型融合方法對第3 階模態的測組3 和6 的識別結果的“抑制作用”更為明顯。其主要原因是，對第3階模態，由測組3和6所識別的局部振型不確定性相對其余測組很大，其被分配權重更小，因而對整體振型識別的影響更小。

表3 各測組局部振型及整體振型不確定性評估（×10-4）

圖6 前4階豎向模態整體振型圖

圖7 前4階豎向模態各測組局部振型不確定性變化圖

為研究頻帶寬度對振型識別結果的影響，設定頻帶寬度參數的變化范圍為κ=1,2,…,6。整體振型不確定性隨頻帶寬度增加的變化趨勢如圖8 所示，虹橋前2 階豎向模態的整體振型不確定性隨寬度的增加而減少。而第3、4階豎向模態的整體振型不確定性為先減少而后有所增加，在頻帶寬度為3時，其不確定性最小。整體振型不確定性在初期隨頻帶寬度增大而減小，主要原因是頻帶寬度越大，數據信息越完備，識別結果的不確定性就越小。而第3階和第4階模態整體振型不確定性在后期隨頻帶寬度增大呈增大的趨勢，主要原因為頻帶越寬，模態激勵和預測誤差的平譜假定及頻帶內單個模態假定等越不滿足，建模誤差越大。因而在橋梁工程實際應用中，需要在兩者之間進行權衡，以選擇最優的頻帶寬度。

圖8 頻帶寬度分析

6 結語

本文提出了橋梁結構多測組振型融合的兩階段快速貝葉斯方法。第1階段采用快速貝葉斯FFT方法求解各測組局部振型的最優值及協方差；第2 階段利用局部振型的協方差，基于快速貝葉斯振型融合方法有效量化整體振型的不確定性。最后通過一座斜拉橋的環境振動實測數據驗證了該方法的有效性。研究所得結論主要包括：

(1)采用快速貝葉斯FFT 方法可以快速識別頻率、阻尼比等模態參數的不確定性。實橋分析表明，頻率變異系數最大為0.12%，表明其不確定性很??；而阻尼比變異系數可達到15%，表明其不確定性非常明顯。

(2)快速貝葉斯振型融合方法通過解析迭代優化算法快速求解整體振型的最優值，通過負對數似然函數關于整體振型的Hessian矩陣的解析表達式，可實現快速計算并量化整體振型不確定性。該方法利用局部振型協方差矩陣自動分配不同測組的權重，有效地考慮了不同測組數據中包含的不確定性信息。

(3)實橋結果表明，快速貝葉斯振型融合方法在振型融合過程中能有效抑制數據質量較差測組的影響，且對數據質量相對越差的測組，其“抑制作用”越明顯。所識別的整體振型及其不確定性可進一步應用于考慮不確定性影響的橋梁結構損傷識別及模型修正等方面。

(4)實橋分析表明，存在若干階模態的整體振型不確定性隨頻帶寬度增大的變化趨勢為先減小而后增大。整體振型不確定性在初期逐步減小，主要原因為頻帶寬度越大數據信息越完備；整體振型不確定性后期逐漸增大，主要原因為頻帶越寬使得模態激勵和預測誤差的平譜假定及頻帶內僅包含單個模態的假定等越不滿足，導致建模誤差增大。致謝：

感謝任偉新教授課題組參與白馬湖虹橋動力測試工作的所有研究生，感謝李志剛和魏錦輝博士在數據分析中提供的幫助和支持。