基于奇異值分解的ERA改進算法及模態(tài)定階

朱 銳， 杭曉晨， 姜 東,3， 費慶國， 靳文冰

(1.東南大學空天機械動力學研究所 南京，211189) (2.東南大學機械工程學院 南京，210096)(3.南京林業(yè)大學機械電子工程學院 南京，210037) (4.國家汽車質(zhì)量監(jiān)督檢驗中心 襄陽，441001)

引 言

結構模態(tài)參數(shù)識別是動力學的核心內(nèi)容之一，可以為動響應預示、故障診斷、安全性評估和結構優(yōu)化設計等提供參考。經(jīng)過幾十年的發(fā)展，模態(tài)參數(shù)識別理論體系已經(jīng)較為成熟，廣泛應用于航空、航天、土木和機械等工程領域。模態(tài)參數(shù)識別與定階問題在理論分析和工程應用方面有重要意義。

特征系統(tǒng)實現(xiàn)算法(eigensystem realization algorithm, 簡稱ERA)是多輸入多輸出的時域模態(tài)參數(shù)辨識方法。該方法最初是由Juang等[1]于1984年提出，核心是利用線性系統(tǒng)中Ho-Kalman的最小實現(xiàn)理論來識別系統(tǒng)模態(tài)參數(shù)。該方法便于確定模態(tài)階次，識別速度較快，同時對具有低頻、模態(tài)密集特征的結構具有較強的辨識能力。Bazan等[2]研究了特征系統(tǒng)算法中噪聲因素對系統(tǒng)極點影響。Zhang等[3]提出了自動特征系統(tǒng)實現(xiàn)算法，并將其應用在朝天門大橋模態(tài)試驗中。Velazquez等[4]根據(jù)線性參數(shù)時變技術改進了特性實現(xiàn)算法，將環(huán)境振動信號作為算法輸入，成功辨識了旋轉系統(tǒng)的模態(tài)參數(shù)。同時，基于奇異值分解[5]的模型定階方法也是近年來國內(nèi)外學者關心的問題，一般的定階方法是從通過奇異值曲線突變特點的角度來考慮。易偉健等[6]提出了殘差期望值法和殘差期望比法，假設噪聲為白噪聲，根據(jù)白噪聲的均值為零的特點確定階次，但是該方法需多次試算來確定Hankel矩陣維數(shù)，計算量較大。周幫友等[7]提出了奇異值插值法，但未深入探討矩陣維數(shù)的確定，同時結果不夠直觀。趙學智等[8]提出了使用奇異值差分譜的概率來反映信號奇異值曲線的轉折點現(xiàn)象。王樹青等[9]提出了利用奇異值的相對變換率的概念來確定系統(tǒng)的階次，該模態(tài)定階方法易受到噪聲干擾。

筆者研究基于奇異值分解的改進的特征系統(tǒng)實現(xiàn)算法(improved eigensystem realization algorithm, 簡稱IERA)，利用結構脈沖響應信號構造初始的Hankel矩陣，根據(jù)Cadzow算法[10]重構新的Hankel矩陣。同時針對奇異值歸一化指標難以確定模態(tài)階次的現(xiàn)象，采用奇異值百分比(singular value percentage，簡稱SVP)反映選定階次增加對系統(tǒng)貢獻量，并給出了模態(tài)的定價指標(RSVP)。通過仿真算例，研究噪聲因素對辨識的模態(tài)參數(shù)及模態(tài)階次確定的影響，最后將所提方法應用在某三廂車排氣系統(tǒng)模態(tài)試驗中。

1 理論基礎

1.1 特性系統(tǒng)實現(xiàn)算法

對于有限、離散時間的線性定常系統(tǒng),其狀態(tài)方程可以表示為

其中:x∈Rn，u∈Rm，y∈Rp分別為狀態(tài)向量、輸入向量和輸出向量;A，B，C分別為系統(tǒng)的狀態(tài)矩陣、輸入矩陣和輸出矩陣。

利用脈沖響應函數(shù)構造Hankel矩陣，即

Hrs(k-1)=

(3)

其中:Y(k) ∈RL×P為脈沖響應函數(shù)矩陣。

Y(k)=

(4)

其中：hij(k)為k時刻在激勵點j和響應點i之間的脈沖響應函數(shù)。

由式(1)和(4)可得Y(k)表達式為

Y(k)=CAk-1B

(5)

將式(3)可改寫成

Hrs(k)=VrAkWs

(6)

(7)

(8)

其中：Vr∈Rrp×n為可觀測矩陣；Ws∈Rn×ms為可控矩陣。

對矩陣Hrs(0)進行奇異值分解

Hrs(0)=PDQT

(9)

其中：P∈Rrp×n；D∈Rn×n；Q∈Rms×n。

由奇異值分解確定的階次獲得系統(tǒng)的最小實現(xiàn)

對矩陣A進行特征值分解，求得特征值矩陣，然后求出特征向量矩陣

Ψ-1AΨ=Z

(12)

其中：Z為特征值矩陣；Ψ為特征向量矩陣

(13)

根據(jù)矩陣A的特征值與系統(tǒng)特征值ζr之間的關系確定模態(tài)頻率ωr和模態(tài)阻尼ζr分別為

1.2 改進后算法

傳統(tǒng)特征系統(tǒng)實現(xiàn)算法在模態(tài)分析過程中，基于奇異值分解后奇異值曲線的突變特點來確定階次，但是有時候突變并不是很明顯，特別在實際工程應用過程中，當噪聲污染嚴重或者模態(tài)階次過高時，大量的噪聲模態(tài)和計算模態(tài)會出現(xiàn)在穩(wěn)定圖[11]上，易出現(xiàn)虛假模態(tài)。為剔除原始數(shù)據(jù)中的噪聲[12]，在特征系統(tǒng)實現(xiàn)算法的基礎上研究一種根據(jù)基于奇異值原理和Cadzow算法改進的特征系統(tǒng)實現(xiàn)算法。Hankel矩陣H由脈沖響應函數(shù)構造，包含無噪聲的數(shù)據(jù)信號和噪聲信號，用奇異值分解進行去噪。

(17)

(18)

(19)

分析流程如圖1所示。通過流程圖可以看出，改進后的算法較傳統(tǒng)的ERA多進行了一次奇異值分解并進行Hankel矩陣重構。為了保證第1次奇異值分解模態(tài)的定階包含真實的階次,l1可適當取值稍大同時滿足l1≥l2。

圖1 改進ERA流程圖Fig.1 Flowchart of improved ERA

1.3 模態(tài)定階：奇異值百分比

利用特征系統(tǒng)實現(xiàn)算法進行模態(tài)參數(shù)識別，需要確定模態(tài)階次[14]。如果選用的階次低于系統(tǒng)真實的模態(tài)階次，就有可能遺漏真實模態(tài)；如果選用的階次高于系統(tǒng)真實的模態(tài)階次，則識別結果中會呈現(xiàn)較多的虛假模態(tài)，對真實模態(tài)造成干擾。針對傳統(tǒng)奇異值曲線可能出現(xiàn)的突變不明顯的特點，以及隨著階次的增加曲線趨向水平漸近線的現(xiàn)象，筆者提出一種奇異值指標——奇異值百分比，能有效地反映所取模態(tài)階次在系統(tǒng)模態(tài)參數(shù)中的比重，進而有效確定模態(tài)的階次

(20)

從奇異值所占百分比的角度考慮，R為計算矩陣維數(shù)，j為選定的階次，通過百分比反映選定階次增加對系統(tǒng)貢獻量，該指標能有效反應選定階次在結構中比重。根據(jù)奇異值百分比曲線的特點，先確定一個最小預估值，再利用奇異值百分比相鄰增量的比值變化的特點作為模態(tài)的定階指標，公式為

(21)

其中:ΔPj+1,j表示為Pj+1-Pj。

因為奇異值百分比反映所選階次在模態(tài)參數(shù)中的比重，相鄰增量比值會在真實模型階次附近出現(xiàn)巨大波動。這是因為真實模態(tài)以后階次是噪聲成分，其增量在同一量級上，模型的定階指標(RSVP)將趨于平緩，接近為1，而在真實模態(tài)處由于增量的量級不同，會出現(xiàn)一個最小值，在曲線將出現(xiàn)最低處即為模態(tài)階次。其優(yōu)點是針對受噪聲影響奇異值突變不明顯、不易確定階次的現(xiàn)象，通過定義的RSVP，可以有效反應模態(tài)階次值，提高模態(tài)的定階精度。

2 算例研究

仿真模型采用六自由度彈簧-振子系統(tǒng)，驗證改進方法在不同噪聲工況下的識別能力，同時也驗證了在不同噪聲工況下采用RSVP指標來確定模態(tài)階次的適用性,并將改進方法應用在某三廂車排氣系統(tǒng)模態(tài)試驗中。

2.1 六自由度彈簧-振子模型

為了驗證改進型的算法識別效果的優(yōu)勢及RSVP指標的模態(tài)定階適用性，建立如圖2所示的六自由度系統(tǒng)。m1=m6=20 kg,mi=10 kg(i=2,3,4，5)，kj=8×106kN/m,(j=1,2，…，7)。通過模態(tài)疊加法確定六階頻率，ci(i=1,2，…，7)的取值使第1階模態(tài)阻尼比為3.00%，其余模態(tài)阻尼比為1.00%，如表1所示。

圖2 六自由度彈簧-振子系統(tǒng)Fig.2 Six-DOF spring-mass system

表1 系統(tǒng)模態(tài)參數(shù)Tab.1 System modal parameters

圖3 脈沖響應函數(shù)Fig.3 Impulse response function

圖4 5%白噪聲Fig.4 5% white noise simulation

圖5 六自由系統(tǒng)穩(wěn)定圖比較Fig.5 Comparison of two stability diagrams of six-DoF system

在質(zhì)量塊m4上施加單位脈沖激勵，在mi(i=1,2,…，6)上采集振動位移響應。模型的分析頻率f=300 Hz，采樣數(shù)N=6 000。以m5點輸出響應為例，加入5%高斯白噪聲后信號曲線與初始曲線不易區(qū)分，圖3,4為無噪聲信號脈沖響應函數(shù)(impulse response function，簡稱IFR)曲線及5%的高斯白噪聲曲線。采用傳統(tǒng)ERA及改進后算法對含有5%噪聲FRF信號進行模態(tài)參數(shù)識別。圖5為兩種辨識算法的穩(wěn)定圖。曲線為5%噪聲下FRF點表示取不同階次時識別的頻率值，點的顏色反映辨識的阻尼比與理論值的誤差。分析可知，改進后算法比傳統(tǒng)ERA的穩(wěn)定圖更加簡潔，能有效剔除虛假模態(tài)，可直觀看出改進后算法的辨識阻尼精度更高。

實心點表示歸一化奇異值；空心點表示奇異值百分比圖6 不同噪聲情況下彈簧-振子系統(tǒng)奇異值曲線及百分比Fig.6 Normalized singular values and SVP of spring-mass system under different noise situations

筆者研究5種不同工況下，無噪聲以及不同噪聲情況下(5%，10%，15%，20%)的模態(tài)辨識結果。給出了不同工況下Hankel矩陣奇異值歸一化及奇異值百分比，如圖6所示。由理論解可知，六自由度系統(tǒng)奇異值定階指標應為12，在無噪聲情況下，奇異值指標在階次為12發(fā)生明顯突變，此時奇異值百分比接近為100%，與理論相符。虛線框為突變區(qū)域，隨著噪聲不斷增加，突變不斷減緩，而百分比不斷減低，但在階次12時，不同工況下百分比仍大于80%。此時，在噪聲情況下難以根據(jù)奇異值突變的特點確定模態(tài)的階次。給出了模型在4種噪聲工況下RSVP值的變換曲線，如圖7所示。可以看出，4種不同噪聲工況下模態(tài)定階指標(RSVP)都在階次12后出現(xiàn)最低，故系統(tǒng)的模態(tài)階次為12，與理論值相符，說明在噪聲工況下，該定階方法可以準確地判定系統(tǒng)真實的模態(tài)階次。

表2和表3給出了在4種噪聲工況下(工況1,2,3,4分別為5%,10%,15%,20%噪聲工況)，傳統(tǒng)ERA和改進后算法辨識得到的頻率與阻尼比誤差值。

由表2可以看出，傳統(tǒng)ERA在剔除虛假模態(tài)的情況下，ERA與改進后的算法結果的頻率誤差都較小。表3反映兩種算法在不同工況下阻尼比誤差。分析可知，高頻段的阻尼辨識誤差相對較大，阻尼比受噪聲的干擾較大。一般表現(xiàn)隨著噪聲的增加，阻尼辨識誤差增大。同時，改進后的算法比傳統(tǒng)ERA辨識阻尼比誤差較小，識別精度有顯著提高。

采用六自由度系統(tǒng)，通過模態(tài)疊加法，ci(i=1,2,…，7)的取值使第1階模態(tài)阻尼比為0.30%，其余模態(tài)阻尼比均為0.20%。仍以m5點輸出響應為例，加入5%噪聲信號。采用傳統(tǒng)ERA及改進后算法對含有5%噪聲FRF信號進行模態(tài)參數(shù)識別，二種辨識算法的穩(wěn)定圖如圖8所示。曲線為5%噪聲下FRF通過穩(wěn)定圖結果對比。分析可知，改進后算法比傳統(tǒng)ERA的穩(wěn)定圖更簡潔，能有效剔除虛假模態(tài)。改進后算法的辨識阻尼比與理論值的誤差較小，說明阻尼比識別精度更高。

圖7 不同噪聲情況下模態(tài)階次指標Fig.7 Mode order indicators under different noise situations

表2 模態(tài)頻率辨識結果誤差Tab.2 Error of Identified modal frequencies

表3 模態(tài)阻尼比識別誤差Tab.3 Error of Identified modal damping ratio

圖8 5%噪聲下穩(wěn)定圖比較Fig.8 Comparison of two stability diagrams with 5% noise

該工況下Hankel矩陣奇異值歸一化及奇異值百分比如圖9所示。RSVP值的變換曲線如圖10所示?？梢钥闯?，該工況下模態(tài)定階指標(RSVP)在階次12后出現(xiàn)最低，故系統(tǒng)的模態(tài)階次為12，與理論值相符，說明在結構呈小阻尼特性下，該定階方法仍能準確的判定系統(tǒng)真實的模態(tài)階次。傳統(tǒng)ERA已改進后算法辨識結果，如表4所示。分析可知，傳統(tǒng)ERA在剔除虛假模態(tài)的情況下，針對小阻尼特性結構的情況，ERA與改進后的算法結果的頻率誤差都較小。同時，改進后的算法比傳統(tǒng)ERA辨識阻尼比誤差較小，識別精度更高。

圖9 5%噪聲下奇異值及百分比Fig.9 Normalized singular values and SVP with 5% noise

圖10 5%噪聲下模態(tài)定階指標Fig.10 Mode order indicators with 5% noise

2.2 某三廂車排氣系統(tǒng)模態(tài)試驗

研究某三廂車排氣系統(tǒng)，主要由前消聲器、后消聲器兩部分組成，如圖11所示。將模型用彈性繩懸掛以模擬自由—自由邊界條件。模態(tài)試驗的激勵控制、響應測量、數(shù)據(jù)處理和模態(tài)分析，都在LMS Test Lab 11B軟件下完成。

圖11 某三廂車排氣系統(tǒng)傳感器分布Fig.11 The sensor placement of a sedan car exhaust system

圖12 某三廂車排氣系統(tǒng)物理模型Fig.12 Physical model of a sedan car exhaust system

根據(jù)文獻[15]中的方法和經(jīng)驗選取傳感器安裝位置，采取單點激勵多點測量移動傳感器固定錘擊點的測量方式，選取了21個測點作為三向傳感器的安裝位置，以編號20點作為激勵點，z方向進行敲擊，如圖12所示。利用LMS系統(tǒng)測得結構的加速度響應，取14點信號為例，分別給出兩種辨識算法的穩(wěn)定圖，如圖13所示。

表4 模態(tài)參數(shù)識別結果Tab.4 Identification results of modal parameters

圖13 排氣系統(tǒng)穩(wěn)定圖比較Fig.13 Comparison of two stability diagrams of a exhaust system

圖14 排氣系統(tǒng)奇異值及百分比Fig.14 Normalized singular values and SVP of a exhaust system

分析可知，改進后算法比傳統(tǒng)ERA穩(wěn)定圖更加簡潔，可以有效剔除虛假模態(tài)，便于確定模態(tài)階次；以LMS系統(tǒng)試驗值作為參考值，頻率點顏色表征識別阻尼比誤差的大小。通過對比分析，傳統(tǒng)算法識別阻尼比誤差整體范圍在10%～20%，誤差較大，而改進后算法阻尼比誤差整體范圍在0.0%～6.0%，誤差相對較小，故改進后的算法阻尼比誤差更小，精度更高。14點處歸一化奇異值分百分比曲線如圖14所示。根據(jù)奇異值百分比大于80%確定一個模態(tài)的階次預設最小值為18。奇異值模態(tài)階次指標如圖15所示。在階次為22時指標將要達到最小值，故系統(tǒng)的模態(tài)的階次22。同時給出了傳統(tǒng)ERA與改進后算法辨識結果，表5為表示LMS系統(tǒng)試驗測量值與兩種算法辨識結果。

圖15 排氣系統(tǒng)模態(tài)階次指標Fig.15 Mode order indicators of a exhaust system

表5 模態(tài)參數(shù)辨識結果Tab.5 Identification results of modal parameters

3 結 論

1) 根據(jù)奇異值分解原理對特征系統(tǒng)實現(xiàn)算法進行改進，并提出一種判定奇異值階次的指標——奇異值百分比，根據(jù)奇異值百分比曲線的特點，首先確定一個最小預估值，再利用其增量比值變化的特點來確定模態(tài)的階次，在曲線將出現(xiàn)最低處對應著模態(tài)真實的階次。

2) 通過仿真算法辨識驗證了改進后ERA的穩(wěn)定圖更加簡潔，且阻尼比識別精度更高,研究了噪聲因素對辨識參數(shù)及模態(tài)階次的影響。結果表明，改進后算法具有更強的抗噪能力，特別是阻尼識別精度更高，同時利用RSVP方法定階效果更好。

3) 將上述方法應用在某三廂車排氣系統(tǒng)模態(tài)實驗中，識別結果良好，驗證了改進ERA和RSVP定階方法的有效性。

[1] 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.

[2] Bazan F S V. Eigensystem realization algorithm (ERA): reformulation and system pole perturbation analysis[J]. Journal of Sound and Vibration, 2004, 274(1): 433-444.

[3] Zhang Guowen, Ma Jinghua, Chen Zhuo, et al. Automated eigensystem realization algorithm for operational modal analysis[J]. Journal of Sound and Vibration, 2014, 333(15): 3550-3563.

[4] Velazquez A, Swartz R A. Output-only cyclo-stationary linear-parameter time-varying stochastic subspace identification method for rotating machinery and spinning structures[J]. Journal of Sound and Vibration, 2015, 337: 45-70.

[5] Jiang Huiming, Chen Jin, Dong Guangming, et al. Study on Hankel matrix-based SVD and its application in rolling element bearing fault diagnosis[J]. Mechanical Systems and Signal Processing, 2015, 52: 338-359.

[6] 易偉建, 劉翔. 動力系統(tǒng)模型階次的確定[J]. 振動與沖擊, 2008, 27(11): 12-16.

Yi Weijian, Liu Xiang. Order identification of dynamic system model[J]. Journal of Vibration and Shock, 2008, 27(11): 12-16. (in Chinese)

[7] 周幫友, 胡紹全, 杜強. 特征系統(tǒng)實現(xiàn)算法中的模型定階方法研究[J]. 科學技術與工程, 2009 (10): 2715-2716.

Zhou Bangyou, Hu Shaoquan, Du qiang. Study about calculation the order of model in eigensystem realization algorithm[J]. Science Technology and Engineering, 2009 (10): 2715-2716. ( in Chinese)

[8] 趙學智, 葉邦彥, 陳統(tǒng)堅. 基于奇異值曲率譜的有效奇異值選擇術[J]. 華南理工大學學報:自然科學版, 2010, 38(6):11-18.

Zhao Xuezhi, Ye Bangyan, Chen Tongjian. Selection of effective singular values based on curvature spectrum of singular values[J]. Journal of South China University of Technology: Natural Science Edition, 2010, 38(6):11-18. (in Chinese)

[9] 王樹青, 林裕裕, 孟元棟, 等. 一種基于奇異值分解技術的模型定階方法[J]. 振動與沖擊, 2012, 31(15): 87-91.

Wang Shuqing, Lin Yuyu, Meng Yuandong, et al. Model order determination based on singular value decomposition[J]. Journal of Vibration and Shock, 2012, 31(15): 87-91. (in Chinese)

[10] Cadzow J A. Signal enhancement-a composite property mapping algorithm[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1988, 36(1): 49-62.

[11] Mohanty P, Rixen D J. Modified ERA method for operational modal analysis in the presence of harmonic excitations[J]. Mechanical Systems and Signal Processing, 2006, 20(1): 114-130.

[12] Sanliturk K Y, Cakar O. Noise elimination from measured frequency response functions[J]. Mechanical Systems and Signal Processing, 2005, 19(3): 615-631.

[13] Golafshan R, Sanliturk K Y. SVD and Hankel matrix based de-noising approach for ball bearing fault detection and its assessment using artificial faults[J]. Mechanical Systems and Signal Processing, 2016, 70: 36-50.

[14] Sau-Lon James Hu, Bao Xingxian, Li Huajun. Model order determination and noise removal for modal parameter estimation[J]. Mechanical Systems and Signal Processing, 2010, 24(6): 1605-1620.

[15] 費慶國, 李愛群, 繆長青, 等. 基于主列篩選的動態(tài)測試傳感器配置方法研究[J]. 力學學報, 2008, 40(4): 543-549.

Fei Qingguo, Li Aiqun, Miao Changqing, et al. Vibration sensor placement method based on principal subset selection[J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(4): 543-549. (in Chinese)