快速特征系統實現算法用于環境激勵下的結構模態參數識別

蒲黔輝， 洪 彧， 王高新， 李曉斌

(1.西南交通大學 土木工程學院，成都 610031； 2.東南大學 土木工程學院，南京 210096)

試驗模態參數識別仍是目前研究結構動力學特性的主要方法，識別的結構模態參數可以用于結構有限元模型修正、損傷識別和結構的實時安全監測等。然而土木結構往往都具有結構形式復雜、體型巨大、周圍環境復雜多變等特點，例如高樓和橋梁，如果運用傳統的錘擊或者激振器方式來獲得結構強迫振動往往效果不佳，可能僅有激勵點附近的布設的傳感器才能采集到結構的響應信號，不利于進行整體結構的模態分析。而如果采用大能量的激勵會增加試驗成本，甚至導致結構產生永久性的損傷，所以一般情況下運營中的土木結構不允許采用此種模態試驗方法。此外，運營中的土木結構，特別是橋梁，一般不能長時間完全封鎖路段，而車輛行駛造成結構振動會對模態測試造成誤差。因此，基于環境激勵的模態測試技術成為了土木結構中常用的模態試驗方法。

基于環境激勵的模態測試是指采集結構在地脈動、風荷載、移動車輛荷載等未知荷載作用下產生的振動響應，并提取響應中的自由衰減成分進行模態參數識別。這種方法不需要專門的激勵設備，也不影響結構的正常運營，是一種操作簡單、成本低、安全性高的測試方法。James等[1-2]提出了自然激勵技術(Natural Excitation Technique, NExT)，證明了結構上任意兩個自由度上的位移響應之間的互相關函數滿足結構自由振動方程。Ibrahim等[3-4]提出了ITD(Ibrahim Time Domain)算法，該算法利用自由衰減響應信號進行三次不同延時采樣，構造自由響應數據的增廣矩陣，并由響應與特征值之間的復指數關系，建立特征矩陣，求解特征值和特征向量，再根據模型特征值與振動系統特征值的關系，求解出系統的模態參數。但該方法只適用于單自由度激勵情況。后來，Juang等[5-6]提出了特征系統實現算法(Eigensystem Realization Algorithm, ERA)，該算法是在系統最小實現理論的基礎上發展起來的時域方法，可適用于多輸入多輸出(Multiple Input Multiple Ouput, MIMO)結構系統。該方法需要利用脈沖響應構建Hankel矩陣并對其進行奇異值分解，再利用分解結果構建系統矩陣，從而實現模態參數識別。祁泉泉等[7]提出了擴展特征實現算法(Extended Eigensystem Realization Algorithm, EERA)，把傳統的ERA算法使用范圍擴展至了任何隨機信號。崔定宇等[8-9]提出了一種改進的ERA算法，該方法利用奇異值截斷后重建Hankel矩陣來對信號降噪，使得結構的高階模態更容易識別。包興先等[10]也改進了ERA方法，提出了利用低秩Hankel矩陣逼近方法對響應信號降噪，通過降噪信號后進行模態參數識別。目前采用ERA來進行模態參數識別方法已被國內外的學者廣泛使用[11-12]。

由于土木結構復雜幾何形式，往往需要較多的傳感器采集結構振動信號才能得到完整的模態信息。并且對于環境激勵而言，充足的采樣時間才能保證結構模態分析的可靠性[13]。因此，采用ERA算法構建的Hankel矩陣往往尺寸巨大，對于大型矩陣進行奇異值分解需要消耗大量的運算和儲存時間。針對這一問題，本文提出的快速特征系統實現算法，采用特征值分解來代替奇異值分解，僅利用非零的特征值和對應的特征向量來表達系統矩陣，從而實現模態參數識別。本文采用了一個四層框架仿真模擬來驗證該方法的實用性，并且把該方法運用到了一座人行橋的模態參數識別上。

1 環境激勵下的結構模態參數識別理論

1.1 自然激勵技術

對于一個具有n個自由度的線性結構，在時刻t的振動微分方程為

(1)

假設結構的輸入外力和響應都是平穩的隨機過程，在式(1)左右兩端同時乘上s時刻自由度i上的位移響應qi(s)，并對等式的左右兩端取期望，整理可得

(2)

式中：τ=t-s>0，Rqqi(τ)為結構所有自由度的位移響應與i自由度位移響應的互相關函數。由式(2)可以看出結構位移響應之間的互相關函數是與結構自由振動有類似的表達公式，James等也詳細地證明了利用采集的結構位移響應之間的互相關函數可以進行模態參數識別。同樣的，按照類似的推導方式也可以得出：加速度響應之間的互相關函數也可以進行模態參數識別。

1.2 傳統特征系統實現算法

為了避免處理二階求導問題，式(1)可以在狀態空間里表示為一階微分方程。而試驗模態分析是要利用結構振動響應信號反推其模態參數，但是結構響應是傳感器以一定的采樣頻率間斷采集得到的。所以式(1)需要在狀態空間中表示為離散的線性時不變系統，如下所示

x[k+1]=Adx[k]+Bdp[k]

(3a)

y[k]=Cdx[k]+Ddp[k]

(3b)

式中：Ad∈R2n×2n,Bd∈R2n×np,Cd∈Rny×2n,Dd∈Rny×np分別為離散系統的系統矩陣、控制矩陣、觀測矩陣和輸入觀測矩陣；x∈R2n×1為狀態空間變量；y∈Rny×1為系統輸出，ny為結構中傳感器的個數，系統輸出可以為加速度、速度和位移，取決于觀測矩陣和輸入觀測矩陣的改變；p∈Rnp×1為系統輸入，np為輸入的個數。

假設結構沒有初始位移和速度，即x[0]=0，把式(3a)代入式(3b)，經過多次迭代可得

(4)

定義Hankel矩陣

(5)

其中，

(6)

H[k-1],H[k]∈Rnyr×nps且s,r≥2n；式(6)是將頻響函數矩陣表達式代入式(5)整理后得到。劉宇飛等提出精確的模態分析要求Hankel矩陣尺寸需要足夠大，至少要包括信號振動的幾個周期。

當k=0時，有H[0]=ZW，對H[0]進行奇異值分解

H[0]=U∑VT

(7)

定義H[0]的廣義逆矩陣為H[0]+

H[0]+=V∑-1UT

(8)

為了找到Hankel矩陣與頻響函數(包含結構系統矩陣)之間的關系，定義Ey=[I0 … 0]T∈Rnyr×ny，Ep=[I0 … 0]T∈Rnps×np，它們有以下關系存在

(9)

把式(5)～式(8)代入式(9)的左邊進行計算分解，可以得到Ad，Bd，Cd的表達式

(10)

1.3 快速特征系統實現算法

通常情況下，由于H[0]矩陣的維數較大，對其進行奇異值分解需要占用計算機相當部分內存，因此本文提出了一種快速特征系統實現算法。利用特征值分解(一種更經濟的計算方法)來代替奇異值分解，僅采用非零的特征值與相應的特征向量來構建系統矩陣，相比傳統的ERA算法，既提高了運算速度，又減少了數據儲存。

定義一個新的矩陣S=H[0]TH[0]，把式(7)代入S后整理可得

S=V∑T∑VT

(11)

在式(11)的左部分右乘V，又由左右奇異向量矩陣的正交性可知，VTV=I，UTU=I，經過整理后可得

SV=V( ∑T∑)=VΛ

(12)

式(12)為S矩陣的特征方程，其中Λ為對角矩陣，對角線上的元素為S矩陣的特征值。由式(12)可以看出，S矩陣的特征值為Hankel矩陣H[0]特征值的平方，并且S矩陣的特征向量矩陣等同于H[0]的右奇異向量矩陣。

矩陣H[0]與矩陣H[1]存在以下關系，定義傳遞矩陣T

H[1]=H[0]T

(13)

把式(7)、式(12)和式(13)代入式(10)便可利用S矩陣特征分解的結果重新表達系統矩陣Ad，控制矩陣Bd，觀測矩陣Cd，

(14)

值得注意：式(13)中的傳遞矩陣T不能直接按照廣義求逆的方法得到(H[0]和H[1]都有可能不是滿秩矩陣)。把T矩陣分塊表示

(15)

(16)

T2=VΛ-1VTH[0]T

(17)

得到了離散系統的最小實現矩陣Ad后，對其進行特征值分解

AdΦd=ΦdΔd

(18)

式中：Δd=diag(σd1σd2…σd2n)∈C2n×2n為特征值矩陣；Φd=[Φd1Φd2…Φd2n]∈C2n×2n為特征向量矩陣。

根據離散系統與連續時間系統之間的關系，可以求得結構的復模態參數。

(19)

(20)

式中：“abs”為絕對值；“real”為實部。

1.4 快速特征系統實現算法用于環境激勵下的結構模態參數識別的步驟

步驟1 對采集信號進行濾波，降噪等預處理；

步驟2 選定參考點，并計算其它響應與參考點的互相關函數；

步驟3 利用互相關函數構建Hankel矩陣H[0]與對稱矩陣S；

步驟4 對S矩陣進行特征值分解求得特征值對角陣Λ和特征向量矩陣V；

步驟5 求解T2矩陣，從而求得T矩陣；

步驟6 組建離散系統的系統矩陣Ad，控制矩陣Bd，觀測矩陣Cd，并對Ad進行特征值分解；

步驟7 把離散系統的模態參數轉換為連續系統(結構)的模態參數。

2 數值仿真

采用文獻[14]中的四層框架模型，仿真其在縮小比例的El Centro地震加速度下的各種結構振動，框架模型見圖1。仿真模型每層樓板質量均為4.99 kg；層間剛度k1=1 401.1 N/M、k2=1 576.1 N/M、k3=1 225.9 N/m、k4=1 050.8 N/m，假設每層有黏滯阻尼：c1=1.40 Ns/m、c2=1.31 Ns/m、c3=1.49 Ns/m、c4=1.23 Ns/m，其它結構信息參可參見該文。

圖1 用于仿真的框架模型Fig.1 Simulated frame model

假設在框架每層樓板中心位置均安裝可測得沿圖中地震方向的位移計、速度計和加速度計，采樣時間為40 s，采樣頻率為100 Hz。地震加速度與仿真得到的#4測點的位移、速度、加速度響應如圖2所示。為了避開模態節點，選擇#1測點為參考點，可以得到#1測點和其它點間的互相關函數，現以#1測點與#4測點的互相關函數為例，見圖3。

對于自然激勵試驗而言，信號時長的選用可直接根據互相關函數圖像考慮，選用時長應盡量包含振動衰減的完整過程，在本次仿真算例中選擇時長0～8 s的數據。利用互相關函數組建Hankel矩陣，再按照1.4節中步驟進行模態參數識別，結果見表1。表1中的理論值是按照質量、剛度、阻尼矩陣直接求解得出的模態參數精確值。從表1可以看出，對于該四層結構，FERA算法能采用位移、速度、和加速度進行模態參數識別(利用相同的MATLAB程序，只改變響應的類型)，且均能得到較好結果，可計算出頻率的最大誤差僅為0.5%；阻尼比誤差相對較大，最大誤差約為25%；模態振型誤差幾乎為0，MAC值全為1，MAC值算法見文獻[15]。與傳統的ERA算法相比，FERA算法的精度相當，并均顯示出利用加速度響應進行模態參數識別結果與理論值最為接近。僅以此小結構為例，FERA算法就較ERA算法速度提高了約一半。

圖2 地震激勵和#4測點結構響應Fig.2 Earthquake excitation and structural responses on #4

圖3 #1測點與#4測點的互相關函數Fig.3 Cross correlation functions of responses on #1 and #4

模態1模態2模態3模態4時間/s理論值a0．92262．47793．84404．9461b0．00290．00870．01330．0152c1．00001．00001．00001．0000快速特征實現算法位移a0．92652．47393．83914．9509b0．00230．00860．01020．01260．28601c1．00001．00001．00001．0000速度a0．92792．47363．83964．9509b0．00240．00890．00910．01150．28292c1．00001．00001．00001．0000加速度a0．92332．47413．84094．9508b0．00330．00910．01230．01210．27863c1．00001．00001．00001．0000特征實現算法位移a0．92652．47393．83914．9511b0．00230．00860．01020．01200．38067c1．00001．00001．00001．0000速度a0．92812．47343．84294．9545b0．00260．0090．00970．01160．41137c1．00001．00001．00001．0000加速度a0．92392．47243．84064．9551b0．00310．00980．01110．01200．41140c1．00001．00001．00001．0000注：a為頻率/Hz；b為阻尼比；c為MAC值

3 人行橋環境激勵試驗

為了驗證FERA算法可以用于實際結構的模態參數識別，本文對美國佐治亞理工校園內的一座人行橋進行了模態參數識別，如圖4所示。該橋為簡支梁橋，長30.18 m、寬2.13 m、高2.74 m，其它結構信息、傳感器信息等可參見文獻[16-17]。該橋測試數據來自佐治亞理工智能結構實驗室，為了獲得較為完整的結構模態振型，在橋的幾何節點處都布置了加速度測點(共44個)，采集每個測點的豎向和橫向加速度，傳感器布置參見圖5。為了解決測點數量多而傳感器數量受限制的問題，對該橋進行了測試分段。首先在測點1～14安裝雙向無線加速度傳感器，采集自然激勵下的加速度響應，然后把加速度傳感器移動至測點11～26，再采集加速度響應，然后按照相同方法依次采集測點23～34和31～44的加速度響應。每次環境激勵測試的采樣時間為150 s，采樣頻率為100 Hz。以#8測點為例，加速度響應如圖6所示。

圖4 實驗人行橋Fig.4 Experimental pedestrian bridge

圖5 傳感器布置Fig.5 The allocation of accelerometers

圖6 #8測點橫向加速度Fig.6 Transverse accelerations on #8

測試結束后需要對各個分段結構單獨進行模態參數識別，然后利用平均的方法得出整體橋梁的頻率和阻尼比，并且利用共節點拼接分段模態振型來得到結構的完整模態振型。以第一分段為例，選擇#8測點橫向自由度為參考點，計算其他所有自由度上加速度響應與該自由度加速度響應的互相關函數，示例如圖7所示。利用互相關函數組建Hankel矩陣，再按照1.4節中步驟對人行橋進行模態參數識別。最終識別出前4階振型，結果如圖8所示，第1階振型為橫向模態、第2階振型為豎向模態、第3階振型為橫向模態、第4階振型為扭轉模態。為了與模態參數識別結果對比，利用ANSYS對人行橋進行了仿真模態分析，結果顯示試驗模態與理論模態的頻率值有較小差距，但是二者的

模態振型能很好對應，說明FERA算法同樣適用于真實結構的模態參數識別(見圖9)。本研究同樣使用NExT+ERA方法對該橋進行了模態參數識別，與NExT+FERA方法的識別結果很接近，所以不再贅述。但總體上看，本文提出的FERA算法較ERA方法快2～3倍。

圖7 #8測點橫向與豎向加速度的互相關函數Fig.7 Cross correlation functions of the transverseand vertical accelerations on #8

圖8 模態參數識別結果Fig.8 Modal identification results

圖9 ANSYS模態參數識別結果Fig.9 Modal identification results from ANSYS analysis

4 結 論

本文提出的FERA模態參數識別算法是對傳統ERA算法的一種改進，利用特征值分解代替了復雜的奇異值分解。本文通過四層框架的數值仿真試驗和人行橋的自然激勵試驗驗證了FERA算法的準確性和優勢。有如下結論：

(1)四層框架數值模擬表明FERA算法較ERA算法有著相同的精度，適用于各種動力響應。但較其他兩種響應而言，利用加速度響應進行模態參數分析能得到更精確的結果，所以推薦使用加速度響應進行模態參數識別。

(2)人行橋的自然激勵試驗表明了FERA算法同樣適用于真實結構，并且速度統計顯示，結構越大，FERA算法較ERA算法的速度優勢更明顯。

(3)FERA算法可以提高結構安全監測、結構模型修正、損傷識別的效率。

[ 1 ] JAMES G H, CARNE T G, LAUFFER J P. The natural excitation technique (NExT) for modal parameter extraction from operating wind turbines: SAND-92-1666[R]. 1993.

[ 2 ] JAMES G H, CARNE T G, LAUFFER J P, et al. Modal testing using natural excitation[C]∥Proceedings of the International Modal Analysis Conference. Albuquerque: SEM Society for Experimental Mechanics Inc., 1992: 1208.

[ 3 ] IBRAHIM S R, MIKULCIK E C. A method for the direct identification of vibration parameters from the free response[J]. The Shock and Vibration Bulletin, 1977, 47(4): 183-

198.

[ 4 ] 楊佑發，李帥，李海龍. 環境激勵下結構模態參數識別的改進ITD法[J]. 振動與沖擊, 2014, 33(1): 194-199.

YANG Youfa, LI Shuai, LI Hailong. Improved ITD method for structural modal parameter identification under ambient excitation [J]. Journal of Vibration and Shock, 2014, 33(1): 194-199.

[ 5 ] JUANG J N, PAPPA R S. An eigensystem realization algorithm for modal parameter identification and model reduction [J]. Journal of Guidance, Control, and Dynamics, 1985, 8(5): 620-627.

[ 6 ] JUANG J N, PAPPA R S. Effects of noise on modal parameters identified by the eigensystem realization algorithm [J]. Journal of Guidance Control and Dynamics, 1986, 9(3): 294-303.

[ 7 ] 祁泉泉，辛克貴，崔定宇.擴展特征系統實現算法在結構模態參數識別中的應用[J]. 工程力學, 2011, 28(3): 29-

34.

QI Quanquan, XIN Kegui, CUI Dingyu. The application of the extended eigensystem realization algorithm for structural modal parameter identification [J]. Engineering Mechanics, 2011, 28(3): 29-34.

[ 8 ] 崔定宇，辛克貴，祁泉泉.擴展特征系統實現算法的模態參數識別特性研究 [J]. 工程力學, 2013, 30(8): 49-53.

CUI Dingyu,XIN Kegui,QI Quanquan. Modal parameter identification research for the extended eigensystem realization algorithm [J]. Engineering Mechanics, 2013, 30(8): 49-

53.

[ 9 ] 辛峻峰, 張永波, 盛進路. 一種改進的特征系統實現算法及其在海洋平臺上的應用[J]. 振動與沖擊, 2015, 34(24): 197-201.

XIN Junfeng, ZHANG Yongbo, SHENG Jinlu. Improved eigen system realization algorithm and its application in the analysis of ocean platform [J]. Journal of Vibration and Shock, 2015, 34(24): 197-201.

[10] 包興先，李昌良，劉志慧. 基于低秩Hankel矩陣逼近的模態參數識別方法 [J]. 振動與沖擊,2014, 33(20): 57-

62.

BAO Xingxian, LI Changliang, LIU Zhihui. Modal parameters identification based on low rank approximation of a Hankel matrix [J]. Journal of Vibration and Shock, 2014, 33(20): 57-62.

[11] REYNDERS E. System identification methods for (operational) modal analysis: review and comparison[J]. Archives of Computational Methods in Engineering, 2012, 19(1): 51-124.

[12] 劉宇飛, 辛克貴, 樊健生, 等. 環境激勵下結構模態參數識別方法綜述[J]. 工程力學, 2014, 31(4):46-53.

LIU Yufei , XIN Kegui, FAN Jiansheng, et al. A review of structure modal identification methods through ambient excitation [J]. Engineering Mechanics, 2014, 31(4): 46-

53.

[13] CAICEDO J M. Practical guidelines for the natural excitation technique (NExT) and the eigensystem realization algorithm (ERA) for modal identification using ambient vibration[J]. Experimental Techniques, 2011, 35(4): 52-58.

[14] HONG Y, LIU X, DONG X, et al. Experimental model updating using frequency response functions[C]∥Proceedings of SPIE: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems. Las Vegas: SPIE， 2016.

[15] ALLEMANG R J, BROWN D L. A correlation coefficient for modal vector analysis[C]∥ Proceedings of the 1st International Modal Analysis Conference. Orlando: IMAC， 1982： 110-116.

[16] ZHU D, GUO J, WANG Y, et al. Field validation of flexure-based mobile sensing nodes on a space frame bridge[C]∥ Proceedings of the 8th International Workshop on Structural Health Monitoring.Stanford: IWSHM, 2011: 13-15.

[17] ZHU D, CHO C, WANG Y. Finite element model updating of a space frame bridge with mobile sensing data[C]∥ Proceedings of SPIE, Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems. San Diego: SPIE， 2012.