优化递归变分模态分解及其在非线性信号处理中的应用*

许子非 岳敏楠 李春†

1) (上海理工大学能源与动力工程学院,上海 200093)

2) (上海市动力工程多相流动与传热重点实验室,上海 200093)

经验模态分解一类的递归算法所产生的模态混淆和端点效应将导致所获物理信息失真,变分模态分解可改善这些问题.但因其需预设参数,对信号分解精度影响显著,为此,提出采用目标信号功率谱峰值所对应的频率以初始化变分模态分解所需中心频率,借鉴经验模态分解递归模型,基于能量截止法将变分模态分解改进为递归模式算法,并采用粒子群优化算法对具有带宽约束能力的惩罚因子进行最优取值,构成优化递归变分模态分解.通过对比分析经验模态分解,集成经验模态分解及优化递归变分模态分解在分解信号时的计算精度；研究传统变分模态分解与优化递归变分模态分解在处理实际振动信号时计算速率.结果表明：优化递归变分模态分解在处理目标信号时精度最高,与原分量相关性达99.9%；与集成经验模态分解对比,可由低至高将信号分解至不同频段,物理意义更加清晰且不产生虚假模态；处理实际非线性信号时,优化递归变分模态分解无需预设分解模态个数,相比于传统变分模态分解,计算速率高12.5%—18.5%.

1 引 言

非平稳信号具有随机性强、稳定性差与物理信息混叠等特点[1].为探寻非稳定信号中所含信息,在结构分析、故障诊断乃至医学领域[2-5],相关学者对信号(时间序列)分解技术已展开诸多研究.其中,应用较为广泛的有小波分解(wavelet transform,WT)、经验模态分解(empirical mode decomposition,EMD)、集合经验模态分解(ensemble empirical mode decomposition,EEMD)和局部均值分解(local mean decomposition,LMD)等[6-9].WT方法采用缩放、平移窗函数将原信号分解为若干小波叠加的形式,但易产生无物理意义的虚假谐波,影响结果分析[10].EMD与EEMD均属于递归模式分解,EEMD在处理信号时,通过对信号加入高斯噪声以抑制EMD处理信号时易产生的模态混淆的问题,可有效处理非线性、非平稳信号分解；但在EEMD迭代过程中,由于包络线估计误差的叠加,导致引发端点效应及模态混淆[11].LMD作为递归算法,可自适应分解无规则信号至多个含物理意义的乘积函数分量,较小波分析有更好的自适应性,但无法分离频率相近的信号.以上方法因求解本质属于递归分解,均存在包络线估计误差随迭代而累计,从而导致模态混淆及端点效应,影响对物理信息的提取及分析.

为此,Dragomiretskiy和Zosso[12]提出基于Hilbert变换、Wiener滤波及频率混合等概念的新型分解方法—变分模态分解(viriational mode decomposition,VMD),该方法因具有完备的数学理论基础,通过求解约束变分模型代替原有的递归模式,从而抑制模态混淆、避免端点效应及解决相近频率难分离的问题[13].

采用VMD对信号进行求解时,涉及模态分解数、惩罚因子、保真系数及收敛条件等参数的预设.研究表明,其中惩罚因子α与模态分解数K对信号分解精度影响最为显著[14].已有学者针对惩罚因子α与模态分解数K的何种取值对分解效果如何影响做了相关研究.文献[15]以VMD分解所获分量间的相关系数作为约束,从而确定模态分解数K的取值,并将改进VMD与深度置信网络相结合,实现故障预警.文献[16]以包络熵作为VMD参数优化的目标函数进行寻优取值,结果具有良好的分解精度,并将其应用于滚动轴承故障诊断中.文献[17]将VMD与多尺度排列熵相结合,在生物组织变性识别中获得较优的聚类效果与分类性能.Baldini等[18]比较VMD与EMD在分解无线电信号的准确性,并采用机器学习对信号进行分类,得出基于VMD分解的信号具有更高的准确性.Chen等[19]分别采用VMD与EMD对轴承振动信号进行分解,以能量熵作为特征值构建特征向量,并采用支持向量机对特征向量进行分类,结果表明以VMD-能量熵构建的特征样本,使得分类结果具有更高的准确性,但未考虑VMD分解时预设参数对VMD分解精度的影响.为此,Cui等[20]基于目标信号瞬时频率,对VMD算法K值设定进行改进,指出该方法可提升VMD分解精度,并成功用于智能负载建模与去噪,但忽略了VMD分解时以约束带宽的参数—惩罚因子对结果的影响.

以上研究采用VMD处理信号时,或考虑模态分解数对分解效果的影响,而忽略了惩罚因子对带宽约束,致使分解效果不佳；或采用局部寻优的方法对惩罚因子与模态分解数进行选择,未考虑预设参数间的互交性；又或基于分量特征值作为分解效果判别标准,以获取全局优化的参数组合,但缺乏一定的物理意义.

因此,本文针对VMD处理信号前需预设惩罚因子α与模态分解数K,且α与K取值组合对分解精度有显著影响的问题,提出一种优化递归VMD算法(optimized recursive variational mode decomposition,ORVMD),该算法借鉴EMD递归分解时可自适应获取模态分解数K的优点,以目标信号功率谱最大值对应的频率作为初始频率,并基于能量差追踪法将传统VMD改进为递归VMD,并采用粒子群优化算法获取惩罚因子α的最优取值,得到最佳α,K组合,以保证分解所得分量具有更多的物理信息.对比EMD,EEMD及ORVMD在分解信号时的精确性,分析EEMD与ORVMD在处理非线性调制信号时的可靠性,研究传统VMD与ORVMD解决实际振动信号时的时效性,发现ORVMD具有更高的计算精度与效率,亦可保留更多的物理信息,以期为信号处理、结构分析及故障诊断等方面奠定理论基础,提供处理方法.

2 VMD算法

VMD算法是基于Wiener滤波、Hilbert变换与外差解调所形成的一种分解算法,因其具有完备的数学基础,采用VMD处理信号时可有效避免EMD类算法导致的模态混淆及端点效应.与属递归性质的EMD类算法不同,VMD属约束变分性问题,通过假设模态分量的中心频率,将约束模态带宽的过程转化为约束变分问题,求解变分模型以实现模态分解.VMD算法中涉及的约束变分模型为

式中 {uk}={u1,u2,···,uk}与{ωk}={ω1,ω2,···,ωk}分别为第k 个模态分量及其相应中心频率,共K个模态；δ(t) 为单位脉冲函数；j表示虚数单位；* 表示卷积运算；∂t为偏导运算；f为目标信号.

引入惩罚因子α和Lagrange乘子Λ以求解变分约束问题.所得增广Lagrange表达式如下：

采用交替方向乘子算法(alternate direction method of multiplers,ADMM)更新迭代求解(2)式的鞍点,在频域内迭代更新 uk,ωk及Λ.

VMD将信号分解为K个模态分量,步骤如下：

2) uk和 ωk分别由(3)式和(4)式迭代更新,

式中τ为保真系数；∧表示傅里叶变换；n 为迭代次数；

3) 通过(5)式更新Λ,

4) 重复步骤2和3,至满足迭代终止条件,终止条件由(6)式给出,

式中ε为判别精度,且 ε＞ 0；上标n 为迭代步数,下标k 表示当前模态数；

5) 输出K个模态分量.

3 优化递归VMD算法

与EMD等递归分解算法相比,VMD具有较高的准确性与稳定性.但VMD算法需给出预设模态分解数K,且K的取值影响分解精度.因此,借鉴EMD递归思想,提出一种递归VMD算法.借鉴能量差追踪法以设定停止条件,实现递归VMD算法.该方法既拥有VMD算法中可抑制求解包络线迭达误差而导致失真的优点,也可实现自适应分解,无需预设模态分解数K.

3.1 递归算法

3.1.1 改进变分模型

基于传统VMD算法,给定初始模态分解数K=1的约束变分模型,获得有限带宽本征模态函数(bandwidth intrinsic mode function,BIMF),将(1)式转化为(7)式.

引入惩罚因子α和Lagrange乘子Λ后,(2)式将变为(8)式,(5)式和(6)式被转换为(9)式和(10)式,其余步骤与传统VMD算法相同.

采用ADMM算法时,中心频率 ωk需预计算.而Dragomiretskiy等[12]所提出三种初始化中心频率的方法(分别为：赋予初值零、线性化及随机初值化)均有一定的局限性.其中线性化以获取初值具有较好的准确性与稳定性,但当 K=1 时,中心频率 ωk=0 ,与赋予初值零情况相同,从而导致线性化失败.因此,提出以目标信号功率密度谱(power density spectrum,PSD)最大值所对应的频率作为 ωk初值.

3.1.2 停止准则

Huang等[21]与Damerval等[22]所提出的停止准则均不适用于VMD算法,故采用Cheng等[23]所提出的能量差追踪法应用于本文所提出的改进变分模型.基于能量差追踪法的VMD停止准则如下.

假设目标信号 f (t) 分解后所得BIMF分量uk(t)具有正交性,由(11)式给出,其总能量表达式由(12)式给出.

式中E为原信号能量.

因各BIMF的正交性,(12)式可表示为

式中 E1+E2+···+En为各BIMF分量的能量.

当BIMF分量完全正交时,各BIMF能量之和 Etotal应与原信号能量E相等,

但当BIMF分量不完全正交时,Etotal与E之间存在误差 Eerr,

能量误差 Eerr越接近零,表明分解效果越好,与原分量符合度更高,所得模态分量颇含物理意义.因此,令 Eerr的绝对值 |Eerr|为停止准则,当|Eerr|小于收敛阈值时停止迭代,以获取最佳模态分解数K.

3.1.3 递归VMD算法

基于以上提出的改进变分模型,结合能量差追踪法确定收敛条件,递归VMD算法流程图如图1所示,具体计算步骤如下：

1) 输入目标信号,通过计算目标信号的PSD,以PSD最大值对应的频率最为中心频率 ωinitial；

图1 递归VMD流程Fig.1.The recursive VMD diagram.

2) 采用步骤1获取的中心频率进行迭代,并预设模态分解数 K=1 ,获取(7)式的变分模型,通过(9)式和(10)式分别获取 uk与 ωk；

3) 采用步骤2获取的 uk作为BIMF分量,将f-uk作为新目标信号并重复步骤1与2；

4)对于分解所得BIMF分量,通过(15)式计算能量误差 |Eerr|,当 |Eerr|小于收敛阈值时停止分解并获取全部BIMF分量.根据计算经验,采用|Eerr|≤[(0.7～ 2.0)E]/100为收敛阈值；

5) 输出全部BIMF分量.

3.2 PSO算法

粒子群优化算法(particle swarm optimization,PSO)是由Kennedy 和Eberhart[24]受鸟群觅食行为启发所提出的.在PSO算法中,每个粒子具有相应的速度与位置以调整自身的状态,粒子的位置代表待优化问题的潜在解.PSO算法的数学描述如下.

在D维空间中,存在由m 个粒子组成的种群,第i 个粒子在D维空间中的位置向量为Xi=(xi1,xi2,···,xiD)T,其速度向量为Vi=(vi1,vi2,···,viD)T,第i 个粒子的最佳位置由向量Pi=(pi1,pi2,···,piD)T表示,种群最佳位置由向量Pg=(pg1,pg2,···,pgD)T表示,速度与位置的更新由(16)式和(17)式所示.

式中 i =1, 2,···,m ；d =1, 2,···,D；k 为当前进化代数；c1与 c2为学习因子；ω为惯性权重；η为[0,1]内随机数.

3.3 ORVMD算法

以上所提出的递归VMD方法虽可自适应获取分解模态数K,但惩罚因子α作为约束带宽的参数,其取值也严重影响了信号处理的准确性.带宽的变化影响分解所获各模态能量的变化,从而进一步影响递归模态分解数K或信号分解失效[25].因此,以能量误差 |Eerr|为寻优过程的适应度函数,将适应度函数最小取值作为寻优目标,基于PSO算法对改进递归VMD算法进行参数寻优,获取最佳分解参数组合,其流程如图2所示.

1) 初始化PSO算法各项变量并确定寻优过程中的适应度函数 |Eerr|,以惩罚因子α及模态分解数K作为粒子位置,并随机初始化粒子移动速度；

图2 基于PSO优化改进递归VMD参数流程Fig.2.The process of using PSO to optimize recursive VMD parameter.

2) 在不同粒子位置条件下对信号进行递归VMD处理,获该位置下信号能量误差 |Eerr|,以此作为对应粒子的适应度函数；

3) 对比各粒子适应度函数大小(优劣),若有粒子的适应度函数 |Eerr|小于当前最小适应度函数,则对粒子进行更新；

4) 通过(16)式与(17)式更新粒子属性；

5) 判断粒子是否满足种群进化停止条件,若不满足则重复步骤2继续寻优,直到满足最大种群进化预设值,输出最佳粒子,即为最优参数组合 [α,K].

基于PSO优化算法及改进迭代VMD算法,称为优化递归VMD算法(optimized recursive variational mode decomposition,ORVMD),该改进算法能获取较精确的惩罚因子α以约束带宽,并自动获得分解模态个数K,确定最佳参数组合 [α,K].

4 ORVMD算法验证

为验证ORVMD分解信号的可靠性及有效性,以两个低频信号与间断高频冲击信号组成仿真复合信号.仿真复合信号 f (t) 与其分量表达式由(18)式所示：

各分量及其合成信号的波形图如图3所示.为验证算法ORVMD的有效性,分别采用EMD,EEMD及ORVMD对信号 f (t) 进行分解,对比其结果.

首先采用EMD对仿真信号进行分解,分解结果如图4所示.由图4可知,EMD将信号分解为5个IMF分量及1个残余分量,分解效果较差,IMF1与IMF2均出现模态混淆现象,IMF3勉强具有一定的物理意义.

EEMD引入白噪声以克服EMD分解时易出现对极值点选判出错的影响.选取噪声幅值系数0.02,集合平均次数50,采用EEMD分解结果如图5所示.

由图5可知,EEMD分解出7个IMF分量及1个残余分量.其中IMF1,IMF3与IMF4分别为原信号中的间断高频信号及两个低频信号,与原信号分量的相关系数分别为0.981,0.942和0.944,表明EEMD对分解信号的物理意义还原性较好,较EMD具有更好的准确性,但仍产生许多虚假模态不具有物理性.

图4 EMD分解结果 (a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) resFig.4.The results of EMD：(a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) res.

图5 EEMD分解结果 (a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) resFig.5.The results of EEMD：(a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) res.

采用本文提出的ORVMD对信号进行分解,其中分解模态数 K=3,α=500 ,分解结果如图6所示.

如图6所示,原信号3个分量被较好地还原至3个IMF分量,其中IMF3为间断高频信号,两个低频信号分别为IMF1与IMF2.各IMF分量与对应原分量的相关系数分别为0.998,0.991与0.999,较EEMD具有更高的精度,可获取更具物理意义的信号分量.

图6 ORVMD分解结果 (a) IMF1；(b) IMF2；(c) IMF3Fig.6.The results of ORVMD：(a) IMF1；(b) IMF2；(c) IMF3.

5 非线性信号分解

处理实际信号时,目标信号多为调制信号,故为验证ORVMD在处理实际信号也具有高效性与准确性,以如下数学模型模拟滚动轴承早期局部损伤信号 x (t)[26]：

式中 s (t) 为周期性冲击信号.信号衰减指数C=1000,系统共振频率 fn=4000Hz.τk表示第k 次冲击相对于特征周期的小波动,该随机波动服从均值为零,标准差为0.5%转频的正态分布.特征频率 fr=140Hz.U(tk) 为单位阶跃函数.加入高斯噪声信号后得信噪比为—5 dB的仿真信号.采样频率为16000 Hz,采样点N= 4096.信号时域及频谱如图7所示.

5.1 EEMD分解

采用EEMD对信号 x (t) 进行分解,其分解所得时域与频域如图8所示,共11个IMF分量及1个残余分量.

如图8所示,IMF1—IMF6五个分解分量包含良好的物理信息,将重叠在全频段的信号分解在不同频段,但低频信号仍互相重叠,无法体现轴承损伤的物理信息.且IMF7—IMF11基本属于无效信息,蕴含的物理信息较少,这说明EEMD对非线性信号的分解效果不佳.

图7 早期轴承内圈故障信号 (a) 时域；(b) 频谱Fig.7.Early inner race fault diagnosis signal.

图8 故障信号EEMD分解结果 (a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) IMF8；(i) IMF9；(j) IMF10；(k) IMF11；(l) IMF12Fig.8.The results of EEMD for fault signal：(a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) IMF8；(i) IMF9；(j) IMF10；(k) IMF11；(l) IMF12.

5.2 VMD分解

采用所提出的ORVMD方法对信号进行处理,其中模态分解数 K=9 ,优化惩罚因子α=16100,所得结果如图9所示.

由图9可知,每一阶IMF都包含一定的物理意义,原本冗杂在各个频段的信号已从低频至高频准确区分,且未出现过分解或者虚假模态.IMF1中蕴含内圈损伤特征频率,IMF5因损伤导致轴承系统发生共振,主频与谐波之间的频段间隔清晰可见,包含一定的物理信息.相比于EEMD,ORVMD具有更高的精度,不存在虚假模态,且未分解出不具物理意义的IMF分量以干扰对信号的分析.

5.3 ORVMD与传统VMD对比

因采用VMD分解时,已有研究表明模态分解数K与惩罚因子α对分解结果及精度有显著影响,在此不做分析.除模态分解数K可自适应获得外,配合PSO优化算法获取最佳参数组合,较之于传统VMD,本文所提出的ORVMD因改变模态分量的初始中心频率获取方式,以PSD最大值对应的频率作为初始中心频率,代替传统VMD线性初始化的方法,可减少迭代次数以提升计算速度.

采用西储大学轴承振动信号为例[27],对风扇端轴承内圈故障的25组信号分别进行ORVMD以及传统VMD进行分解,其中实验台转速为1730 rmp,采样频率12000 Hz,负载为3 Hp.为保证对比有意义,采用与ORVMD分解时所需参数[α,K]组合作为传统VMD分解时的输入参数.

采用ORVMD方法对25组信号进行分解,获得25组α及K组合,采用PSO算法时,惩罚因子α寻优的取值范围为[5000,20000],取值步长为100.ORVMD参数组合如表1所列.

由于实验采集信号含有不同程度的噪声,噪声对分解产生一定的干扰,导致最优模态分解数会稍有波动.由表1可得,80%的模态分解数均为12,故选取模态分解数 K=12 为传统VMD分解时的模态预设定.惩罚因子看似毫无规律,但对应实验数据采样频率可知,惩罚因子作为约束带宽的参数,均在12000附近波动,与采样频率相近(12000 Hz),ORVMD获取的25个最优惩罚因子均值为12600,但其中7次为12000,6次为11900.因此,确定α= 12000为传统VMD分解时的惩罚因子.

以 α=12000,K=12 作为传统VMD分解预设参数,分别采用传统VMD与ORVMD对上述25组实验振动信号进行分解,两种方法分解信号耗时如图10所示.

图9 故障信号ORVMD分解结果 (a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) IMF8；(i) IMF9Fig.9.The results of ORVMD for fault diagnosis：(a) IMF1；(b) IMF2；(c) IMF3；(d) IMF4；(e) IMF5；(f) IMF6；(g) IMF7；(h) IMF8；(i) IMF9.

图10 VMD与ORVMD计算耗时Fig.10.The duration of calculation for VMD and ORVMD.

表1 ORVMD参数组合Table 1.Parameter set of ORVMD.

由图10可知,改进的ORVMD算法在计算时间上比传统VMD所需时间少,虽然平均耗时仅减少12.5%—18.5%,但由于ORVMD在处理信号时调用PSO算法以获取最优惩罚因子的过程占用较大计算资源,因此,认为ORVMD比传统VMD具有更高分解精度与计算效率.

6 结 论

对于传统VMD算法需预设分解个数及惩罚因子,提出ORVMD算法,该算法可对目标信号递归并获取确定模态分解数K,并基于PSO算法对惩罚因子α进行寻优取值,获取最佳预设参数组合.同时,对比分析EMD,EEMD,VMD与ORVMD对非线性信号分解,表明所提出的ORVMD较EMD,EEMD具有更高的精准性.该方法可将物理信息分解至不同频段,且不产生虚假模态,保证了信号分解的有效性及准确性.并且,与传统VMD方法相比,ORVMD具有更高的计算效率,效率升幅在12.5%—18.5%之间.此外,还可开展以复信号、分布式信号[28-31]为研究对象的变分模态分解处理效果研究.