改進VMD算法在顫振試驗信號模態參數辨識中的應用

顧文景　周麗

摘要： 提出了一種基于改進變分模態分解（Variational Mode Decomposition， VMD）的模態參數辨識算法，用于顫振試驗信號的數據處理。采用自然激勵技術提取脈沖響應信號;利用信號的先驗信息結合本文提出的適應度函數，求解最優分解參數;用參數優化后的VMD算法將信號分解為指定個數的信號分量，每個分量僅含單一頻率的振動模態;用矩陣束法識別模態參數。數值仿真和風洞試驗研究表明：改進的VMD算法可以有效分離顫振試驗信號中的密集模態，提高模態參數辨識的精度;結合顫振裕度法，有助于顫振邊界的預測。

關鍵詞： 顫振試驗; 模態參數辨識; 變分模態分解; 參數優化; 顫振邊界預測

中圖分類號： V216.2+4 ? ?文獻標志碼： A ? ?文章編號： 1004-4523（2021）02-0292-09

DOI：10.16385/j.cnki.issn.1004-4523.2021.02.009

引 言

顫振是結構在空氣動力、慣性力和彈性力耦合作用下產生的一種具有破壞性的自激振動，因而在整個飛行包線內都不允許出現顫振。由于理論分析和試驗模型不足以模擬真實的飛行環境，顫振試驗仍是飛機設計中必不可少的一個環節，以確保飛機在設計的飛行包線內不會發生顫振。

顫振試驗通常采用環境激勵的形式，利用大氣紊流對飛機結構的擾動力進行激勵，不需要額外的附加裝置，相比其他激勵方式更經濟方便。但環境激勵下的響應屬于輸入未知的振動響應信號，無法根據系統的輸入、輸出估計頻響函數或脈沖響應函數。且顫振試驗響應信號存在信噪比低、模態分布密集等特點，對信號處理方法提出了更高的要求。

對于環境激勵下的模態識別[1?3]，為簡化問題，通常認為激勵信號是高斯白噪聲，而響應信號則為平穩的隨機信號。在此假設下，可以從頻域或時域的角度，利用信號的統計特征進行系統辨識。頻域的模態識別方法通常采用經典譜估計，利用輸入、輸出的功率譜密度求解頻響函數。由于白噪聲的功率譜密度為常數，因而可以將響應信號的功率譜密度函數近似地代替頻響函數，以進行后續的模態識別。常用的時域方法大致可分為直接法和間接法兩類。直接法，如隨機子空間法[4]、ARMA分析法[5]等，通過建立參數化模型直接求解信號的模態參數;間接法，首先對信號進行處理，利用隨機減量法[6]或自然激勵技術[7]，得到其脈沖響應或相關函數后采用時域的模態識別算法進行計算。

希爾伯特?黃變換（Hilbert?Huang Transform，HHT）[8]作為一種最常用的時域信號處理技術，屬于基于經驗的數據分析方法。信號由經驗模態分解算法（Empirical Mode Decomposition， EMD）分解成一系列自適應的IMF（Intrinsic Mode Function）分量后，經希爾伯特變換即可得到信號的模態瞬時頻率和阻尼等信息。其展開基是自適應的，可以對非線性和非平穩過程產生的數據，獲得具有物理意義的表示[9]。基于HHT的數據分析算法在模態參數辨識中得到了廣泛應用[10?12];然而由于缺乏嚴格的數學背景，EMD及其相關的改進算法仍存在一些固有缺陷[13]，包括端點效應、易受噪聲干擾、存在虛假分量等問題。

VMD是一種完全非遞歸的分解算法，通過構造并求解約束變分問題，將信號分解成K個中心頻率為{ωk}的調幅?調頻信號分量[14]。相比EMD算法，VMD算法具有更嚴格的數學模型，克服了EMD算法的缺陷，可以有效分離密集模態，在信號分析、故障診斷、時間序列預測等領域取得了廣泛應用[15?17]。但是，VMD的分解參數如分解層數K和懲罰因子α的選擇對分解效果影響較大，其取值尚未有明確的理論指導。

對此，本文提出了新的評價函數，利用智能優化算法如遺傳算法、粒子群算法等進行參數優化，得到最優的分解層數K和懲罰因子α;然后，利用參數優化后的VMD算法對信號進行分解，得到一系列僅含單一模態的信號分量;最后，利用矩陣束法[18]識別各信號分量的模態參數。本文結合數值仿真算例和風洞顫振試驗，對比傳統的EMD算法和頻域內的PolyMAX法[19]，驗證所提出方法的準確性和有效性。

1.2 改進的VMD算法

針對VMD的參數選擇問題，國內外學者以包絡熵、正交系數、相關系數等參量構造評價函數，采用智能優化算法同時搜索分解層數K和懲罰因子α的最優值，在故障診斷等領域取得了一定成果[22?23]。但是，由于VMD算法計算效率的限制，導致大規模的超平面參數尋優效率低下，針對VMD分解參數的兩參甚至多參優化將耗費大量計算資源，且現有的單一指標的評價函數并不適用于振動信號的模態分解。對此，本文提出了新的VMD優化算法，對K和α的取值進行單獨尋優，并對中心頻率{ω_k^1}的初始化過程進行優化，進而提高參數優化效率以充分發揮VMD的分解性能。

在模態參數識別過程中，頻率的識別精度要高于阻尼，甚至在大多數情況下，根據信號的頻譜分析結果就能較準確地估計頻率范圍，并判斷主要模態數。而VMD算法的目標便是將信號分解成K個中心頻率為{ωk}的調幅?調頻信號，本質上便是以{ωk}為中心頻率的窄帶濾波器組。因此，本文采用簡單的峰值法預先確定信號的主要模態數M及其對應的中心頻率{fi， i=1， 2，…， M}，并將分解參數K賦值為K = M，中心頻率初始化為{ω_k^1} = {fi}。從而將原本多參優化問題簡化為僅對懲罰因子α的單參優化，加快ADMM收斂進程以提高VMD的計算效率。

針對懲罰因子α的單參優化，首先要建立評價函數，基于VMD分解結果的后驗信息對α的取值進行修正。但由于α與最終分解得到的信號分量之間沒有明確的函數關系，無法通過建立的評價函數直接求得α最優值的解析解。此外，若采用傳統優化方法（如牛頓法、單純形法等），需要遍歷整個搜索空間，加上VMD計算效率的限制，無法在短時間內完成搜索。因此，本文同樣采用智能優化算法（如遺傳算法、粒子群算法等），在超平面內搜索懲罰因子α的最優值。

本文所建立的評價函數應能準確反映理想狀態下每個分解得到的信號分量僅包含單一振動模態且沒有虛假分量及冗余模態的特征。現有的單一指標的評價函數，如包絡熵、正交系數、相關性等，雖能表征信號的稀疏性，但在實際應用過程中極易發生負優化的現象，導致α數值過大，而分解得到一組簡諧信號分量的情況。這是由于VMD窄帶濾波的特性，懲罰因子α數值越大，信號分量的帶寬越小，最終便退化成簡諧信號。簡諧信號分量相比真實的目標信號分量，在上述單一指標的評價函數中卻能得到更高的評價，從而導致負優化的情況。但這類簡諧信號屬于虛假分量，僅占原信號能量的極少部分。

因此，本文首先引入能量評價指標E=∑_（k=1）^K?（u_k^2）/x^2 ∈（0，1）表征各信號分量的能量占比之和。目標信號分量應保留原信號的大部分能量，以避免上述負優化情況的發生。針對參數優化過程中可能出現的虛假分量或冗余模態，本文再次引入相關性評價指標R=min{r（u_k，x）}∈（0，1），其中r（u_k，x）表示信號分量u_k和原始信號x的相關系數。這兩個評價指標具有相同的數量級，因此本文采用乘積運算構造聯合評價函數P = E·R，用以同時約束目標信號分量的能量及與原信號之間的相關性。為便于優化算法的求解，將其改寫成如下適應度函數

本文給出的優化目標本質是在避免出現虛假分量或冗余模態的前提下，使分解得到的信號分量具有較大的能量占比，且與原信號保持較高的相關性。在給定適應度函數的條件下，具體采用何種優化算法求解，對最后的優化結果影響不大，并不在本文討論范疇內。參數優化的具體流程如圖1所示。

2 顫振試驗信號的模態參數辨識

對于顫振試驗信號的處理，本文采用時域法。首先利用自然激勵技術提取信號的脈沖響應，對于平穩的隨機響應信號，其自相關函數與脈沖響應具有相同的數學表達式，因而可以用相關函數近似代替脈沖響應函數進行模態識別;然后利用本文提出的改進VMD算法將其分解成僅含單一模態的信號分量;最后采用矩陣束法對每個分量進行模態參數辨識。

矩陣束法的推導過程如下：

3 數值仿真研究

3.1 平板機翼模型算例

建立一個平板機翼模型，如圖3所示。機翼的后掠角為15°，半展長為140.06 mm，順氣流方向弦長為51.76 mm，機翼厚度為1.0 mm，前后緣處厚度為0，翼根AB處固支。利用Nastran計算機翼在連續大氣紊流激勵下的響應，空氣密度取為1.226 kg/m3，飛行速度為50 m/s，馬赫數為0.45，突風均方根值為1.0 m/s，采用Von Karman突風功率譜密度，采樣頻率為4096 Hz。數值計算中，取系統前4階結構模態數。

對仿真信號疊加信噪比SNR=10 dB的高斯白噪聲用以模擬測量噪聲，圖4給出了含測量噪聲的翼尖加速度響應信號的時間歷程及頻譜圖。采用自然激勵技術提取的脈沖響應信號如圖5所示。

用傳統的EMD算法對脈沖響應信號進行分解，結果如圖6所示。其中第1個IMF分量主要包含700 Hz的高頻成分;第2個IMF分量則同時包含了230和260 Hz的頻率成分，發生了模態混疊現象;第3個IMF分量包含1個低頻模態。在此算例中，EMD算法在分解過程中發生了模態混疊現象，其希爾伯特譜圖（如圖7所示）存在嚴重的鋸齒線，分解效果較差。

根據本文提出的改進VMD算法，對脈沖響應信號進行分解。首先，采用峰值法對含噪加速度響應信號的頻譜進行分析，確定共4個主要模態數，其頻率中心{fi}初步定為{40 ，230， 270， 700}Hz。用遺傳算法求解得懲罰因子的最優值αopt=2.5×104，接著令K = 4，α = 2.5×104，{ω_k^1} = {fi}，利用參數優化后的VMD算法分解脈沖響應信號，得到的4個分量及其頻譜如圖8所示。

改進的VMD算法成功分離出了4個頻率成分的信號分量，且每個分量僅包含單一頻率的振動模態，并沒有出現模態混疊現象。其希爾伯特譜圖（如圖9所示）相比EMD分解結果有了明顯改善。

用矩陣束法對EMD和改進VMD算法得到的信號分量分別進行模態參數識別。同時，采用文獻[19]提出的顫振試驗信號模態參數辨識的頻域法對翼尖加速度響應信號進行處理，將Welch法[24]估計的自功率譜密度作為頻響函數，然后用頻域內的PloyMAX法識別模態參數，最終得到的穩定圖如圖10所示。

三種算法識別的模態參數與Nastran計算的理論值對比結果如表1所示。三種算法均準確識別出了前4階模態頻率，且誤差在2%以內。EMD算法的阻尼識別誤差最大;而PolyMAX算法對二、三階模態的阻尼識別出現了大的偏差;本文提出的改進VMD算法由于成功分離出了4個單一頻率的振動模態，阻尼識別精度最高。

3.2 計算效率分析

值得注意的是，利用本文提出的適應度函數同樣可以進行多參優化，同時搜索分解層數K和懲罰因子α的最優值。為對比分析不同優化參數設置下的計算效率，本文開展了進一步的數值仿真研究，考慮下列三種優化問題：（1）同時優化K和α;（2）僅優化α，中心頻率初始化為0，即{ω_k^1}=0;（3）僅優化α，根據本文提出的方法對中心頻率進行初始化，即{ω_k^1}={fi}。將上一小節得到的翼尖加速度響應作為輸入信號，對上述每種優化問題進行多次重復求解，得到平均計算時間，以最終適應度的平均值作為評價優化效果的指標（適應度的值越小則優化效果越好）。采用遺傳算法進行尋優，其中K的取值范圍為[2， 10]，α的取值范圍為[102， 106]，除第1種優化問題因計算時間的限制僅重復計算10次外，其余兩種問題均重復計算了100次，最終結果如表2所示。

由表2可知，用本文提出的適應度函數對VMD分解參數進行單參或多參優化的效果相差無幾。在同時優化K和α的情況下，最終得到的適應度值最小，但相應的計算時間相比其余兩種情況增加了一個數量級。在約束K的情況下，對α做單參優化能極大縮短計算時間。其中，利用本文提出的方法進行參數優化所需的計算時間最短，平均為38.05 s，對于實際工程應用也在可接受的范圍內。同時優化K和α的全局尋優方式，雖然能充分發揮智能優化算法在多參超平面優化方面的優勢，但相較于優化效果的細微提升，額外增加的大量計算消耗反而顯得得不償失。相反，約束K和初始中心頻率{ω1k}的局部尋優方式，在極大縮短計算時間的同時，仍具有較好的優化效果。

4 風洞顫振試驗驗證

針對某低速顫振試驗模型，布置了28個加速度傳感器，模型示意圖及傳感器分布如圖11所示。試驗風速段為28?36 m/s，試驗顫振速度為36/m，顫振類型為小阻尼顫振型。某風速下通道404（右翼肋后）加速度響應信號的時間歷程及其頻譜如圖12所示，采樣頻率為256 Hz，采樣時間為16 s。

利用峰值法確定該通道包含三個主要模態，其頻率中心初略估計為[3.7， 7.4， 10.6] Hz。帶通濾波的頻帶設置為2?20 Hz，用于濾去低頻剛體模態和高頻噪聲信號。采用自然激勵技術提取脈沖響應信號，分別用EMD算法和本文提出的改進VMD算法對其模態分解。

如圖13所示，由于EMD算法的缺陷，仍不可避免地出現了模態混疊現象，分解效果較差。而改進的VMD算法則成功地分離出了三階低頻密集模態，如圖14所示，每個信號分量僅包含單一頻率的振動模態，提高了后續模態識別的精度。

根據本文提出的改進算法，取404通道在28?35 m/s風速下加速度響應信號識別的前三階主要模態，模態參數的識別結果如圖15和16所示。整體而言，頻率的識別結果相對穩定，阻尼的識別結果波動較大，且阻尼隨風速有逐步衰減的趨勢，符合小阻尼顫振型的特點。利用識別的模態參數，結合顫振裕度法[25]做顫振邊界預測，預測結果如圖17所示。

顫振裕度法利用結構的模態參數構造顫振預測判據，線性擬合預測判據關于動壓的變化曲線后，外推該曲線得到判據為零時的動壓即為預測的顫振點。本例中以速度的平方代替動壓，對預測的顫振點經簡單換算后即可得到預測顫振速度。采用傳統的最小二乘法進行線性擬合易受異常值的影響，導致預測結果失真。因此，本文采用穩健擬合[26]的方法，在回歸分析中自動剔除異常值，得到更為穩健的擬合結果。如圖17所示，剔除其中兩個異常點后，最終預測顫振速度為36.26 m/s。

5 結 論

（1）傳統的EMD算法由于本身的缺陷，在分離顫振試驗信號中的密集模態時，不可避免地會產生模態混疊現象，影響參數識別精度。

（2）本文對影響VMD算法分解效果的兩個關鍵參數：分解層數K和懲罰因子α進行單獨優化。利用峰值法初步確定試驗信號的主要模態數及其中心頻率范圍，將分解層數設定為模態數，初始中心頻率設置為模態的中心頻率進行后續迭代;并利用智能優化算法結合本文提出的適應度函數求解懲罰因子的最優值。數值仿真和風洞試驗算例表明，本文提出的改進VMD算法能有效分離密集模態，提高了后續模態參數識別的精度。

（3）本文提出的適應度函數同樣可用于多參優化，同時搜索分解層數K和懲罰因子α的最優值。多參優化的全局尋優方式雖能提升一定的優化效果，但由于VMD計算效率的限制，需要消耗大量的計算時間。而本文提出的局部尋優方法對K和初始中心頻率{ω_k^1}進行約束，能極大縮短計算時間，同時具有較好的優化效果。

（4）本文將改進的VMD算法結合矩陣束法對顫振試驗信號進行模態參數識別，識別結果具有較高的精度，結合顫振裕度法，有助于顫振邊界的預測。

參考文獻：

[1] 陳永高， 鐘振宇. 環境激勵下橋梁結構信號分解與模態參數識別[J]. 振動、測試與診斷， 2018， 38（6）： 1267-1274.

Chen Y G， Zhong Z Y. Signal decomposition and modal parameter identification for bridge structural under environmental excitation[J]. Journal of Vibration， Measurement & Diagnosis， 2018， 38（6）： 1267-1274.

[2] 包興先. 環境激勵下基于信號降噪的模態參數識別研究[J]. 振動與沖擊， 2014， 33（21）： 67-72.

Bao X X. Modal parameters identification based on signals denoised under ambient excitation[J]. Journal of Vibration and Shock， 2014， 33（21）： 67-72.

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

Liu Y F， Xin K G， Fan J S， et al. A review of structure modal identification methods through ambient excitation[J]. Engineering Mechanics， 2014，31（4）： 46-53.

[4] 張家濱， 陳國平. 基于隨機子空間的遞推在線模態識別算法[J]. 振動與沖擊，2009，28（8）：42-45.

Zhang J B， Chen G P. Stochastic subspace based on-line recursive modal identification method[J]. Journal of Vibration and Shock， 2009， 28（8）： 42-45.

[5] Torii J， Matsuzaki Y. Flutter margin evaluation for discrete-time systems[J]. Journal of Aircraft， 2001， 38（1）： 42-47.

[6] 陳太聰， 沈文杰. 模態辨識中隨機減量技術的實用改進[J]. 振動.測試與診斷， 2019， 39（6）： 1153-1159+1355.

Chen T C， Shen W J. Practical improvement of the random decrement technique in modal identification[J].Journal of Vibration， Measurement & Diagnosis， 2019， 39（6）： 1153-1159+1355.

[7] 鐘軍軍， 董 聰. 環境激勵下識別結構模態自然激勵-時域分解法[J]. 振動與沖擊， 2013， 32（18）： 121-125.

Zhong J J， Dong C. Natural excitation technique-time domain decomposition algorithm for structural modal identification[J].Journal of Vibration and Shock，2013， 32（18）： 121-125.

[8] Huang N E， Shen Z， Long S R， et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings A， 1998， 454 （1971）： 903-995.

[9] 楊永鋒， 吳亞峰. 經驗模態分解在振動分析中的應用[M]. 北京：國防工業出版社， 2013.

Yang Y F， Wu Y F. Applications of Empirical Mode Decomposition in Vibration Analysis[M]. Beijing：National Defense Industry Press， 2013.

[10] 練繼建， 榮欽彪， 董霄峰， 等. 抑制模態混疊的HHT結構模態參數識別方法研究[J]. 振動與沖擊， 2018， 37（18）： 1-8.

Lian J J， Rong Q B， Dong X F， et al. Structural model parameter identification method based on an improved HHT for suppressing mode mixing[J]. Journal of Vibration and Shock， 2018， 37（18）： 1-8.

[11] 李 揚， 周 麗， 楊秉才. 飛機紊流激勵響應的模態參數識別[J]. 振動工程學報， 2016， 29（6）： 963-970.

Li Y， Zhou L， Yang B C. The modal parameters identification upon atmospheric turbulence excitation response[J]. Journal of Vibration Engineering， 2016， 29（6）： 963-970.

[12] Peng Z K， Tse P W， Chu F L. An improved Hilbert?Huang transform and its application in vibration signal analysis[J]. Journal of Sound and Vibration， 2005， 286（1-2）：187-205.

[13] Huang N E， Shen S S P. Hilbert-Huang Transform and Its Applications[M]. World Scientific， 2005.

[14] Dragomiretskiy K， Zosso D. Variational mode decomposition[J]. IEEE Transactions on Signal Processing， 2014， 62（3）：531-544.

[15] 唐貴基， 王曉龍. 變分模態分解方法及其在滾動軸承早期故障診斷中的應用[J]. 振動工程學報， 2016， 29（4）： 638-648.

Tang G J， Wang X L. Variational mode decomposition method and its application on incipient fault diagnosis of rolling bearing[J]. Journal of Vibration Engineering， 2016， 29（4）： 638-648.

[16] 江春冬， 王景玉， 杜太行， 等. 基于變分模態分解算法的單通道無線電混合信號分離[J]. 上海交通大學學報， 2018， 52（12）： 1618-1626.

Jiang C D， Wang Y J， Du T H， et al. Separation of single channel radio mixed signal based on variational mode decomposition[J]. Journal of Shanghai Jiao Tong University， 2018， 52（12）： 1618-1626.

[17] Lahmiri S. Comparing variational and empirical mode decomposition in forecasting day-ahead energy prices[J]. IEEE Systems Journal， 2017， 3（11）： 1907-1910.

[18] Hua Y， Sarkar T K. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[J]. IEEE Trans. on ASSP， 1990， 38（5）： 814-824.

[19] Zeng J， Kukreja S L. Flutter prediction for flight/wind-tunnel flutter test under atmospheric turbulence excitation[J]. Journal of Aircraft， 2013， 50（6）： 1696-1709.

[20] Paternina M R A， Tripathy R K， Mendez A Z， et al. Identification of electromechanical oscillatory modes based on variational mode decomposition[J]. Electric Power Systems Research， 2019， 167： 71-85.

[21] Wang Y H， Yeh C H， Young H V， et al. On the computational complexity of the empirical mode decomposition algorithm[J]. Physical A， 2014， 400： 159-167.

[22] Wang X B， Yang Z X， Yan X A. Novel particle swarm optimization-based variational mode decomposition method for the fault diagnosis of complex rotating machinery[J]. IEEE/ASME Transactions on Mechatronics， 2018， 23（1）： 68-79.

[23] 李舒適， 王豐華， 耿俊秋， 等. 基于優化VMD的高壓斷路器機械狀態檢測[J]. 電力自動化設備， 2018， 38（11）： 148-154.

Li S S， Wang F H， Geng J Q， et al. Mechanical state detection of high voltage circuit breaker based on optimized VMD algorithm[J]. Electric Power Automation Equipment， 2018， 38（11）： 148-154.

[24] 宋 寧， 關 華. 經典功率譜估計及其仿真[J]. 現代電子技術， 2008，31（11）： 159-161.

Song N， Guan H. Classical power spectrum density estimation and its simulation[J]. Modern Electronics Technique， 2008，31（11）： 159-161.

[25] Zimmerman N H， Weissenburger J T. Prediction of flutter onset speed based on flight testing at subcritical speeds[J]. Journal of Aircraft， 1964， 1（4）： 190-202.

[26] 胡清華， 汪 運. 考慮數據噪聲的魯棒回歸建模方法綜述[J]. 西北大學學報（自然科學版）， 2019， 49（04）： 496-507.

Hu Q H， Wang Y. A review of robust regression modeling approaches with noise[J]. Journal of Northwest University （Natural Science Edition）， 2019， 49（04）： 496-507.

Modal parameter identification based on optimized variational mode decomposition and its application in signal processing of flutter test

GU Wen-jing， ZHOU Li

（State Key Laboratory of Mechanics and Control of Mechanical Structures， Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China）

Abstract： A modal parameter identification method applicable to flutter test data is proposed based on optimized variational mode decomposition （VMD）. Firstly， the natural excitation technique （NExT） is employed to extract impulse response signal from the test data. Then， the decomposition parameters are optimized by using the prior information of the test data combined with the proposed new fitness function. Finally， the target signal is decomposed into multiple monocomponents that each contains an independent oscillation mode. The matrix pencil method is adopted to identify the modal parameters. Numerical simulations and the wind-tunnel flutter test demonstrate the effectiveness of the proposed algorism in separating close modes of flutter test data. While associated with the flutter margin method， the optimized VMD can help provide an accurate flutter boundary prediction.

Key words： flutter test; modal parameter identification; variational mode decomposition; parameter optimization; flutter boundary prediction

作者簡介： 顧文景（1994?），男，博士研究生。電話：（025）84891722;E-mail：wenjinggu@nuaa.edu.cn