基于蝙蝠算法优化VMD参数的滚动轴承复合故障分离方法

2022-10-27李军霞陈维望

振动与冲击 2022年20期

张伟，李军霞，陈维望

(1. 太原理工大学机械与运载工程学院，太原 030024；2. 太原理工大学矿山流体控制国家地方联合工程实验室，太原 030024)

滚动轴承是现代机械设备的重要组成部分，在工业生产系统中发挥着重要作用，广泛用于矿山、交通运输和医药制造等行业。由于负荷大、运行环境恶劣，滚动轴承等关键部件可能发生复合故障，与任何单一故障相比，其危害性更大。复合故障的特征提取一直以来都是故障诊断领域面临的难题[1]。因此，提供一种能够从原始振动信号中提取微弱周期特征信息，准确识别复合故障位置和类型的方法，对机械设备安全高效运行具有重要意义。

在工业生产中，机械振动响应是多个频率特征信息的叠加，需要通过信号分解和滤波等操作来提取故障特征，以便最终进行故障识别[2]。针对此问题，国外学者已经设计了多种方法，如小波分解(wavelet transform，WT)、小波包分解(wavelet packet decomposition, WPT)、经验模态分解(empirical mode decomposition，EMD)和局部平均分解(local mean decomposition，LMD)等。小波分解和小波包分解是非自适应信号分析方法，需要提前确定小波基函数。EMD可以将复杂信号自适应分解为有限个本征模函数(intrinsic mode function，IMF)，但由于模态混叠和边界效应等缺陷限制了其应用范围。为了进一步提高EMD的性能，相继提出了如集成经验模态分解(ensemble empirical mode decomposition, EEMD)[3]、互补集合经验模态分解(complementary ensemble empirical mode decomposition，CEEMD)和自适应噪声完备集合经验模态分解(complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)等方法，并用于机械故障诊断。这些方法在一定程度上缓解了模式混合等问题，但增加的白噪声不能有效消除，降低了计算效率。

受EMD的启发，Dragomiretskiy等[4]在2014年提出了一种具有坚实理论基础的完全非递归方法，称为变分模态分解(variational modal decomposition, VMD)。VMD在信号处理方面具有优异的性能，在音调检测、音调分离和噪声鲁棒性方面优于EMD，避免了模态混合、端点效应和噪声抑制等不利影响。此外，VMD已广泛应用于工程实践中。胡爱军等[5]采用参数优化VMD结合1.5维谱的方式，实现了滚动轴承复合故障的有效分离。万书亭等[6]提出基于VMD及最大相关峭度解卷积(maximum correlated kurtosis deconvolution，MCKD)的复合故障分离方法，仿真及实测信号表明该方法能够从复合故障中分离出单一故障特征。以上研究表明，采用VMD方法能够实现滚动轴承复合故障的诊断，但强噪声干扰下的复合故障诊断仍需要进一步研究。

然而，VMD的优越性能严重依赖模态数量和惩罚因子等参数的选取，这些参数对分解效果有着显著影响。利用经验或先验准则选取参数极易导致分解结果不准确，降低了VMD的分解效率。郑圆等[7]根据峭度最大原则确定分解层数K，再通过鲸鱼算法优化选择惩罚因子α，但忽略了两个参数之间的关联性。一些学者利用灰狼优化算法(grey wolf optimizer，GWO)、蚱蜢优化算法(grasshopper optimization algorithm，GOA)、布谷鸟搜索 (cuckoo search，CS)算法等智能优化算法同时确定模态数量和惩罚因子的最优值，与以上仿生算法相比，蝙蝠算法(bat algorithm, BA)采用频率调谐手段，增加了搜索到最优参数的概率，能够在全局寻优和局部寻优过程之间自由切换，避免陷入局部收敛，提高了算法的收敛性。Meng等[8]通过多个基准问题和工程设计表明，蝙蝠算法在全局优化方面具有高效性和稳定性。

在采用仿生算法优化VMD时，选取合适的目标函数至关重要，目标函数将直接影响最终分解的准确性和效率。目前常用的目标函数有局部极小熵值[9]、集合峭度、加权峭度[10]、平均包络熵[11]等，采用这些目标函数取得了不错的分解效果，但仍存在信息丢失或模式冗余等问题。针对此问题，本文提出一种新的复合影响指数，该指数能够准确识别故障并抑制噪声干扰。

为了自适应获取变分模态分解的参数，本文将最小平均复合影响指数作为目标函数，采用蝙蝠算法对参数空间进行全局寻优。对分解后的信号进行包络解调分析，提取故障特征。通过仿真信号和工程数据表明，该方法能够准确识别强噪声干扰下的滚动轴承复合故障类型。

1 基本原理

1.1 VMD算法

VMD算法通过迭代搜寻变分模型的最优解，将原始信号分解为带宽之和最小的K个IMF分量uk(t)。

其变分模型为

(1)

式中：f(t)为初始输入信号； {uk}，{ωk}为各模态信号以及中心频率；*为卷积运算。

引入二次惩罚因子α和Lagrange算子式λ(t)，将约束性变分问题变为非约束变分问题，其拉格朗日表达式[12]为

(2)

在变分模态分解中，对中心频率分别进行零值、均匀分布、随机分布初始化，采用交替方向乘子法进行多次迭代求取式(2)的最优解，不断更新其中心频率，结果发现得到的中心频率基本保持不变。在误差允许范围内，为了减少计算量，初始模态设为0矩阵，VMD算法流程如下：

步骤2进行迭代n=n+1，执行整个循环；

步骤3初始k=0，k=k+1，对所有的ω≥0，更新uk，ωk为

(3)

(4)

步骤4更新λ为

(5)

1.2 蝙蝠算法

蝙蝠算法通过模拟蝙蝠的觅食行为对优化问题的解空间进行搜索，通过自适应调节声波脉冲响度、频率，以达到寻优的目的[13]。在蝙蝠算法中，每一个蝙蝠代表一个可行解，而蝙蝠的猎物则代表优化问题的最优解。在整个解空间里，蝙蝠频率、位置和速度的更新计算式分别为

fi=fmin+(fmax-fmin)β

(6)

(7)

(8)

式中：β∈[0,1]的随机数；x*为t时刻全局搜索过程中蝙蝠的最优位置。

局部搜索时，从当前最优解集中随机选取一个解并对其施加一个随机扰动，然后在该解的邻域内进行搜索。此过程中，蝙蝠位置更新计算式为

xnew=xold+ηAt

(9)

式中：xold为所选取的最优解；η∈[-1,1]的随机数；At为所有蝙蝠在时刻为t时发出声波的平均响度。

觅食过程中，为了平衡全局搜索和局部搜索过程，蝙蝠发射的声波响度和频率需要跟随搜索猎物的进程自动调节。蝙蝠发出声波的响度Ai和发射速率ri的调节计算式分别为

(10)

(11)

式中：α为声波响度增加系数；γ为脉冲频度增强因子。

BA算法的伪代码如下：

目标方程f(x),x=(xi,…,xd)T

初始化蝙蝠种群位置xi和速度vi(i=1,2,…,n)

定义位置xi的发射频率fi

初始化脉冲发射率ri和响度Ai

while (t<最大迭代次数)

通过式(6)～式(8)更新蝙蝠的频率、速度和位置形成新的解

if(rand>ri)

选出全局最优解

围绕选定的最佳解决方案生成局部解决方案

end if

通过随机飞行产生一个新的解(见式(9))

if(rand

接受新的解

增大ri并且减小Ai(式(10)和式(11))

end if

对所有蝙蝠排序，找到目前最好的解x*

end while

后处理结果和可视化

2 方法步骤

2.1 复合影响指数

选取合适的评估标准对机械故障诊断非常重要，峭度、包络谱、功率谱、相关系数，熵等单一指标已被广泛研究。考虑到各个单一指标关注重点不同，目前已经提出了包络谱峭度(kurtosis of envelope spectrum, KES)、平方包络谱峭度(kurtosis of squared envelope spectrum, KSES)和包络熵值(envelope entropy, EE)等混合指数。然而，在工程实践中，除了背景噪声外，采集到的振动信号还包含有随机冲击和周期性谐波等干扰信号，这对目前的测量指标是一个相当大的挑战。

本文基于峭度、功率谱、相关系数等指标对故障诊断的影响，提出了一种新的复合影响指数(composite impact index, CII)，以便及时准确检测到故障冲击，其计算过程如下：

对于采集到的信号f(n)，其希尔伯特变换后的信号fH(n)计算为

fH(n)=Hilbert{f(n)}

(12)

对希尔伯特变换后的信号fH(n)进行N次采样，并进行离散傅里叶变换，进而得到信号的功率谱SH(k)

(13)

(k=0,1,…,N-1)

(14)

计算功率谱峭度(kurtosis of power spectrum, KPS)

(15)

考虑到模式与主信号之间的联系，引入相关系数，其定义为

(16)

式中：f和uk分别为原始信号和分解后的模态分量； E[·]为数学期望。

由此可得，CII的表达式为

CII=|CC|·KPS

(17)

本文提出的指数，与文献[14]提出的综合影响指数相比，该参数计算过程相对简单，省去了构造解析函数和计算包络谱等过程，加快了计算速度。该指数通过计算信号的功率谱峭度，将脉冲故障特征放大，有利于检测到信号的故障信息。考虑到相关系数的优点，该指数能够实现滤波功能，减少噪声干扰。

为了评估该指标的性能，根据第3章使用的仿真信号，构建了一个包含轴承故障冲击、随机冲击、周期谐波和高斯白噪声的模拟信号。本节将CII与一些现有指标进行比较，如包络谱峭度、平方包络谱峭度、整体峭度(ensemble kurtosis, EK)和包络熵值，为了公平的比较这些指标的性能，对结果进行归一化处理。首先，每个指标的振幅除以所有振幅的总和，使它们在[0,1]内；然后，将第一步结果的平方根作为最终的归一化振幅，结果如图1所示。

为了扩大复合影响指数的应用范围，以下证明该指数对其余周期性冲击依然有效。现构建一个包含周期冲击、随机冲击、周期谐波和高斯白噪声的模拟信号，分别计算KES、KSES、EK和EE等指标，并进行归一化处理，结果如图2所示。

2.2 自适应变分模态分解

对复杂环境下采集的故障信号进行VMD分解时，需要确定模态数K、惩罚因子α、噪声容限τ和收敛误差ε等参数。噪声容忍度和收敛误差对分解结果的影响较小，一般采用默认值。模态数K和惩罚因子α需要通过蝙蝠算法对全局进行有效的搜索。在选取目标函数时，应最大程度上考虑分解后的模态分量的故障信息。因此，本文提出最小平均复合影响指数(minimum average composite impact index, MACII)作为目标函数，可表示为

(18)

针对变分模态分解参数取值难以确定，严重依赖经验选取而导致解析力低的问题。本文提出以最小平均复合影响指数为目标函数，将蝙蝠算法对VMD参数进行全局寻优，自适应地实现信号的频域剖分及各分量的有效分离，具体流程如图3所示。

自适应变分模态分解实施步骤为[15]：

步骤1初始化蝙蝠算法参数，设置VMD参数K和α的范围，确定适应度函数；

步骤2将信号进行变分模态分解，根据式(18)计算信号的平均复合影响指数；

步骤3将平均复合影响指数作为适应度函数进行全局搜索，寻找最小值；

步骤4根据式子更新蝙蝠个体的位置和速度；

步骤5重复步骤2～步骤4，直至确定最小复合影响指数或达到所设定的迭代循环次数，输出最佳蝙蝠个体(K，α)；

步骤6使用具有最佳参数的VMD分解原始信号。

3 仿真分析

为了验证该方法的可行性，建立滚动轴承复合故障特征模型，进行自适应变分模态分解，对分解后的信号进行包络解调分析，从包络谱中识别出复合故障特征。

考虑到机械系统工况复杂，采集到的振动信号信噪比低，存在各种干扰成分，特别是由轴旋转或齿轮啮合引起的随机冲击，严重阻碍了故障脉冲的提取。为了模拟复杂工况下采集到的真实振动信号，在振动信号中加入背景噪声、随机冲击和周期性谐波。以滚动轴承为例，其复合故障模拟信号由以下5部分组成[16]

(19)

在式(19)中：第一项和第二项分别代表轴承外圈故障和内圈故障造成的周期性冲击；第三项表示由背景和电磁干扰等产生的随机冲击；第四项表示周期谐波，最后一项表示高斯白噪声。在该振动模型中，A，B，R和P是不同项的振幅，T是两个相邻脉冲之间的时间间隔，τ是由于滚珠随机滑动效应而引起的时间滞后，它通常占T的1%～2%，用于模式振幅调制。fh表示谐波干扰的频率，φh为相位。s(t)为脉冲响应的函数，其表达式为

s(t)=e-βntsin(2πfnt+φ)

(20)

式中：βn和fn分别为阻尼系数和共振频率；φ为相位。

通常，轴承外圈故障冲击的振幅(式(19)中的Ai)为常数，轴承内圈故障冲击振幅计算如下

B=1+Bbsin(2πfrt+φb)

(21)

本文构建了含有高斯噪声的的复合故障模拟信号，如图4所示。外圈故障、内圈故障、随机冲击、周期谐波、高斯噪声和混合信号的时域波形如图4(a)～图4(f)所示。关键参数值如下：采样速率为12 kHz，采样数量为7 200。轴旋转频率fr为30 Hz，内外圈故障特征频率fi和fo分别为110 Hz和80 Hz(Ta=1/fo，Tb=1/fi)。外圈故障、内圈故障和随机冲击共振频率分别为3 000 Hz，4 000 Hz，4 500 Hz。其阻尼系数分别为1 000，1 200，1 500。Ai和Bb分别为1.2和1.4，τi和τj分别为0.01Ta和0.01Tb。Rm和Tr由MATLAB中randn函数随机选择产生，Rm的最大值为3。周期谐波部分写为Psin(2πf1t)cos(2πfrt)，其中振幅P为0.01，频率f1为10 Hz，高斯白噪声的信噪比为5 dB。

如图5所示，由于背景噪声的干扰，频谱图中谱线较为混乱，在3 000 Hz和4 000 Hz出现共振频带，无法识别出轴承内圈和外圈的故障频率，得不到所需的有效信息。包络谱中虽然能够找到内、外圈故障特征频率，但是这些脉冲相当随意，没有呈现规律变化，干扰谱线较多，加大复合故障分离难度。

为了能够精确识别出轴承复合故障的类型，通过蝙蝠算法对变分模态分解进行参数寻优，设置蝙蝠迭代次数为10，种群数量为20。平均复合影响指数随蝙蝠更新代数变化趋势如图6所示，从图6中可以看出，随着迭代次数增加，适应度值逐渐减小，这说明采用优化后的参数后，分解结果逐渐接近最优值，最后，通过不变的适应度表明，找到了最优参数。当迭代次数为3时，适应度值最小，此时模态数量和惩罚因子分别为(5,100)。选取最优参数对仿真信号进行变分模态分解，结果如图7所示，从频谱图中，频率从低到高分布，没有出现频率混叠的现象。从包络谱中可以看出，故障脉冲信号分别位于3模式和4模式，其他模式相对于故障诊断是冗余的。从模式3中，可以清楚地找到外圈故障的频率以及多个倍频。模式4中同样可以找到内圈故障的频率及多个倍频，在其倍频两侧均匀分布着谱线，与倍频间隔都是fr，可见转频信号对外圈故障信号起着调制作用。为了消除转频对诊断结果的影响，需要结合其他降噪方法进行进一步的研究。

4 案例研究

采用凯斯西储大学的滚动轴承试验数据验证本文所提方法的有效性，试验台如图8所示。试验台由驱动电机、扭矩传感器、控制器和故障轴承等组成，电机驱动端安装有SKF6205型号的深沟球轴承，加速度传感器安装在轴承的正上方[17]。测试轴承通过电火花加工的方式在其内圈，外圈加工直径为0.177 8 mm，深度为0.279 4 mm的单点模拟故障，试验条件如表1所示。在本试验中，采样频率为12 kHz，采样时间为1 s，采用Miao等的试验方法，将轴承内圈故障信号和外圈故障信号混合得到复合故障信号。考虑到现场噪声干扰的随机性，在试验信号中加入了5 dB的高斯白噪声。

根据表1中的参数，内圈、外圈故障特征频率由式(22)和式(23)计算得到

表1 试验条件Tab.1 Test conditions

fbpfi=5.415 2×fr=162.185 Hz

(22)

fbpfo=3.584 8×fr=107.365 Hz

(23)

如图9所示，周期脉冲信号杂乱无章难以识别，埋没在噪声之中。图10为试验信号的包络谱图，包络谱中虽然能找到内外圈故障频率，由于受到背景噪声的干扰，周围存在众多谱线，影响轴承故障类型的诊断。因此，需要从噪声信号中进一步提取故障特征。

首先需要对VMD的模态数量和二次惩罚因子进行优化。设置种群数量为20，迭代次数为10。图11显示了搜索结果，随着迭代次数的增加，平均复合影响指数逐渐减小。当迭代次数为4时，得到最优参数，此时K和α分别为(5,429)。采用最佳参数对信号进行变分模态分解，分解信号如图12所示，对分解后的信号进行包络解调分析，如图12(b)所示，频谱图中没有出现模态混合和复制的现象。从图12(c)中可以看出，故障脉冲信号分别位于模式3和模式4，其他模式相对于故障诊断是冗余的。从模式3中，可以清楚地找到外圈故障的频率、2倍频以及3倍频，模式4中同样可以找到内圈故障的频率及2倍频，同时混杂着转频及其倍频，需要进行一步的研究，去除转频的影响。