基于Welch 法的協方差隨機子空間方法的模態參數識別1)

李雪艷 官宇航 羅銘濤 吳博宇

(暨南大學力學與建筑工程學院,重大工程災害與控制教育部重點實驗室,廣州 510632)

引言

環境激勵下的模態參數分析[1-3]是在風荷載、車輛荷載、地震和沖擊荷載的作用下,對僅有輸出(響應)的工程結構進行模態參數識別的一類結構動力學反問題.相較于傳統的模態參數識別方法,這種方法無需測定未知的激勵,不影響結構的正常運行,不僅費用低,而且可以識別出結構正常運行時的模態,因而該類方法具有廣闊的應用前景[4].

隨機子空間法(stochastic subspace method,SSI)[5-8]作為一種適用于環境激勵下的模態參數識別的時域識別方法,基于系統的離散時間狀態空間模型,并假定環境激勵為平穩白噪聲,再經過矩陣QR 分解或者奇異值分解之后得到系統狀態矩陣,再利用特征值和特征向量計算出頻率、振型和阻尼比等模態參數[9].由于系統的真實階次難以確定,Peeters等[10]提出穩定圖法來確定系統階次,以此來識別出更加真實準確的模態.隨機子空間法可以分為協方差驅動隨機子空間法[11-12]和數據驅動隨機子空間法[13-14],前者需要由Hankel 矩陣計算輸出響應的協方差矩陣,再經過奇異值分解得到系統矩陣,后者由Hankel 矩陣得到將來行空間投影到過去行空間的投影矩陣,再進行卡爾曼濾波和QR 分解來求解出系統矩陣,避免了協方差的計算.然而,實際工程結構受到的環境激勵是較為復雜的,因此把環境激勵假定為平穩白噪聲可能會使得識別結果中出現虛假模態或者真實頻率出現誤差[15].

為了改進隨機子空間法,目前已有許多嘗試.秦世強等[16]利用大橋頻率低、模態密集等特點,提出引入屏蔽信號限制經驗模態分解(EMD)過程中的帶寬以消除模態混疊.Lu 等[17]提出了一種能實現非同步數據的模態參數識別隨機子空間法,該方法的核心是使用模態的平均相位偏移(mean phase deviation,MPD)來找出實際的時滯以及同時計算實際振型分量.Lin 等[18]把協方差驅動隨機子空間法(COV-SSI)和數據驅動隨機子空間法(DATA-SSI)相結合,通過對改進投影矩陣進行奇異值分解,避免使用原始大維數據矩陣,來求解系統矩陣,可以有效提高隨機子空間法的計算效率.陳永高等[19]嘗試采用滑動窗口和相似度來確定系統的真實階次,引入改進集成經驗模式分解算法消除噪聲影響,剔除虛假模態并避免模態遺漏.周筱航等[20]使用滑窗技術改進隨機子空間法.Jin 等[21]提出一種短時隨機子空間識別(ST-SSI)框架,通過分析駛過橋梁的車輛的動態響應來識別具有時變性的車橋耦合系統的頻率.張永祥等[22]構建了一種基于譜系聚類的隨機子空間模態參數自動識別算法,其基本原理是將模式樣本按距離準則逐步聚類.Wen 等[23]提出了一種基于滑動窗口模糊C 均值聚類算法和確定性隨機子空間辨識(SC-CDSI)相結合的時域識別方法,實現非線性時變結構模態參數的在線跟蹤和識別.首先,對非線性時變結構的輸入輸出信號進行加窗分割;其次,為了滿足穩定圖中有效模態的自動識別,以模糊C 均值聚類算法(fuzzy C-means clustering analysis,FCM) 為基礎,以頻率、阻尼比和振型為聚類元素,實現了有效模態的自動智能選擇.

盡管通過信號預處理、加窗函數、奇異值截斷法和聚類分析等技術的引入,能部分改善隨機子空間法的效果,但是隨機子空間法在實際應用中仍然存在模態遺漏、虛假模態、計算量大和自動定階困難等問題.實際工程結構承受著多樣的環境激勵,隨機子空間法要準確地識別出模態參數,則需要較長時間段的輸入數據,但這樣會導致計算效率低下.所以本論文提出使用Welch 方法[24]和Hanning 窗,對信號在頻域進行去除激勵和噪聲影響的處理,然后生成包含結構更多固有模態和較少激勵頻率的Toeplitz矩陣,再進行特征值分析,并使用模糊C 均值聚類算法[25-26]和模態平均相位偏移(MPD)[27-28]指標來識別結構的模態參數.并把本文所提出的方法應用于一座大跨懸索橋的實測加速度響應和一座高層建筑結構的實測加速度響應數據分析和模態參數識別,演示了所提方法的可行性和有效性.

1 基于Welch 法的協方差驅動隨機子空間方法

1.1 識別狀態矩陣

假設yk=是結構的m個測點的tk時刻的加速度響應,由t0到t2i+j-2時刻的加速度響應構造如下Hankel 矩陣

把矩陣分為兩部分,“過去”和“將來”輸出數據分別為Yp∈Ri×j和Yf∈Ri×j,由Hankel 矩陣可以得到Toeplitz 矩陣

其中

式中,C為跟測量位置有關的矩陣,A為離散狀態矩陣,包含了結構的狀態信息,如頻率、振型和阻尼等參數,G為狀態-輸出協方差矩陣,i和j為數據長度,由于j為有限數,所以導致式(3)的左邊不能嚴格等于右邊,使得結構固有頻率不能從眾多激勵和噪聲頻率中突顯出來,為后續的子空間識別帶來虛假頻率.

而結構健康監測系統測試的加速度響應信號具有連續不間斷和數據海量的特點,另一方面環境不確定因素多和激勵復雜,激勵數據未知.為了消除這些不利因素,可把時域的海量加速度響應信號通過加窗函數和傅里葉變換到頻域進行平均計算.

式中,ω 是圓頻率變量,j*是復數符號,Mr1,Cr1,Kr1分別為第r1階模態質量,模態阻尼和模態剛度,F(ω)是激勵力向量的傅里葉變換,Φl1r1是第r1階振型Φr1的第l1個分量,N為結構的總自由度數.實際運營結構,如橋梁承受著多種激勵,例如車流、風和地脈動等,但是這些外部激勵具有隨機性和不確定性,所以F(ω)具有隨機性,而結構的固有頻率不受外部環境影響,對事先給定時間長度的加速度響應信號的l段功率譜進行平均,可以降低激勵頻率的功率,而結構固有頻率的功率將更突出,成為功率譜曲線上的少數極值點.再對進行逆傅里葉變換得到互相關函數

并把式(7)代入式(2)中則可以提高Toeplitz 矩陣的精度.

把式(3)代進Toeplitz 矩陣,則有

其中,Oi為可觀矩陣,Γi為可控反轉矩陣,再對由式(7)的相關函數構造的Toeplitz 矩陣進行奇異值分解

其中,U和V是正交矩陣,S是奇異值對角矩陣,S1和S2分別為非零奇異值和零奇異值對角矩陣,則有

式中rs為非零奇異值的個數.

離散狀態矩陣A可由下式求得

其中 (*)+表示廣義逆.

1.2 頻率識別

由離散狀態矩陣A進行模態參數識別,對A進行特征值分解可得

1.3 模糊C 均值聚類(FCM)分析

其中,α 是加權指數,J(·) 為目標函數,uij和dij分別為xj對pi的隸屬度和歐式距離.

引入 ρ=[ρ1,ρ2,···,]算子則有

對所有輸入參量求偏導可知式(20)取最小值的必要條件為

由離散狀態矩陣A可以得到一組頻率 [ω1,ω2,···,],通過變換式(10)中的r∈[rmin,rmax],可以得到多組[ω1,ω2,···,]r,把所有頻率作為樣本進行FCM 聚類分析.給定加權指數 α=2[29],選擇最多的一組頻率 [ω1,ω2,···,ωc] 作為初始聚類中心P,初始隸屬度矩陣U按下式計算

1.4 平均相位偏移(MPD)定階

當結構具有低阻尼和對稱正定的剛度矩陣時,結構振型為實數振型,由式(18)得到的各階振型的各個分量的相位角應當相同或者相近,所以其對相位平均值的偏移較少.或者說相應模態振型的相位分散程度較少.

設第i階模態的平均相位為[30]

其中,φir為第r階振型 Φr的第i個分量的相角,m為每一階振型分量總數,也即測點數,Wi為每一階振型分量的加權系數,可取1 或者 |Φir|,本文取后者即第r階振型的第i分量的模.

第r階振型第i分量的相角 φir按下式計算

式中,Re(·) 和 Im(·) 分別為變量的實部和虛部,若arctan(Re(Φir)/Im(Φir)) ≥0,φir原值不變,為負時后者加 π/2 .

各階振型向量的平均相位偏移為

計算結果為偏移的弧度值,對于實模態,換算成角度值之后的MPDr計算結果應趨近于0,一般倘若MPD≥20°則認為相位偏移較大,對應的模態為虛假模態[30].

最后通過聚類的樣本數和平均相位偏移來共同確定模態階數.基于Welch 法,FCM 和MPD 分析的隨機子空間法在模態參數識別時的流程圖如圖1 所示.

圖1 基于Welch 法的隨機子空間法的流程圖Fig.1 Flow of stochastic subspace identification based on Welch method

2 實際工程結構模態參數識別

2.1 實橋模態分析

對位于廣東省的一座大跨懸索橋的實測數據進行分析和演示本文所提出的方法.該橋的監測系統包括SLJ-100 力平衡加速度計、MSM-48 多通道數據采集器、信號電纜、GL-3 B 直流電以及用于控制和采集測量數據的計算機.加速度計的測量范圍為 ± 2g,頻率范圍為0.0～80.0 Hz,靈敏度為1.25(V)/g,噪聲小于0.5g.

以200 Hz 的采樣頻率采集加速度響應.傳感器編號系統及其布置如圖2 所示.包括3 個縱向(x方向) 加速度計(傳感器8,32,20)、5 個橫向(y方向)加速度計(傳感器7,10,13,16,19)和8 個豎向(z方向)加速度計(傳感器9,11,12,14,17,18,21).2020 年3 月2 日早上9 時錄得的5 s 時間長度的加速度響應,如圖3 所示.可以看到豎向加速度響應幅值明顯大于縱向和橫向的,說明主要是車輛載荷激振.

圖2 懸索橋上加速度計的布置Fig.2 Arrangement of accelerometers on suspension bridge

圖3 懸索橋的加速度響應Fig.3 Acceleration response of suspension bridge

對橋板上13 個傳感器(傳感器8,32 和20;7,10,13 和19;9,11 和12,14 和15,17)的加速度響應數據進行分析和模態參數識別.首先對30 s 共6000 個數據點的加速度響應進行協方差驅動的隨機子空間分析,得到的穩定圖如圖4(a)所示.對1 小時的加速度響應數據進行Welch 法加頻域分解法(frequency-domain decomposition,FDD)[31],得到的功率譜曲線如圖4(a)中所示的藍線所示,文中后面圖中的該曲線都為FDD 法計算得到.可以看出由傳統SSI 法得到的穩定圖中的優勢頻率較少,這主要是由于6000 個數據點長度的環境激勵加速度響應數據在生成Hankel 矩陣和Toeplitz 矩陣時,無法有效消除環境激勵等虛假頻率.但是直接使用更長時間的加速度響應,會使得隨機子空間法耗時更長,效率降低.

把1 小時長的各通道的加速度響應數據進行自相關分析,得到3000 數據點長度的相關函數數據代替加速度響應數據來生成Hankel 矩陣,并進而計算Toeplitz 矩陣和進行特征值分析,得到各個狀態矩陣的頻率,進行FCM 聚類分析,并使用MPD 進行模態參數的篩選確定,得到結果如圖4(b)所示.可以看到已經能識別出13 個結構固有頻率.之所以比圖4(a)所示情況得以改善,是因為1 小時長度的自相關分析,顯著降低了噪聲和環境激勵的影響,使得結構固有頻率在相關函數中突顯了出來,但是跟頻域分解法的功率譜曲線相比,有模態遺漏現象.即第1,2,3,4,8 個峰值對應的頻率都未能識別出來.

對1 小時長度的加速度響應數據,進行如式(4)～式(6)所給出的加窗和平均功率譜計算,對功率譜進行如式(7) 所給出的逆傅里葉變換,并直接構建Toeplitz 矩陣,然后進行狀態矩陣計算和特征值分析,再進行FCM 聚類分析,聚類分析以后的穩定圖如圖4(c)中的紅色部分所示,此時得到的類數目為241,遠大于實際模態階數,所以需要去除虛假模態.對每個類所對應的振型進行平均相位偏移(MPD)分析,即通過隨機子空間法得到的每個復數振型的各個分量要具有相同或者相近的相位角,以使得實際振型為實數振型.對聚類分析以后各個類的振型進行MPD分析,得到如圖5 所示的MPD分布圖,對MPD＞ 20°所對應的類進行剔除,留下的模態對應的頻率穩定圖如圖4(c)的藍色豎線部分所示.可以看出第1,2,3 個峰值對應的頻率都能識別出來,說明了基于Welch 法的隨機子空間法在改善模態遺漏方面有顯著提高.

圖4 頻域分解法FDD、傳統SSI、基于相關分析的SSI 和基于Welch 法的SSI 的某大跨懸索橋的頻率識別Fig.4 Frequency identification of a long-span suspension bridge by FDD,traditional SSI,SSI based on correlation analysis and SSI based on Welch method

圖5 對各聚類所對應的振型進行MPD分析得到的MPD分布圖Fig.5 the distribution of MPDobtained by MPDanalysis of vibration modes corresponding to each cluster

在識別出結構的頻率以后,再提取結構的振型,由于該實橋橋板上分三個方向僅安裝了16 個傳感器,且只有13 個傳感器在該時段測得有效數據,所以只能識別得到部分振型分量.在圖6 中顯示了三個頻率各自對應的振型分量,可以看到在相同截面和相同方向的振型分量有相近的值,符合理論結果.也跟用頻域分解法識別得到的對應振型分量進行了比較,發現具有較好的匹配性,再一次證明了所提方法的正確性.

圖6 使用基于Welch 法的SSI 和FDD 法得到的懸索橋的三階振型分量對比Fig.6 Comparison of three orders mode shape components of suspension bridges obtained by using SSI based on Welch method and FDD

使用Welch 法+頻域分解的峰值法、傳統協方差隨機子空間法、基于相關分析的隨機子空間法和基于Welch 法的隨機子空間法的識別結果如表1所示,可以看到基于Welch 法的隨機子空間法,在FCM 和MPD 的協助下,在模態遺漏、虛假模態剔除和自動識別方面具有明顯優勢.

表1 大跨懸索橋的頻率識別結果(Hz)Table 1 The identified frequency results of long-span suspension bridge (Hz)

2.2 環境激勵下的高層建筑模態分析

第二個研究對象是如圖7 所示的七十層高層建筑,在該建筑物的地下和第68 層安裝了三方向加速度傳感器,采用INV3062V 采集儀和91B 加速度傳感器,采樣頻率為100 Hz.分析第68 層的3 個方向的加速度響應,該加速度響應如圖8 所示.

圖7 某高層建筑的傳感器布置圖Fig.7 Sensor layout of a high-rise building

圖8 某大廈的第68 層處測點的加速度響應Fig.8 Acceleration response of a certain measuring point in a high-rise building

首先使用20 s 共2000 個數據點長度的加速度響應,直接進行協方差驅動的隨機子空間分析,如圖9(a)所示,紅色部分數據是傳統SSI 分析的穩定圖結果,藍色部分數據是使用FCM 聚類分析和MPD篩選以后所確定的模態階數,所計算結果也跟頻域分解法的功率譜曲線進行了比較,可以看到模態遺漏問題非常嚴重,這是因為結構固有模態的振動能量跟環境激勵頻率的振動能量相比較小,在隨機子空間法中不能成為主要特征值被識別出來.

圖9 某大廈的頻率識別Fig.9 Frequency identification of a building

對1 小時的加速度響應進行自相關分析,取前3000 個數據點來形成Hankel 矩陣和Toeplitz 矩陣,最后的識別結果如圖9(b)所示,同樣藍色部分為經過FCM 和MPD 分析后的模態階數,相比于圖9(a)的結果,模態階數顯著增加了,但是相對于頻域分解法,仍然有模態遺漏的問題.

對1 小時的加速度響應使用本文中提出的Welch 法分析和直接形成Toeplitz 矩陣,最后的識別結果如圖9(c)所示,可以看到在通過FCM 和MPD分析以后,模態階數顯著增加,跟頻域分解法能識別出來的峰值基本對應,但是能避免頻域分解法中的虛假峰值,從而識別出更準確的結構模態.

頻域分解法、傳統協方差SSI、基于相關分析的SSI 和基于Welch 法的SSI 識別的某高層大廈的模態頻率如表2 所示,從表中看出基于Welch 法的SSI 法在避免模態遺漏,剔除虛假頻率,自動識別方面具有顯著優勢.雖然頻域分解法也能獲得比較全面的結構模態,但是由于所有頻率都會以峰值的形式呈現出來,所以虛假峰值問題也難以避免.

表2 某大廈的頻率識別結果(Hz)Table 2 The identified frequency results of a building (Hz)

3 結論

由于環境激勵復雜,并不完全符合平穩白噪聲的特征,使得傳統的隨機子空間方法在計算效率、剔除虛假模態和避免模態遺漏方面存在不可忽視的問題,限制了該方法在實際工程結構的自動模態參數識別追蹤方面的應用.因此嘗試在頻域對測試數據的模態頻率振動能量進行加強,而對環境激勵頻率和噪聲頻率的振動能量進行降低,所以提出了基于Welch 方法的隨機子空間法.由Welch 方法得到經過頻域平均的振動響應協方差參數,并形成Topelitz 矩陣,然后進行奇異值分析,由不同階數的奇異值構造的狀態矩陣進行特征值分析,得到頻率數據樣本,然后進行模糊C 均值聚類分析,和模態的平均相位偏移分析,剔除MPD 值過大的虛假模態,對SSI 識別的模態進行定階.通過對一座大跨懸索橋的橋板上的13 個加速度傳感器的1 小時數據的分析,和一座高層大廈的3 個加速度傳感器的1 小時數據分析,發現基于Welch 方法的SSI 相比于傳統SSI 和基于相關分析的SSI,在避免模態遺漏和計算效率方面有顯著提高,而相對于頻域分解法則在自動識別和剔除虛假模態方面有明顯優勢.所以本文提出的基于Welch 方法的SSI,在環境激勵下工程結構的模態參數識別追蹤方面具有較好的應用前景.