杭 曉晨，蔣 令 聞，顧明 華，吳 邵 慶，費慶 國

(1.江蘇省工程力學(xué)分析重點實驗室，江蘇南京210096；2.東南大學(xué)工程力學(xué)系，江蘇南京210096；3.中國航空工業(yè)集團公司中航飛機研發(fā)中心，陜西西安710089)

頻域分解模態(tài)識別方法的阻尼識別精度研究

杭 曉晨1，2，蔣 令 聞1，2，顧明 華3，吳 邵 慶1，2，費慶 國1，2

(1.江蘇省工程力學(xué)分析重點實驗室，江蘇南京210096；2.東南大學(xué)工程力學(xué)系，江蘇南京210096；3.中國航空工業(yè)集團公司中航飛機研發(fā)中心，陜西西安710089)

在模態(tài)參數(shù)識別中，模態(tài)阻尼比往往比固有頻率、振型參數(shù)更難被準確識別。為了提高頻域分解模態(tài)識別方法的阻尼識別精度，分析了算法中產(chǎn)生誤差的主要因素，提出了采用自回歸模型譜估計代替原有周期圖譜估計的改進方法。開展了仿真算例研究，以白噪聲激勵下結(jié)構(gòu)的振動響應(yīng)為輸入條件，使用頻域空間域分解法進行模態(tài)參數(shù)識別。算例結(jié)果驗證了采用基于自回歸模型的譜估計代替原有的周期圖法譜估計，既保留了對密頻模態(tài)的分辨能力，又提高了阻尼參數(shù)的識別精度，且在采樣數(shù)據(jù)長度較短的情況下具有優(yōu)越性。

參數(shù)識別；阻尼識別；頻域分解法；譜估計

引 言

頻域分解(Frequency Domain Decomposition，F(xiàn)DD)類方法是一類運行模態(tài)分析［1］(Operational Modal Analysis，OMA)方法，通過測量結(jié)構(gòu)在環(huán)境激勵下的振動響應(yīng)，得到能夠反映真實動力學(xué)特性的模態(tài)參數(shù)。在模態(tài)參數(shù)中，阻尼比往往會比固有頻率、振型更難被準確識別［2-4］。因此有必要對此類方法的阻尼識別精度展開研究，分析阻尼識別誤差的來源，改善識別效果。FDD方法最早由Brincker等在復(fù)模態(tài)指示函數(shù)的基礎(chǔ)上提出［5］，在滿足白噪聲激勵假設(shè)和小阻尼假設(shè)時，能夠很好地識別模態(tài)頻率和振型，且對密集模態(tài)的分辨率較高，但存在無法識別模態(tài)阻 尼的缺 點［6］；為 了解決 這一問 題，Brincker等隨后提出了增強頻域分解法［7］(Enhance Frequency Domain Decomposition，EFDD)，在FDD的基礎(chǔ)上將頻域信息轉(zhuǎn)化到時域，然后在時域中使用對數(shù)衰減(Logarithmic Decrement，Logdec)法得到模態(tài)阻尼比。然而EFDD方法在時-頻域轉(zhuǎn)化的過程中必然帶來數(shù)據(jù)截斷，造成了阻尼識別的不精確；王彤等［8］提出頻域空間域分解法(Frequency and Spatial Domain Decomposition，F(xiàn)SDD)，引入增強功率譜密度(Enhance Power Spectrum Destiny，EPSD)的概念，直接在頻域內(nèi)對EPSD進行最小二乘擬合后獲得模態(tài)頻率和阻尼。Zhang等［9］利用FSDD法實現(xiàn)了對幾個實際工程結(jié)構(gòu)的模態(tài)識別。

總的來說，F(xiàn)DD類方法使用奇異值分解［10］(Singular Value Decomposition，SVD)技術(shù)，將信號空間與噪聲空間分離開來，抗噪聲干擾能力較強，且具有較好的密頻分辨能力。本文重點關(guān)注FDD類方法的阻尼識別精度問題，首先從算法上分析帶來誤差的幾個因素，并針對其中影響較大的可控因素——功率譜估計方法進行探討，提出使用基于自回歸(Autoregressive，AR)模型的譜估計方法進行改進，開展了兩個仿真算例驗證，得到了一些結(jié)論。

1 頻域分解模態(tài)識別方法的理論分析

FDD類方法都是基于如下的輸入、輸出功率譜密度函數(shù)(Power Spectrum Destiny，PSD)關(guān)系

式中 Gxx(jω)，Gyy(jω)分別為輸入和輸出PSD矩陣；H(jω)為頻響函數(shù)(Frequency Response Function，F(xiàn)RF)矩陣，F(xiàn)RF可以用部分分式形式表示式中 N為模態(tài)階數(shù)，λr為第r階極點，Rr為第r階留數(shù)矩陣，R*r為第r階模態(tài)振型與模態(tài)參與向量轉(zhuǎn)置的乘積，即Rr=φrγTr。假設(shè)輸入為一白噪聲，則Gxx(jω)為一實常數(shù)對角陣，將式(2)代入到式(1)，化簡得到

式中 Ar為相關(guān)的留數(shù)項，在小阻尼情況下，r階模態(tài)附近Ar可以被簡化為

式中 dr=Cγr為一實數(shù)。在某一頻率處，當只有少數(shù)幾個模態(tài)(一般只有一個或兩個)貢獻顯著時，定義這些模態(tài)集合為Sub(ω) ，響應(yīng)PSD矩陣可以寫成

式中 Ui為奇異值向量矩陣，Si為奇異值組成的對角矩陣。

在輸出的PSD圖中，如果譜峰附近僅有第r階模態(tài)起主導(dǎo)作用，那么方程(5)可以近似為

對比式(6)和(7)，第一個奇異值向量ur即為第r階模態(tài)振型，對應(yīng)的奇異值具有單自由度特性，通過比較奇異值曲線峰值附近模態(tài)振型的模態(tài)置信度(Modal Assurance Criterion，MAC)確定此單自由度特性曲線的范圍，進而確定固有頻率。

EFDD法是FDD法的改進方法，它解決了FDD法無法識別阻尼參數(shù)的問題。EFDD選取了奇異值曲線峰值附近MAC值較大的頻段，對其作逆傅立葉變換運算(Inverse Fast Fourier Transform，IFFT)轉(zhuǎn)換到時域，獲得了近似的單自由度相關(guān)函數(shù)曲線，然后用對數(shù)衰減法識別模態(tài)頻率和阻尼。然而，EFDD法使用IFFT時會使數(shù)據(jù)截斷，得到的相關(guān)函數(shù)是不精確的，導(dǎo)致阻尼識別產(chǎn)生偏度誤差。

FSDD法直接在頻域中識別模態(tài)參數(shù)，因此避免了對奇異值曲線作IFFT帶來的誤差。將奇異值曲線作為指示函數(shù)，并定義第r階模態(tài)的增強PSD

聯(lián)系式(7)，在第r階模態(tài)主導(dǎo)的譜峰附近窄帶，EPSD具單自由度特性，且

用最小二乘法擬合此EPSD曲線，求出極點λr，進而得到模態(tài)頻率和阻尼。

FDD類方法的基本原理可以概括為圖1所示。分析以上3種方法的基本理論，可以得到FDD類方法參數(shù)識別的誤差的來源：

1)輸出PSD估計?Gyy與真實PSD間存在誤差。響應(yīng)PSD估計作為FDD類識別算法的最初輸入，對最終識別結(jié)果的影響較大；

2)使用激勵白噪聲假設(shè)以及小阻尼假設(shè)進行近似產(chǎn)生誤差。這些假設(shè)是整個識別算法成立的前提，但在激勵時間足夠長，模態(tài)阻尼較小時，這些假設(shè)引起的誤差可以忽略；

3)由單自由度特性曲線得到模態(tài)參數(shù)時的數(shù)值誤差。EFDD法中IFFT運算、Logdec法均會產(chǎn)生此類誤差。而FSDD直接在頻域中進行識別，數(shù)值誤差僅存在于最小二乘擬合過程中，當擬合數(shù)據(jù)較為平滑時，可以認為這種誤差很小。

圖1 FDD類方法基本流程圖和誤差分析Fig.1 Flowchart and error analysis of FDD methods

2 功率譜估計理論

FDD類方法均采用周期圖(Periodogram)法來估計輸出響應(yīng)的PSD［5］。假設(shè)N點的離散輸出響應(yīng)信號為x(n) ，n=0，1，2，…，N-1，周期圖法PSD估計［11］

周期圖法PSD估計的均值

式中 W(ω)為數(shù)據(jù)截斷時采用的窗函數(shù)w( n) 的頻域描述。可見，周期圖法PSD估計為真實譜S(ω)與窗譜W(ω) 的卷積，只有當N→∞時，E((ω))→S(ω) ，才是無偏估計。而實際情況N總是一個有限值，所以周期圖法PSD估計與真實譜之間存在偏差

改進的周期圖法對總數(shù)據(jù)進行有重疊分段，然后對每一段數(shù)據(jù)選用合適的窗函數(shù)以減少譜泄漏，最后對各分段估計PSD后進行平均。假設(shè)總數(shù)據(jù)量為N，分成L段，每段有M個數(shù)據(jù)，由于有數(shù)據(jù)重疊，N≠M×L，運算過程可以表示為：

改進的周期圖法最大的優(yōu)點在于減小方差，當N→∞時，Var→0，為漸進一致估計?？偟膩碚f，周期圖譜估計方法具有較高的頻域分辨率，但估計方差大，譜泄漏嚴重，其改進方法采用分段、加窗等措施，降低了頻域分辨率，提高了估計效果，但在數(shù)據(jù)長度一定的條件下，依然不滿足無偏性、一致性估計條件。FDD類方法以此PSD估計為基礎(chǔ)，模態(tài)參數(shù)的識別精度不能保證。尤其是在顫振試飛試驗等工程應(yīng)用領(lǐng)域，由于試驗條件限制，采樣數(shù)據(jù)不可能很長，PSD估計誤差愈加明顯。

本文提出基于AR模型譜估計代替周期圖譜估計的改進方法。AR模型譜估計是一種參數(shù)化模型估計方法，它不存在對離散序列進行時-頻域轉(zhuǎn)化的過程，避免了周期圖法中存在的泄漏、方差等問題，在合適的模型階數(shù)下，PSD估計效果更好［11］。AR模型可用如下差分方程表示

式中u( n) 是均值為零方差為σ2的白噪聲序列；p為AR模型階數(shù)；ak，k=1，2，…，p為對應(yīng)的p階模型參數(shù)。通過AR系統(tǒng)的傳遞函數(shù) H( z) ，可以推得信號 x( n) 的PSD估計

可見基于AR模型的PSD估計關(guān)鍵在于求解模型參數(shù)。式(19)經(jīng)過相關(guān)函數(shù)的變換可以得到Y(jié)ule-Walker方程［12］

式中 r(p)為信號 x ( n) 的相關(guān)函數(shù)。求解Yule-Walker方程即可得到模型參數(shù)，Burg法［13］可以避免對離散信號進行自相關(guān)估計，因此本文之后的AR模型參數(shù)估計均使用Burg法。

3 算例研究

采用兩個算例進行仿真研究，第1個為二自由度彈簧-振子系統(tǒng)，用于說明更換譜估計方法后能否保持FDD類方法對密頻模態(tài)的識別能力；第2個為大展弦比機翼仿真模型，用于比較兩種PSD估計方法對阻尼識別精度的影響。算例均采用FSDD法識別模態(tài)參數(shù)。

3.1 二自由度彈簧-振子模型

為了驗證基于AR模型的PSD估計方法對密頻模態(tài)的識別能力，建立如圖2所示的兩階頻率靠近的二自由度系統(tǒng)。m1=10 kg，m2=1.03 g，k1=1 000 k N/m，k2=100 N/m，通過模態(tài)疊加法確定兩階密頻模態(tài)頻率為f1=49.61 Hz，f2=50.40 Hz，c1和c2的取值使兩階模態(tài)阻尼比均為2%。

圖2 二自由度彈簧-振子系統(tǒng)Fig.2 Two-DOFspring oscillator

在質(zhì)量塊m1上施加白噪聲激勵力，在m1，m2上采集振動位移響應(yīng)。響應(yīng)信號的采樣頻率fs=200 Hz，分別取N1=2 000和N2=10 000個采樣點，來研究兩種采樣長度條件的識別。

取m2點輸出響應(yīng)為例，圖3為使用AR模型和使用周期圖法得到的PSD對比圖?？梢娭芷趫D法譜估計曲線不平滑，且第二階模態(tài)峰不突出，幅值上已經(jīng)失真。而基于AR模型估計的PSD曲線光滑，能夠很好地分辨兩階峰。兩種譜估計方法對應(yīng)的奇異值曲線對比如圖4所示。

圖3 兩種估計方法的響應(yīng)PSD比較Fig.3 Comparison of two PSD estimations

模態(tài)參數(shù)識別結(jié)果如表1和2所示，經(jīng)過分析可得：

1)雖然周期圖法對響應(yīng)信號的PSD估計效果較差，但經(jīng)過FSDD識別后依然能夠區(qū)分兩階密頻模態(tài)，驗證了 FDD類方法易于區(qū)分密頻模態(tài)的優(yōu)點?；贏R模型進行PSD估計，保留了較高的密頻分辨率；

圖4 奇異值曲線比較Fig.4 Comparison of two singular value plots

表1 采樣數(shù)據(jù)N1=2 000時2種方法識別結(jié)果Tab.1 Identification results of the two methods when data length N1=2 000

表2 采樣數(shù)據(jù)N2=10 000時2種方法識別結(jié)果Tab.2 Identification results of the two methods when data length N2=10 000

2)采用周期圖法和AR模型PSD估計方法都能獲得良好的頻率識別精度。然而，在數(shù)據(jù)長度較短時，前者對阻尼的識別誤差很大，后者阻尼識別精度較高。

3.2 大展弦比機翼模型

將一大展弦比機翼進行有限元建模，它的計算模態(tài)分析結(jié)果作為準確值。通過在有限元模型上施加隨機激勵，進行瞬態(tài)響應(yīng)分析，提取時程響應(yīng)數(shù)據(jù)，來模擬運行模態(tài)試驗。最后對響應(yīng)數(shù)據(jù)進行識別得到模態(tài)參數(shù)，與計算模態(tài)分析的結(jié)果進行對比驗證。

機翼具體參數(shù)如表3所示，幾何模型和有限元模型如圖5所示，邊界條件為根部固支，模態(tài)阻尼比設(shè)為2%。機翼有限元模型的前7階模態(tài)分析結(jié)果如表4所示。

表3 大展弦比機翼模型參數(shù)Tab.3 Properties of the wing with high aspect ratio

表4 有限元模型模態(tài)分析結(jié)果Tab.4 Modal analysis results of FEM

在機翼模型上進行多點的均值為零、方差為1的高斯白噪聲激勵，激勵時間50 s。響應(yīng)采集通道設(shè)置為沿機翼表面均勻分布的12個測點(如圖5所示)，采樣頻率fs=200 Hz，采集時間取T1=20 s，T2=50 s，以便分析數(shù)據(jù)長度對識別結(jié)果的影響。

分別使用改進的周期圖法、AR模型法進行輸出譜的估計，固有頻率的識別結(jié)果如表5所示，模態(tài)阻尼的識別結(jié)果如表6所示。

圖5 大展弦比機翼的幾何模型Fig.5 Geometry of the wing with high aspect ratio

表5 模態(tài)頻率識別結(jié)果Tab.5 Identification results of modal frequencies

表6 模態(tài)阻尼識別結(jié)果Tab.6 Identification results of modal damping

分析表5和6數(shù)據(jù)可得到初步的結(jié)論：

1)兩種方法對模態(tài)頻率都具有良好識別精度，誤差都在2%以內(nèi)。而對于阻尼參數(shù)的識別，同數(shù)據(jù)長度條件下，基于AR模型的PSD估計方法識別效果更好；

2)周期圖法在短數(shù)據(jù)阻尼識別精度不理想，增大數(shù)據(jù)強度可以提高識別精度，這是因為更大的數(shù)據(jù)長度能使周期圖PSD估計趨向于無偏估計；

3)對于AR模型PSD估計方法來說，較短的輸入數(shù)據(jù)(4 000個采樣點)也能獲得良好的阻尼識別精度。

由于作為激勵的白噪聲信號的產(chǎn)生具有隨機性，一次仿真試驗的結(jié)果往往不能說明結(jié)論的普遍性。因此重復(fù)進行了25次仿真，每次產(chǎn)生一組隨機激勵信號，響應(yīng)數(shù)據(jù)點定為10 000個。取第三階模態(tài)為例，只考慮識別其模態(tài)阻尼比，識別誤差如圖6所示。

由圖6可知：采用周期圖法譜估計，模態(tài)阻尼識別誤差波動較大，最大誤差可能超過10%；而采用基于AR模型的譜估計方法，模態(tài)阻尼的識別誤差基本在5%以下，證明了采用自回歸模型譜估計方法能獲得更好的阻尼識別效果。

圖6 第三階模態(tài)阻尼識別誤差Fig.6 Damping identification error of the third mode

4 結(jié) 論

本文分析了FDD類模態(tài)參數(shù)識別算法產(chǎn)生誤差的主要因素，提出采用基于自回歸模型的譜估計方法改進響應(yīng)PSD的估計精度，進而提高模態(tài)識別精度，采用兩自由度系統(tǒng)和大展弦比機翼進行了有限元仿真研究，主要結(jié)論如下：

1)功率譜作為FDD類方法的輸入條件，其估計的準確性是影響模態(tài)參數(shù)識別精度主要因素。傳統(tǒng)的周期圖估計不滿足估計的無偏性和一致性，尤其在短數(shù)據(jù)條件下，模態(tài)阻尼的識別精度難以保證；

2)使用基于AR模型的PSD估計方法代替周期圖法，既保留了FDD類方法對密頻模態(tài)的識別能力，又提高了模態(tài)阻尼的識別精度，且識別效果受數(shù)據(jù)長度影響不大。這是因為AR模型估計法避免了對離散信號進行時頻域轉(zhuǎn)換，不存在“譜泄漏”帶來的偏度誤差。但是參數(shù)化方法必須選取合適模型階數(shù)p，本文通過對比周期圖法得到的奇異值曲線，峰值位置相近時確定p的取值。這種方法類似于時域模態(tài)參數(shù)識別中常用的穩(wěn)定圖和頻譜圖聯(lián)合定階的方法［14］，具有較好的定階效果；

3)對于顫振試飛試驗等只能獲得有限數(shù)據(jù)，但又對阻尼識別精度要求較高的情況，本文提出的基于譜估計改進的FDD類方法具很大的優(yōu)越性。

［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］Andrianne T，Dimitriadis G.Damping identification of lightly damped linear dynamic systems using commonbase proper orthogonal decomposition［J］.Mechanical Systems and Signal Processing，2012，28：492—506.

［3］Ding Y，Law S S.Structural damping identification based on an iterative regularization method［J］.Journal of Sound and Vibration，2011，330(10)：2 281—2 298.

［5］Brincker R，Zhang L，Andersen P.Modal identification from ambient responses using frequency domain decomposition［A］.Proceedings of 18th International Modal Analysis Conference［C］.2000：625—630.

［6］Zhang L.An overview of major developments and issues in modal identification［A］.Proc.IMAC XXII［C］.Detroit(USA)，2004：1—8.

［7］Brincker R，Ventura C，Andersen P.Damping estimation by frequency domain decomposition［A］.Proceedings of 19th International Modal Analysis Conference［C］.2001：698—703.

［8］王彤，張令彌.運行模態(tài)分析的頻域空間域分解法及其應(yīng)用［J］.航空學(xué)報，2006，27(1)：62—66.Wang T，Zhang L M.Frequency and spatial domaindecomposition for operational modal analysis and its application［J］.Acta Aeronautica Et Astronautica Sinica，2006，27(1)：62—66.

［9］Zhang L，Wang T，Tamura Y.Afrequency-spatial domain decomposition(FSDD)method for operational modal analysis［J］.Mechanical Systems and Signal Processing，2010，24(5)：1 227—1 239.

［10］Henry E R，Hofrichter J.Singular value decomposition：application to analysis of experimental data［J］.Essential Numerical Computer Methods，2010，210：81—138.

［11］Stoica P，Moses R L.Spectral Analysis of Signals［M］.Upper Saddle River，NJ：Pearson/Prentice Hall，2005.

［12］Marple L.Anew autoregressive spectrum analysis algorithm［J］.IEEE Trans.on Acoustics，Speech and Signal Processing，1980，ASSP-28(4)：441—454.

［13］Brockwell P J，Davis R A.Introduction to Time Series and Forecasting［M］.Springer Science&Business Media，2006.

［14］Rainieri C，F(xiàn)abbrocino G.Automated output-only dynamic identification of civil engineering structures［J］.Mechanical Systems and Signal Processing，2010，24(3)：678—695.

Accuracy of modal damping identification using frequency domain decomposition method

HANG Xiao-chen1，2，JIANG Ling-wen1，2，GU Ming-hua3，WU Shao-qing1，2，F(xiàn)EI Qing-guo1，2
(1.Jiangsu Key Laboratory of Engineering Mechanics，Nanjing 210096，China；2.Department of Engineering Mechanics，Southeast University，Nanjing 210096，China；3.Research and Development Center of AVIC Aircraft Corporation，Ltd.，Aviation Industry Corporation of China，Xi'an 710089，China)

Modal damping ratios are more difficult to be accurately identified compared with nature frequencies and mode shapes in modal parameters identification.To improve accuracy of the damping identification using frequency domain decomposition method，an error analysis is conducted and a spectrum estimation method based on autoregressive model is presented instead of the previous periodogram method.Two simulation examples are implemented in which the measurement response data of structure under white noise excitation is adopted as input to identify the modal parameters using the frequency and spatial domain decomposition.Results indicate that with the proposed method the ability of distinguishing closely spaced modes is retained，and identified damping ratios can be more accurate.Moreover，the proposed method is superior to periodogram method when the sampling data is not long enough.

parameters identification；damping identification；FDD；spectral estimation

TB123；O329

A

1004-4523(2015)04-0518-08

10.16385/j.cnki.issn.1004-4523.2015.04.003

杭曉晨(1990—)，男，博士研究生。電話：(025)83790168；E-mail：hangxioachen@seu.edu.cn

2014-01-11；

：2014-06-17

國家自然科學(xué)基金資助項目(10902024)；教育部新世紀優(yōu)秀人才支持計劃資助項目(NCET-11-0086)；航空科學(xué)基金資助項目(20090869009)；江蘇省普通高校研究生科研創(chuàng)新計劃資助項目(CXZZ13_0084)