模態識別的Bayesian TDD-FFT法及其應用

吳 杰， 徐洪俊， 張其林

(同濟大學 土木工程學院，上海 200092)

不確定度是評價測量結果質量高低的重要指標。隨著科學技術和生產的發展，不確定度分析受到越來越高的重視[1]。但在土木工程中，大部分的模態參數識別方法沒有考慮識別結果的不確定度問題。Yuen等[2-3]提出了貝葉斯快速傅里葉變換(Bayesian Fast Fourier Transform，Bayesian FFT)方法，通過對原始時域數據進行變換，識別結構的模態特性并求解參數的不確定性。但該方法識別過程存在以下三個問題：①由于后驗的概率密度函數(Probability Density Function, PDF)是模態參數的非線性函數，因此模態參數的最佳估計(Most Probable Value, MPV)需要求解一個多維的數值優化問題，計算過程復雜；②Hessian矩陣需要通過有限差分法來確定，導致計算效率低下；③由于模態參數的數量隨著測試自由度數量的增加而增加，目標函數的求解涉及到一個病態矩陣的求逆問題。這些問題的存在嚴重制約了該方法的應用。

Au[4]針對選定的某一共振頻率帶的單個模態，提出了快速貝葉斯FFT(Fast Bayesian FFT)方法，該方法針對某一共振頻率帶的單個模態，不需要考慮所計算的自由度數量，模態參數的最佳估計可通過一個四維的數值優化問題來確定，而且Hessian矩陣可以通過對數似然函數關于模態參數的二階導數快速準確地計算。

Fast Bayesian FFT方法雖然解決了模態參數的MPV和Hessian病態矩陣的求解問題，但并沒有進一步對模態參數的后驗概率分布進行探討[5-9]。鑒于此，本文首先基于貝葉斯理論和結構響應概率模型，提出了一種將多測點多自由度信號轉換成單測點單自由度信號的Bayesian TDD (Time Domain Decomposition)方法。然后，對單測點單自由度響應下的貝葉斯FFT(Bayesian FFT)方法進行了推導，并利用蒙特卡洛抽樣分析得出各個參數的后驗概率分布、模態參數MPV和不確定性評估，從而避免了Hessian矩陣的求解運算。最后，提出Bayesian TDD-FFT方法，并利用該方法對上海中心大廈實測數據進行分析，探討了評估模態參數的后驗概率分布及聯合分布問題，同時將識別參數結果與Fast Bayesian FFT方法進行了對比分析。

1 理論推導

1.1 Bayesian TDD方法

對于線性結構，一般假設輸入為平穩的零均值的Gauss過程，則結構的輸出為0均值平穩Gauss過程，即：

Xn×1～N(0,C)

(1)

式中:C為Xn×1響應輸出的協方差矩陣。

根據多維高斯分布性質[10], 對于n維的高斯分布Xn×1～N(0,C)，必存在正交陣Φn×n，使得：

X～N(0,ΦΛΦT)

(2)

Yn×1=ΦTX～N(0,Λ),其中Λ=diag(λi)

(3)

即Y的各個變量獨立。令θ={C},P(θ)為病態的均勻分布：P(θ)=1，滿足

P(θ|X)∝P(X|θ)P(θ)

(4)

(5)

(6)

將P(X|θ)寫成對數形式，則有：

P(X|θ)=exp[-f(θ)]

(7)

其中，

(8)

由于ΦTCΦ=Λ,Λ=diag(λi)，則：

(9)

(10)

要使得P(θ|X)取得MPV值，需f(λ1,…,λi,…,λn)取最小值，則：

(11)

(12)

(13)

在實際結構中，往往只關心所測結構的前k階振型，測點的維數n>k，由振型疊加原理可得結構的振動響應為：

(14)

1.2 單測點單自由度的Bayesian FFT

對于單測點單自由度振動響應，Bayesian FFT理論如下：

Yk=Fk+iGk=

(15)

考慮到采樣通道的電壓偏置和傅里葉變換的共軛等原因，僅k=1,2,…,Nq(Nq=int[N/2]+1)所對應的FFT數據被應用于結構的模態參數識別，其所對應的頻率值為fk=(k-1)/NΔt。

有關模態參數θ的后驗概率密度函數P(θ|{Zk})正比于似然函數P({Zk}|θ)：

P(θ|{Zk})∝P({Zk}|θ)=

(16)

式中：Zk=[FkGk]T∈R2，k=1,2,3,…,Nq；Ck為矩陣Zk的協方差矩陣，其表達式為：

(17)

式中：σ2為預測誤差的功率譜密度；I2為2階單位矩陣；Hk為模態響應自功率譜，表達式為：

(18)

(19)

式中：βk=f/fk為頻率比，f和ζ分別是單階模態信號特征頻率和阻尼比；S是單階模態響應的自功率譜。式(17)可寫成：

Ck=[(SDk+σ2)/2]I2

(20)

detCk=2-2(SDk+σ2)2

(21)

若將后驗的PDF進一步寫成關于對數似然函數L(θ)的表達式，則有：

P(θ|{Zk})∝exp[-L(θ)]

(22)

式中：

(23)

將Ck和Zk代入式(23)，則有：

L(θ)=-(Nq-1)ln 2+

(24)

(25)

f、ζ、S和σ2的MPV值可通過對式(25)的無約束優化直接求得。求得的f、ζ、S和σ2的MPV值后，利用蒙特卡羅方法對式(22)進行采樣分析，得到f、ζ、S和σ2的后驗概率分布，進行不確定度分析。

2 計算流程

上文提出的Bayesian TDD-FFT方法的計算流程見圖1所示。

圖1 Bayesian TDD-FFT方法計算流程圖

3 數值驗證

3.1 基本參數

本節通過數值模擬對Bayesian TDD-FFT模態參數識別方法進行驗證。利用中心差分法計算懸臂梁在隨機荷載作下的結構響應數據。懸臂梁詳細參數如下：長度12 m，熱軋工字鋼I10，面積為14.345 cm2，慣性矩為245.0 cm4，Q345鋼材，楊氏模量E=206 GPa，泊松比為0.3，密度為7 850 kg/m3，取前3階模態進行分析，每階模態的阻尼比均取0.5%，前三階頻率分別為0.806 8 Hz，5.056 8 Hz，14.160 6 Hz，荷載由隨機白噪聲生成，并為結構響應添加隨機白噪聲。沿梁長度方向均勻設置12個加速度測點，采集頻率50 Hz，時長1 200 s，加載點位于懸臂端，如圖2所示。圖3顯示了測點Ac3的加速度時程曲線。

圖2 懸臂梁測點布置

3.2 Bayesian TDD-FFT參數識別

圖4顯示了測點Ac3的加速度功率譜。利用AMD(Analytical Mode Decomposition)[12]法對原始數據進行濾波處理，帶寬分別取0～2 Hz、4～6 Hz、13～15 Hz。利用Bayesian TDD法計算得到懸臂梁的前3階模態振型如圖5所示，從圖中可看出，計算值與理論值吻合很好，表明Bayesian TDD方法能有效地識別結構振型。

圖3 測點Ac3加速度時程

圖4 測點Ac3加速度功率譜

圖5 Bayesian TDD-FFT法識別出的前3階振型

Fig.5 The first three modes identified by Bayesian TDD-FFT method

表1列出了部分模態參數的識別結果，其中，頻率的識別值與理論值0.806 8 Hz、5.056 8 Hz、14.160 6 Hz很接近，阻尼比的識別值與理論值0.005也基本一致，表明Bayesian TDD-FFT方法能有效地識別結構的動力參數。從變異系數(Coefficient of Variation，COV)可以看出，阻尼比的變異系數是最高的，可見在模態識別中阻尼比的值不易準確測定。

4 工程應用

4.1 上海中心加速度響應采集系統

為了采集自然激勵下的加速度數據，上海中心大廈共沿樓層高度方向共安裝了71個Lance加速度傳感器。加速度傳感器的詳細參數,見表2。

表1 Bayesian TDD-FFT模態參數識別結果

Tab.1 Modal parameters identification based on Bayesian TDD-FFT method

階次參數頻率阻尼比激勵譜密度誤差譜密度第一階MVP0.807 40.004 63.823 71.798 5E-2標準差5.368 4E-34.103 2E-41.727 2E-31.490 5E-4COV/%0.664 98.920 04.517 1E-20.828 7第二階MVP5.050 20.004 92.611 31.262 0E-2標準差0.039 44.790 4E-46.072 6E-41.503 4E-4COV/%0.780 29.855 52.325 5E-21.191 3第三階MVP14.167 80.004 82.168 11.940 8E-3標準差0.064 763.843 8E-45.444 2E-43.655 5E-5COV/%0.457 18.007 12.511 0E-21.883 5注:變異系數COV=標準差/MPV

表2 加速度傳感器參數

4.2 良態風作用下模態參數識別

2018年4月12日～14日上海浦東風速達到4～5級，尤其是4月13日的瞬時最大風速達到30 m/s，實時風速見圖6。考慮到識別效果和儀器的采樣精度，取風速較大時間段(13日7時～8時)的加速度數據進行模態參數識別。考慮到篇幅，本文只選取21層、36層、51層、67層、83層、100層、117層、124層的X(東西)向加速度進行識別。圖7是各層加速度的樣本(600 s)和117層加速度功率譜密度。

圖6 2018年4月13日實時風速

利用AMD(Analytical Mode Decomposition)法對原始數據進行濾波處理，帶寬分別取[0，0.17]、[0.17，0.25]。對濾波后的數據應用Bayesian TDD-FFT方法進行處理，計算結果見圖8～14和表3所示。其中，圖8和圖9顯示了前兩階模態和單自由度響應概率分布。前兩階模態與Fast Bayesian FFT求解結果相似。單自由度響應的概率密度分布與Gauss分布十分吻合，驗證了Bayesian TDD假定和推論的正確性。

表3 Bayesian TDD-FFT前兩階模態參數

Tab.3 The first two order modal parameters based on Bayesian TDD-FFT method

(a) 不同樓層600 s加速度樣本

(b) 117層加速度功率譜

Fig.7 Acceleration sample during the 600 s and power spectral density of 117th floor acceleration

X方向前兩階模態參數抽樣樣本如圖10所示，MPV值見表3。由于抽樣起始點是MPV值，因此抽樣樣本收斂在MPV值附近。從圖11和圖12可以看出，頻率和阻尼比的后驗概率密度分布與Gauss分布相吻合，而激勵譜密度和誤差譜密度的后驗概率密度分布離散且與Gauss分布相差較大。

圖13和圖14為前兩階的模態參數后驗聯合分布直方圖，從圖中可看出，頻率與阻尼比沒有相關性，而激勵譜密度和預測誤差譜密度呈近似的線性相關性：預測誤差譜密度與激勵譜密度呈反比例關系，即噪聲譜密度隨著激勵譜密度的增加而降低，這與信號采集時的常規認知相一致。

(a) X方向第一階陣型

(b) X方向第一階模態響應概率分布

Fig.8 The 1stmode shape and the probability distribution of the 1stmode response inXdirection

(a) X方向第二階陣型

(b) X方向第二階模態響應概率分布

Fig.9 The 2ndmode shape and the probability distribution of the 2ndmode response inXdirection

(a) 第一階模態參數抽樣樣本(b) 第二階模態參數抽樣樣本

圖10 X向第一階模態和第二階模態參數抽樣樣本

圖11 第一階模態參數后驗概率密度分布

Fig.11 The posterior probability density distribution of the 1stmodal parameters

(a) 頻率后驗概率分布圖(b) 阻尼比的后驗概率分布圖(c) 激勵譜的后驗概率分布圖(d) 誤差譜后驗概率分布圖

圖12 第二階模態參數后驗概率密度分布

Fig.12 The posterior probability density distribution of the 2ndmodal parameters

表4為應用Fast Bayesian FFT方法計算得到的前兩階模態結果，Fast Bayesian FFT法采用的濾波帶寬與Bayesian TDD-FFT法相同。從表3和表4可看出：Bayesian TDD-FFT法所求得的頻率和阻尼與Fast Bayesian FFT法接近，但前者的變異系數更低，結果更穩定。Bayesian TDD-FFT法所求激勵譜密度大于Fast Bayesian FFT法所求值，而誤差譜密度小于后者，這是由于在Bayesian TDD法處理時舍棄了部分噪聲分量。表中變異系數表明，兩種方法求解出的頻率、激勵譜密度和誤差譜密度的后驗不確定性小，而阻尼比的不確定性大。

表4 Fast Bayesian FFT前兩階模態參數

Tab.4 The first two order modal parameters based on Bayesian FFT method

階次參數頻率阻尼比激勵譜密度誤差譜密度第一階MVP0.1060.0141.193E-066.114E-05標準差6.546E-42.206E-39.278E-81.190E-8COV/%0.6115.927.771.947 2E-02第二階MVP0.2030.0131.226E-071.241E-06標準差3.854E-41.446E-39.526E-95.479E-12COV/%0.1910.987.764.414 9E-04注:變異系數COV=標準差/MPV

表5列出了SSI(Stochastic Subspace Identification，隨機子空間)法和Bayesian TDD-FFT法識別的模態頻率和阻尼比的對比結果，從表中可看出，兩種方法的頻率識別結果基本相同，但阻尼比的識別結果存在一定差異，進一步說明了阻尼比識別的不確定性相對較大。相比SSI方法，Bayesian TDD-FFT法的優勢在于能夠獲得參數的后驗概率分布，同時可以對參數的不確定性進行評估。

表5 Bayesian TDD-FFT和SSI方法參數識別結果對比

Tab.5 Result comparison of parameter identification using SSI and Bayesian TDD-FFT method

階次Bayesian TDD-FFTSSI頻率/Hz阻尼比頻率/Hz阻尼比第一階 0.1060.012 10.1060.014 2第二階 0.2030.015 60.2040.017 6

5 結 論

Bayesian TDD法將多測點多自由度信號轉換成單測點單自由度信號，有效地降低了監測信號處理的復雜程度。Bayesian TDD和蒙特卡羅方法的聯合應用,避免了Bayesian FFT法中病態Hessian矩陣的求逆運算，能夠有效地對結構模態參數進行評估。實際工程應用結果表明：

(a) 頻率與阻尼比后驗聯合分布直方圖

(b) 頻率與阻尼比聯合分布熱力圖

(c) 激勵譜密度與誤差譜密度后驗聯合分布直方圖

(d) 激勵譜密度與誤差譜密度聯合分布熱力圖

Fig.13 Posterior joint distribution histogram and hot map of the 1stmodal parameters

(1) Bayesian TDD-FFT法能夠有效地識別結構的模態參數，利用蒙特卡羅抽樣分析可以獲得及驗證模態參數的后驗概率分布，此為Bayesian TDD-FFT法優于Fast Bayesian FFT法之處；

(2) 對比分析表明：Bayesian TDD-FFT法雖然是一種數值解法，但仍能夠達到Fast Bayesian FFT法的解析解精度；相比SSI法，Bayesian TDD-FFT法的優勢在于能夠獲得參數的后驗概率分布，同時可以對參數的不確定性進行評估；

(a) 頻率與阻尼比后驗聯合分布直方圖

(b) 頻率與阻尼比聯合分布熱力圖

(c) 激勵譜密度與誤差譜密度后驗聯合分布直方圖

(d) 激勵譜密度與誤差譜密度聯合分布熱力圖

Fig.14 Posterior joint distribution histogram and hot map of the 2ndmodal parameters

(3) 上海中心大廈的實測數據分析結果表明：頻率、激勵譜密度和誤差譜密度的后驗不確定性小，而阻尼比的不確定性大；頻率和阻尼沒有相關性，而激勵譜密度和誤差譜密度存在近似的線性相關性，即反比例關系。