论文部分内容阅读
摘要： 提出了一种基于改进变分模态分解（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
关键词： 颤振试验; 模态参数辨识; 变分模态分解; 参数优化; 颤振边界预测
中图分类号： 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