夏遵平, 王 彤

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

引 言

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

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

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

1 理論背景

1.1 類頻響函數(shù)估計

結(jié)構(gòu)系統(tǒng)同時含有已知激勵和未知激勵的輸入-輸出模型如圖1所示。

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

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

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

(1)

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

Yu=Y-H′Fk

(2)

則得到其功率譜函數(shù)為

(3)

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

(4)

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

(5)

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

(6)

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

(7)

同時注意到，第i點激勵、第j點響應(yīng)的跨點與原點AFRF的關(guān)系為

(8)

(9)

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

1.2 參數(shù)識別

(10)

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

θr=[αrβr]T

(11)

(12)

式中Wo(ωk)為加權(quán)函數(shù)。Xo和Yo分別為：

(13)

(14)

式中 ?表示Kronecker乘法。

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

(15)

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

(16)

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

(17)

2 算例分析

2.1 仿真算例

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

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

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

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

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

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

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

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

表2 估計出的模態(tài)阻尼比及相對誤差

圖5 模態(tài)穩(wěn)定圖Fig.5 Modal stabilization diagrams

2.2 實 驗

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

圖6 實驗裝置Fig.6 Experimental setup

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

圖7 響應(yīng)信號Fig.7 Output signals

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

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

圖9 PolyLSCF識別的穩(wěn)態(tài)圖(基于第一組測試數(shù)據(jù))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 估計出的模態(tài)阻尼比及相對誤差

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

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

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

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

3 結(jié) 論

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

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

(3)該方法可用于機械設(shè)備的工作狀態(tài)或橋梁的交通狀態(tài)的試驗?zāi)B(tài)分析，測試過程中實施較少的人工激勵也能得到較可靠的數(shù)據(jù)，節(jié)省了試驗成本。

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

[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方法中模態(tài)參數(shù)不確定性的計算[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] 張永年. 基于傳遞率函數(shù)的運行狀態(tài)模態(tài)分析方法及軟件實現(xiàn)[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譜變換函數(shù)的運行模態(tài)分析方法[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.