計及噪聲激勵的模態參數識別方法

夏遵平, 王 彤

(南京航空航天大學機械結構力學及控制國家重點實驗室，江蘇 南京，210016)

引 言

振動模態參數識別被廣泛應用于結構動態分析領域，如航空航天、大型船舶、高速列車及其他重型裝備的研制與生產、大型橋梁的驗收與監測等，是結構產品從動態設計到生產和維護過程中不可或缺的技術[1-3]。一般可將基于測試的模態參數識別方法歸納為二大類[4]：一是基于結構激勵和響應的傳統試驗模態分析(Experimental Modal Analysis，EMA)法[5]；二是僅基于結構響應的運行模態分析(Operational Modal Analysis，OMA)法[6-7]。EMA方法雖然參數識別的可靠性較高，但通常需要機械結構處于非工作狀態和較小噪聲環境中測得激勵與響應數據[8-9]。OMA方法不需要測得系統的輸入數據，利用機械結構自身運行或周圍環境噪聲激勵下的響應數據識別結構模態參數[10]。但OMA方法依然存在一些不足，如：由于采用自然激振，不能控制激振力的頻譜成分和施加位置，可能導致某些模態不能被很好地測試出來，存在模態遺漏的風險[11]；激振力無法測得，沒有很好的辦法對模態振型進行歸一化。

近年來，一種機械結構工作狀態下的試驗模態測試技術逐漸發展起來。例如，飛機飛行顫振試驗中，除未知的氣動力激勵外，再施加可測的舵面激勵[12]，將EMA方法與OMA方法結合起來，形成了含已知激勵的運行模態分析(Operational Modal Analysis in presence of eXogenous inputs，OMAX)法[13]。該方法被認為結合了EMA和OMA優點，既可識別出機械結構真實工作狀態下的模態參數，又可彌補OMA方法可靠性差的缺陷。目前，一種實現OMAX的方法是基于頻響函數(Frequency Response Function, FRF)和噪聲響應的半功率譜密度(Half Power Spectral Density，HPSD)函數[14]來識別結構模態參數。該方法通過獲取功率譜密度(Power Spectral Density，PSD)函數的正頻率部分，使得到的HPSD函數與FRF函數具有相同的極點數，認為兩者可以融合。OMAX的另一類實現方法是基于FRF和傳遞率函數(Transmissibility Function，TF)[15]來識別結構模態參數，該方法試圖通過計算各點響應PSD比值之差的倒數來構造可以與FRF融合的函數。目前研究表明，該法虛假極點較多，有時阻尼比和振型識別容易失效[4, 16]。

本文提出一種計及噪聲激勵的試驗模態分析(Experimental Modal Analysis with consideration of noise eXcitation, EMAX)方法。根據含噪聲激勵的無參數模型估計出僅含噪聲響應的功率譜函數，并對其進行Hilbert變換估計出噪聲響應的類頻響函數(Analogous FRF，AFRF)。研究表明，AFRF具有與傳統FRF一致的系統極點和相似的數學表達[17-18]，可以與FRF融合，采用EMA方法識別出模態參數[19]。以彎扭二自由度機翼的數值仿真和GARTEUR飛機模型試驗為例，分析所提AFRF的性質。通過與僅EMA法、僅OMA法和OMAX法對比，驗證本文方法的可靠性和有效性。

1 理論背景

1.1 類頻響函數估計

結構系統同時含有已知激勵和未知激勵的輸入-輸出模型如圖1所示。

圖1 含未知激勵的輸入-輸出模型Fig.1 Input-output model with ambient excitation

圖中，fk和yk分別表示可測激勵及其產生的系統響應，fu和yu分別表示不可測激勵及其產生的系統響應，y表示測得的系統的響應。H′和H″分別表示可測力與不可測力對應的頻響函數，可能只是總體FRF的一個元素或一行或一列。上述各變量的頻譜關系可表示為

Y=Yk+Yu=H′Fk+H″Fu

(1)

式中Y,Yk,Yu,Fk和Fu分別為y,Yk,yu,fk及fu的傅里葉變換，且只有Y和Fk已知(試驗測得)。根據上式，未知響應可表示為

Yu=Y-H′Fk

(2)

則得到其功率譜函數為

(3)

式中 上標H表示復數矩陣的共軛轉置，Sff和Syy分別表示測得的激勵與響應的自功率譜，而Syf表示測得的響應與激勵的互功率譜。對于響應中含噪聲的情況，公式(3)中的頻響函數可采用H1估計[20]方法獲得，即

(4)

將式(4)代入公式(3)得到基于輸入輸出的噪聲響應功率譜函數為

(5)

由于頻響函數的相位與幅值存在Hilbert變換的對應關系[15-16]，進而可得到原點頻響函數的真實相位與該點響應自譜的關系為

(6)

式中 H為Hilbert變換符號。得到未知激勵下的原點AFRF為

(7)

同時注意到，第i點激勵、第j點響應的跨點與原點AFRF的關系為

(8)

(9)

式中 上標T,H及*分別表示轉置、共軛轉置及共軛;Nr為系統的極點數;r為系統極點的階次;Φr為第r階振型;Kr為第r階模態參與向量或含有比例常數的模態參與向量;λr為系統的第r階極點。

1.2 參數識別

(10)

式中Ωr=exp(-jωrΔt)為多項式基本項，r為模型階次，αr和βr為參數向量。設αr和βr的集合為

θr=[αrβr]T

(11)

(12)

式中Wo(ωk)為加權函數。Xo和Yo分別為：

(13)

(14)

式中 ?表示Kronecker乘法。

對公式(12)求平方和并取其實部的跡，得

(15)

式中 Re(·)表示取實部，上標H表示共軛轉置。參數向量α可由最小二乘解出，然后求得其伴隨矩陣的特征值為

(16)

式中Λ即為系統的極點組成的對角矩陣，據此可求出系統的固有頻率和阻尼比;V為參與向量，據此便可求出振型。

(17)

2 算例分析

2.1 仿真算例

圖2所示的彎扭二自由度機翼模型被廣泛應用于飛機飛行顫振分析的仿真驗證。該模型的平動和轉動可模擬真實機翼的彎曲和扭轉兩種振動模態。設機翼質量m=5 kg，轉動慣量JG=3.34 kg·m2，平動彈簧剛度為k1=3.6×104N/m，轉動彈簧剛度為kr=2.0×104N/m；偏心距e= 0.05 m，兩測點分別位于l1=0.8 m,l2=1 m；兩階模態分別加入阻尼比為1.2%和0.86%的結構阻尼。在圖中的1號點和2號點處施加隨機噪聲激勵，同時在1號點處施加0～40 Hz的線性掃頻激振。設掃頻激勵和仿真響應的采樣率為102.4 Hz，采樣時間為800 s，得到時域激勵與響應如圖3所示。根據上節中的Hilbert變換理論，估計出未知噪聲激勵下的AFRF，并將其與理論FRF比較，如圖4所示。

圖2 機翼物理模型Fig.2 The physical model of the aircraft wing

圖3 時域激勵與響應(局部)Fig.3 Input and output signal (zoom in)

圖4 理論頻響函數與類頻響函數Fig.4 The theoretical FRF and analogous FRF

由圖4可知，根據公式(7)和(8)估計出的噪聲激勵的AFRF與理論FRF具有相同的相位，且兩者的幅頻曲線也保持著一致性，僅相當于在數值上相差一個比例常數。因此，可認為僅根據噪聲響應功率譜估計出的AFRF與理論FRF有相同的性質及相似的表達形式。分別基于EMA (僅FRF)法[5,19]、OMA (僅HPSD)法[19]、OMAX (FRF+HPSD)法[13-14]和本文提出的EMAX法識別出結構的模態頻率與阻尼比，如表1和2所示。

表中結果表明，在本算例中，各方法識別出的固有頻率均具有較小的相對誤差，而前3種傳統方法對阻尼比的識別誤差較大。這是因為阻尼參數比固有頻率參數更為敏感，易受信噪比、測試的完備性等因素影響。由誤差對比可知，EMAX方法識別參數的精度明顯高于其他方法。這是因為與EMA方法相比，EMAX方法采用噪聲響應擴充了頻響函數，使其更完備；與HPSD函數相比，基于Hilbert變換的AFRF不僅具有相同的極點數，零點-極點的組成也近似(曲線形狀相似)。各方法計算過程中的穩態圖如圖5所示，圖中曲線為復模態指示因子曲線。其中，EMA方法由于人工激勵的不完備，僅能得到一列頻響函數，即僅有一條復模態指示因子曲線。

盡管我國“保基本、強基層、建機制”的醫療體制改革已使基層醫療衛生機構硬件設施明顯改善，但是，國家衛計委發布的公報顯示，2015年，全國醫療衛生機構總診療人次達77.0億人次。其中，三級醫院的診療人次和入院人數增長分別為7.14%和8.55%，而基層醫療衛生機構出現了負增長，增長率分別為-0.46%和-1.39%[11]。

表1 估計出的固有頻率比及相對誤差

表2 估計出的模態阻尼比及相對誤差

圖5 模態穩定圖Fig.5 Modal stabilization diagrams

2.2 實 驗

如圖6所示的GARTEUR (Group of Aeronautical Research and Technology in EURope)飛機模型是歐洲航空科技研究集團于20世紀90年代中期設計制造的。該模型具有真實飛機的主要振動模態特征，因此，在航空界得到廣泛應用。

圖6 實驗裝置Fig.6 Experimental setup

將基于GARTEUR的振動測試分成兩組：第一組，作為參考，在無噪聲激勵的情況下，利用脈沖激勵測得加速度FRF；第二組，在脈沖激勵的同時，用激振器模擬環境噪聲激勵。測試過程中采樣頻率均設為128 Hz，每個測點實施8次錘擊脈沖激勵，每次脈沖觸發采樣時間為16 s(即頻率分辨為0.0625 Hz)，得到典型的加速度響應信號如圖7所示。

圖7 響應信號Fig.7 Output signals

由公式(6)～(8)估計出測得的噪聲響應的類頻響函數，并將其與真實頻響函數對比，如圖8所示。由圖中兩種曲線對比可知，AFRF的相位與真實FRF的相位一致(部分頻段相位相反)，并且注意到兩者的幅頻曲線形狀也相似，即在數值上僅相差一個比例常數。基于第一組實驗數據，采用多參考點的最小二乘復頻域法(polyLSCF)[19]識別出的模態參數作為參考值，其計算穩態圖如圖9所示；基于第二組實驗數據，分別采用EMA法、OMA法、OMAX法和EMAX法識別出結構50 Hz內的結構模態頻率(如表3所示)和阻尼比(如表4所示)，各方法的計算穩態圖如圖10所示。

圖8 GARTEUR模型的頻響函數與類頻響函數Fig.8 The FRF and AFRF of GARTEUR

圖9 PolyLSCF識別的穩態圖(基于第一組測試數據)Fig.9 Stabilization chart of polyLSCF(based on the first test data)

參考頻率/Hz(第一組實驗)EMA方法/Hz相對誤差%OMA方法/Hz相對誤差/%OMAX法/Hz相對誤差/%EMAX方法/Hz相對誤差/%6.126.120.006.120.006.120.006.130.1616.1216.11-0.0616.09-0.1916.11-0.0616.11-0.0635.7435.740.0035.740.0035.73-0.0335.750.0339.5539.550.0039.550.0039.550.0039.550.0039.94--39.940.00--39.93-0.0347.2347.21-0.0447.230.0047.240.0247.20-0.0648.3848.420.0848.390.0248.410.0648.40-0.08

表4 估計出的模態阻尼比及相對誤差

圖10 4種參數識別方法的穩態圖(基于第二組測試數據)Fig.10 Stabilization charts of four methods (based on the second test data)

對比表3和4中的參數識別結果發現，由于施加的人工激勵不完備，使得激勵出的第5階模態較弱，同時基于HPSD的OMAX方法并未有效改善這一問題，進而導致EMA法和OMAX法僅識別出了6階模態參數，漏掉了第5階模態。此外，注意到OMA方法和OMAX方法識別出的阻尼比誤差較大，這是由于相比于頻率，阻尼參數更為敏感，對識別方法、測試數據質量等有較高要求，而本實驗的響應數據中含有脈沖和隨機兩種不同屬性的成分，而且數據長度有限，因而影響上述方法對模態阻尼比的識別精度。與之相比，EMAX法不僅識別出了全部的結構模態，而且提高了模態阻尼比的識別精度。EMAX法識別出的50 Hz內的模態振型如圖11(a)～(g)所示。

圖11 識別出的模態振型及其MAC值Fig.11 Estimated mode shapes and MAC values

圖11(h)為振型與參考振型的模態置信準則(Modal Assurance Criterion，MAC)矩陣，其對角值均接近于1，非對角值均小于0.3，說明兩者具有較好的匹配度，即采用EMAX方法識別出的振型是可靠的。綜上結果表明，本文提出的EMAX方法是有效且可靠的。

3 結 論

(1)基于噪聲響應估計出的類頻響函數與真實頻響函數具有相同的性質，當噪聲激勵的頻域特性滿足平直譜時，兩者僅相差一個常數，因而更適于聯合。

(2)與OMAX方法依據OMA理論不同，本文提出的參數識別方法立足于EMA，利用結構自身工作或周圍環境激勵，彌補由于人工激勵不完備而可能存在的問題。

(3)該方法可用于機械設備的工作狀態或橋梁的交通狀態的試驗模態分析，測試過程中實施較少的人工激勵也能得到較可靠的數據，節省了試驗成本。

需要說明的是，本文中基于功率譜Hilbert變換的類頻響函數的推導要求環境或結構工作提供的噪聲激勵近似滿足隨機分布。對于多數模態測試情況，考慮較窄的識別頻帶，這種近似條件是合理的。

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

[2] Pintelon R, Schoukens J. Vandersteen G, Frequency domain system identification using arbitrary signals [J]. IEEE Transactions Automatic Control, 1997, 42(12): 1717—1720.

[3] Cheung Y K, Leung A Y T. Finite Element Methods in Dynamics [M]. Kluwer Academic Publishers, Dordrecht, 1991.

[4] Wang T, Celik O, Catbas F N. Damage detection of a bridge model based on operational dynamic strain measurements [J]. Advances in Structural Engineering, 2016, 19(9): 1379—1389.

[5] Verboven P, Cauberghe B, Guillaume P. Improved total least squares estimators for modal analysis [J]. Computers and Structures, 2005, 83: 2077—2085.

[6] Brincker R，Zhang L M, Anderson P. Modal identification from ambient response using frequency domain decomposition[C]. Proceeding of the 18th IMAC, 2000:625—630.

[7] Reynders E, Maes K, Lombaert G, et al. Uncertainty quantification in operational modal analysis with stochastic subspace identification: validation and applications [J]. Mechanical Systems and Signal Processing, 2016, 66-67: 13—30.

[8] 孫鑫暉, 郝木明, 李振濤, 等. CMIF方法中模態參數不確定性的計算[J]. 振動工程學報, 2013, 26(3): 380—386.

Sun X H, Hao M M, Li Z T, et al. Uncertainty calculation in complex modal indictor function identification method [J]. Journal of Vibration Engineering, 2013, 26(3): 380—386.

[9] 王 彤, 張令彌. 多頻段擬合的正交多項式方法及其Matlab工具箱[J]. 振動工程學報, 2003, 16(4):493—497.

Wang T, Zhang L M. FA Multiple frequency band fitting method using orthogonal polynomials and corresponding matlab toolbox [J]. Journal of Vibration Engineering, 2003, 16(4): 493—497.

[10] Schoukens J, Pintelon R. Measurement of frequency response functions in noisy environments [J]. IEEE Transactions on Instrumentation and Measurement 1990, 6(6): 373—377.

[11] Brinker R, Ventura C E. Introduction to Operational Modal Analysis [M]. John Wiley & Sons Ltd., 2015, 10—11.

[12] Zhang L M, Brincker R, Andersen P. An overview of operational modal analysis: major development and issues[C]. Proceedings of the 22nd International Modal Analysis Conference, 2004: 179—190.

[13] Vanlanduit S, Guillaume P, Cauberghe B,et al. On-line identification of operational loads using exogenous inputs [J]. Journal of Sound and Vibration, 2005, 285(1-2): 267—279.

[14] Troyer T D, Runacres M. Frequency-domain modal analysis in the OMAX framework[C]. Proceedings of the 28th International Modal Analysis Conference, Jacksonville, Florida USA, February, 2010: 465—476.

[15] Devriendt C, Troyer T D, Sitter G D, et al. Transmissibility-based operational modal analysis for flight flutter testing using exogenous inputs[J]. Shock and Vibration, 2012, 19(5): 1071—1083.

[16] 張永年. 基于傳遞率函數的運行狀態模態分析方法及軟件實現[D]. 南京: 南京航空航天大學，2014.

Zhang Y N. Transmissibility-based operational modal analysis and software implementation [D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2014.

[17] Agneni A, Crema L B, Coppotelli G. Output-only analysis of structures with closely spaced poles [J]. Mechanical Systems and Signal Processing, 2010, 24(5): 1240—1249.

[18] 夏遵平, 王 彤, 張永年. 基于Hilbert譜變換函數的運行模態分析方法[J]. 地震工程與工程振動, 2014, 34(6): 97—102.

Xia Z P, Wang T, Zhang Y N. Operational modal analysis based on hilbert transformation [J]. Earthquake Engineering and Engineering Dynamics, 2014, 34(6): 97—102.

[19] Peeters P, Van Der Auweraer H, Guillaume P, et al. The PolyMAX frequency-domain method: a new standard for modal parameter estimation [J]. Shock and Vibration, 2004, 11(3-4): 395—409.