基于FVMD的滚动轴承故障特征提取方法

张 爽， 王晓东， 李 祥， 杨创艳

(1.昆明理工大学 信息工程与自动化学院，昆明 650500； 2.昆明理工大学 云南省人工智能重点实验室，昆明 650500)

滚动轴承作为机械设备的关键零部件，常常工作在高温、高压和复杂的力学环境中，极易发生故障[1]。轴承一旦发生故障，就会导致设备停机，生产线难以正常运行，造成重大的经济损失，甚至导致灾难性的人员伤亡[2]。因此，提取滚动轴承的故障特征对及时有效判定其工作状况和保证生产安全具有重大意义。

振动信号分析作为滚动轴承故障诊断中最受欢迎的技术手段。其中，经验模态分解(empirical mode decomposition,EMD)是最具代表性的方法之一，因其在处理非线性、非平稳信号时有着良好的效果，并在滚动轴承故障诊断领域中得到了广泛的应用[3-4]。但同时存在模态混叠，端点效应等问题，严重影响其分解效果[5]。为此，有学者提出了局部均值分解(local mean decomposition,LMD)方法，可很好地改善EMD所遗留的模态混叠，端点效应等问题，保证故障特征提取正确性。但是，LMD本质上仍然是递归模态分解，容易受采样频率的影响，难以保证分解精度[6]。为解决EMD和LMD递归模态分解的问题，Dragomirestkiy等[7]提出变分模态分解(variational mode decomposition,VMD)算法，该算法对信号进行分解时采用的是一种非递归的处理技术，弥补了循环递归筛选的缺陷，从而有效地摆脱LMD与EMD存在的端点效应、模态混叠等束缚[8]。

VMD能克服端点效应和模态混淆的弊端，但涉及到3个参数：模态数K、惩罚因子α、模态初始中心频率ω的预设置。模态初始中心频率ω在传统VMD中给出了3种初始化方式，其中ω等于零和ω在被分解信号的频谱上呈线性分布的方式具有较高的计算效率，另一种ω随机取值的方式效率较低。Bi等[9]将信号频谱中幅值最大点对应的频率作为模态初始中心频率ω，并把模态数K设为1进行分解，有效的减少VMD的模态中心频率迭代次数。另一参数K太小会导致分解结果不完全，故障信息被遗漏，而模态数K太大会导致过分解[10]。惩罚因子α越小，分解后的本征模态函数(intrinsic mode function,IMF)分量衰减越慢、带宽越大，相反，分解后IMF分量衰减越快、带宽越小[11]。Ma等[12]通过尺度空间频谱分割的方法来确定模态数K的值，同时把惩罚因子α默认为2 000，实现了参数K的自适应选择。刘长良等[13]将模态数K设置为2、惩罚因子α设为2 000，并计算VMD分解后IMF的中心频率，若出现了中心频率相近的模态则认为出现过分解，模态数K=K-1,否则K=K+1再次计算IMF的中心频率。Jiang等[14]把模态数K设为1、惩罚因子α设为2 500进行信号的分解，并把残余信号依次分解，计算分解出模态的峭度值，当峭度值开始下降时停止分解，对峭度值最大的模态进行分析。

上述大部分研究把惩罚因子α设为2 000，对模态数K进行优化，但没有考虑到惩罚因子α变化带来的影响。目前，大部分研究采用优化算法对VMD的参数K和α同时进行选择：张俊等[15]以包络谱峰值因子作为优化指标，采用粒子群优化算法寻找最优的模态数K和惩罚因子α；Zhang等[16]把权重峭度作为优化指标，通过蚱蜢优化算法搜寻模态数K和惩罚因子α。利用优化算法和VMD相结合，能够对参数进行有效的设置，但使得整个流程的计算效率降低、忽略了VMD算法的性质。综合上述存在的问题，本文从对模态初始中心频率ω的定位入手，提出了一种快速变分模态分解(fast VMD,FVMD)算法。算法对不同ω的模态计算出对应的α值，能同时对K和α值进行优化。相比于原VMD算法的模态初始中心频率ω定位方式，还能够增加其计算效率。

1 基本理论

1.1 VMD

(1)

式中：δ(t)为冲击函数；f(t)为输入信号； j为虚数单位；∂t为对参数t求偏导。VMD算法可以将输入信号分解为K个以ωk为中心频率的IMF分量uk(t)。VMD算法可以分为两个主要部分：第一个部分是变分问题的构造；第二部分是变分问题的求解。约束变分问题的构造可以用式(1)来进行描述。

然后，引用拉格朗日乘子λ和二次惩罚因子α，将约束变分问题转化为无约束变分问题进行求解，其增广拉格朗日如式(2)所示。

(2)

采用乘子交替方向法求解增广拉格朗日表达式，模态uk(t)及模态中心频率ωk的更新步骤如下所示。

步骤2n=n+1。

(3)

步骤4对拉格朗日算子λ进行了更新如式(4)。

(4)

τ为噪声容限参数，为了获得更好的去噪效果，可设τ=0。

重复步骤2～步骤4，直到满足式(5)则停止迭代。

(5)

式中，ε值通常设置为1×10-6，详细内容请参考Liu等的研究。

1.2 阈值去噪

阈值去噪是小波阈值去噪的主要部分，进行小波阈值去噪时，首先需要设置一个临界阈值λ。当小波系数小于λ，认为该系数主要由噪声引起，去除这部分系数；若小波系数大于λ，则认为此系数主要是由信号引起，予以保留。因此，将此想法应用到频谱趋势的计算中，消除噪声引起的频谱趋势波动。

此处采用硬阈值去噪的方法对频谱趋势进行去噪处理

(6)

阈值的计算公式如式(7)[17]所示

(7)

1.3 频谱趋势分割

快速经验小波变换(fast empirical mode decomposition,FEWT)[18]通过快速傅里叶变换(fast Fourier transform,FFT)和快速傅里叶反变换(inverse FFT,iFFT)计算信号的频谱趋势，对频谱趋势使用阈值去噪进行优化，取优化后频谱趋势的极小值点作为分割边界。本文引用频谱趋势分割的方法，来确定VMD算法中模态数K。同时，VMD优先对频谱中能量大的调频调幅部分进行分解，将故障信号中调频调幅部分的中心频率作为VMD模态初始中心频率能够降低模态迭代次数。而频谱趋势能够凸显此部分，因此，以频谱趋势极大值对应的频率作为模态分量的初始中心频率。其方法流程如图1所示，结合图1和式(8)，方法实现过程阐述如下。

步骤1通过FFT计算得到信号的频谱，如图1(a)所示。对信号的频谱再次使用FFT得到信号频谱的核函数，如图1(b)所示。

步骤2取核函数的前30个点(图1(b)虚线左边的部分)，利用iFFT计算得到频谱趋势(图1(c)频谱上方曲线)，频谱趋势能够刻画故障信号在频谱上的整体走向，并根据其组成成分能量的大小，频谱趋势的高低也随之变化。并把趋势的极小值作为分割的边界，如图1(c)虚线所示。

步骤3利用硬阈值去噪的方法来优化频谱趋势，被噪声干扰的频谱趋势部分呈直线，然后以频谱趋势极小值点和直线趋势的端点作为新分割边界，如图1(d)虚线所示。将图1(d)中的两段直线趋势所在频段视为无效频段(a段和b段)，其余带有趋势极大值点的4段作为有效频段，即模态数K=4。

步骤4通过优化的频谱趋势分割图，可以容易的找到频谱趋势极大值点，将极大值点对应的横坐标频率作为模态初始中心频率即ω=[ω1ω2ω3ω4]，如图1(d)所示。

图1 频谱趋势分割

(8)

式(8)由谐波和周期脉冲组成，同时加入高斯白噪声来模拟实际信号中的噪声干扰，谐波信号y1、y2、y3的频率分别为40、400、900，旋转频率fr为25 Hz，故障频率fm为70 Hz，衰减系因子ζ为1 000，共振频率fn为2 000 Hz，n(t)为零均值1方差的高斯白噪声，采样频率为10 kHz，采样点数N为4 096。

1.4 惩罚因子的计算

根据滚动轴承故障振动信号的频谱分布特征，中低频区域主要由旋转频率及其相关特征频率(如轴承故障特征频率和齿轮啮合频率等)的谐波组成，而故障冲击和噪声干扰大多位于高频区域[19]。同时谐波信号具有时域持续时间长、频域上较为紧凑的特点，而冲击信号具有时域短、频域宽的特点。

因此，以各模态分量的中心频率作为确定相应惩罚因子的依据。通过频谱趋势分割得出的模态初始中心频率较小，说明模态提取的部分主要是谐波，惩罚因子应该大；相反，说明模态提取的部分主要是故障冲击部分，此时的惩罚因子应该选择较小的值。因此，可以建立式(9)的映射关系

(9)

1.5 有效权重峭度指标

从分解的K个IMF分量中筛选含丰富故障信息的分量是另一个需要解决的核心任务。峭度指标通常作为冲击测量指标用于监测机械损伤[20]。但是，峭度指标仅仅考虑冲击信号的分布密度，而忽略振幅较大、分布不均匀的冲击分量。相关系数可以表征两个信号的相似性，能够考虑到冲击分量的不均匀性，但在冲击信号检测中容易受到噪声的影响。

结合两个指标的特点，引入有效权重峭度指标来筛选分量，该指标既考虑到脉冲分量的不均匀性，又有效地滤除了噪声[21]。该指标可以根据模态分量自身的情况，从中选出含冲击信息多、噪声少的有效分量，计算公式如式(10)所示。

(10)

式中，Ki和Ci、KEW,i分别为第i个模态分量的峭度值和相关系数、有效权重峭度指标，i,j=[1,2,3，…，K]。将KEW(k)>0的模态判断为有效模态分量，其余为含噪模态分量。

2 FVMD的故障特征提取方法

根据第1章的理论分析，本文提出了一种FVMD的滚动轴承故障特征提取方法。方法的流程如图2所示，结合图2，方法实现过程阐述如下。

图2 FVMD故障特征提取流程

步骤1收集滚动轴承故障信号，用频谱趋势分割法进行预处理。

步骤2通过频谱趋势分割，确定VMD算法参数ω、K、α，并完成信号的自适应分解。把有效频段中的极大值所对应频率ωk作为模态初始中心频率，ω=[ω1ω2ω3ω4…ωk],1≤k≤K；通过式(9)计算出每个初始中心频率ωk对应的惩罚因子αk，α=[α1α2α3α4…αk]，1≤k≤K；改变模态初始中心频率和惩罚因子初始化设置方式，利用优化后的K、ω和α更新传统VMD参数，并执行信号的VMD分解。

步骤3采用有效权重峭度准则，选取有效模态分量进行重构，并对重构信号进行包络解调完成故障特征提取。计算分解后各模态的峭度值和相关系数；根据式(10)计算各模态的有效权重峭度指标，将有效权重峭度指标大于零的模态分量进行重构；包络解调获取故障特征。

3 仿真试验对比分析

为了验证所提出方法的有效性，对式(8)组成的仿真信号进行分析，其波形如图3所示。

图3 仿真信号

根据图3(c)的频谱趋势分割图可以快速获得FVMD的ωF=[0,460,937,2 031]和K=4,并根据式(9)计算出对应的惩罚因子αF=[2 500,1 191,871,482]，从而分解得到模态，如图4(a)所示。并与传统VMD算法中的两种模态初始中心频率设置方式进行对比，论证FVMD方法动态设置ω和K的有效性。①是ω=0，惩罚因子α设置为2 000、模态数K=4,如图4(b)所示；②是ω在频谱上呈线性分布，惩罚因子α设置为2 000、模态数K=4，如图4(c)所示。此外，图4(d)还给出ω=ωF和惩罚因子α=2 000、K=4的分解结果，论证了惩罚因子α动态设置的必要性。

同时，图4(b)与图4(d)分解出4个故障相关模态，而图4(c)未分解出第二个谐波，验证了模态初始中心频率的合理设置的必要性。

图4 4种VMD分解方式的IMF分量

为了选取有效分量进行后续的频谱分析，分别计算图4(a)～图4(d)中IMF分量的EWK(effective weight kurtosis)，如表1所示。以及对应IMF分量的中心频率迭代次数，如图5所示。根据有效权重峭度准则，选取表1中EWK>0的分量进行重构，并计算包络解调频谱，如图6所示。

表1 4种VMD分解方式下IMF分量的EWK

图6 4种VMD分解方式的重构信号包络解调

分析图5～6可知：①图5(a)FVMD根据频谱趋势分割设置的ωF，相比于图5(b)～图5(c)传统VMD中ω的设置方式，模态中心频率的迭代次数最少，进而提高了VMD的计算效率。在优化ω的基础上，同时对惩罚因子进行优化，相比α=2 000进行分解的方式(见图5(d))增加了迭代次数,但是FVMD方法迭代次数远低于传统VMD方法。②图6(a)所提方法成功提取出故障的基频和2倍～6倍频。而图6(b)～图6(d)中令α=2 000进行VMD分解的方法只提取到故障基频和2倍～4倍频。说明用此方法对α进行优化，有利于取到更丰富的故障特征。③同时计算4种方法重构信号的峭度值，可以看出图6(a)FVMD的峭度值最大，说明了FVMD的有效性。

4 试验数据对比分析

为了验证基于FVMD的滚动轴承故障特征提取方法在实际工程中的有效性和优越性，以下试验使用CWRU(Case Western Reserve University)数据集和NASA(National Aeronautics and Space Administration)数据集并基于表2中的方法进行分析。

表2 试验方法

4.1 CWRU滚动轴承振动数据分析

CWRU滚动轴承故障模拟试验平台如图7所示，具体参数如下：轴承型号为SKF6205，转速为1 750 r/min，采样频率为12 kHz，采样点数为2 048。选取滚动轴承内圈、外圈故障振动信号进行试验分析，滚动轴承外圈理论故障频率BPFO=107.305 Hz、内圈故障频率BPFI=162.185 Hz。

图7 CWRU轴承数据试验平台

4.1.1 轴承外圈故障信号试验对比分析

轴承外圈故障信号与频谱趋势分割如图8所示，可计算得到模态数K=3，初始中心频率ωF=[375，2 813，3 375]、惩罚因子αF=[1 616，489，372]。将其参数代入VMD中分解得到系列IMF分量，如图9(a)所示。同时给出第3章所述的其他3种方法的IMF分量，如图9所示。

图8 外圈故障信号

图9 4种VMD分解方式的IMF分量

从频域上看，4种方法都分解出3个故障相关模态，与频谱趋势分割中的有效频段相对应。但仅通过图9对比分析无法直观看出区别，需要进一步计算图9(a)～图9(d)中各IMF分量的EWK来对比分析，计算值如表3所示。同时，给出各IMF分量的中心频率迭代次数，如图10(a)～图10(d)所示。由表3可知，将EWK>0的IMF分量视为有效分量予以保留，并进一步对有效分量重构后进行包络解调，如图11(a)～图11(d)所示。

表3 4种VMD分解方式下IMF分量的EWK

对图10～11进行对比分析可知：①图10(a)FVMD方法各模态的中心频率迭代次数最少为14，均小于图10(b)～图10(c)传统VMD方法的两种初始化方式，故本文所提方法能有效减少VMD算法的迭代次数。②由图11(a)～图11(d)可知，4种方法都存在与外圈故障频率理论值(107.305 Hz)非常接近的频率105.5 Hz，还存2倍～7倍频。但是图11(a)FVMD方法提取到的故障基频和倍频幅值都大于其他3种方法。③同时，从图11中可知，所提方法的峭度值最大，进一步验证了动态设置α不仅能保留更多的故障信息，而且能突出振动信号的故障冲击分量。

图10 4种VMD分解方式的IMF分量中心频率迭代次数

图11 4种VMD分解方式的重构信号包络解调

4.1.2 轴承内圈故障信号试验对比分析

按照第3章的分析思路，应用频谱趋势分割计算得到内圈故障振动信号的模态数K=4、模态初始中心频率ωF=[562，1 125，2 625，3 563]。根据式(9)计算惩罚因子αF=[1 420，1 045，533，338]，将其参数K、ωF和αF代入进行内圈振动信号的VMD分解，并给出其他3种分解方式的IMF分量。计算分解后各IMF分量的EWK值，将大于零的分量进行重构。并计算其包络解调，如图12所示。

图12 4种VMD分解方式的重构信号包络解调

由图12可知：①通过4种方式进行信号分解，都提取到与内圈故障频率理论值(162.185 Hz)非常接近的频率164.1 Hz，还存2倍～7倍频。但是，图12(a)FVMD方法提取到的故障特征频率幅值最大、故障特征最明显。②进一步计算4种方法重构信号的峭度指标峭度值，图12(a)所提方法的峭度值最大，具有最丰富的故障特征。图11与图12对比可知，外圈与内圈故障振动信号试验取得一致的结论，进一步论证所提方法(FVMD)的有效性和可行性。

4.2 NASA滚动轴承振动数据对比分析

NASA轴承全寿命周期试验平台如图13所示。采样频率为20 kHz，驱动电机转速为2 000 r/min，径向负荷2 721.554 22 kg，共采集984组数据,每组长度为20 480。试验采用外圈第874组数据进行分析。

图13 NASA轴承数据试验平台

同理，采用FVMD方法对故障振动信号进行处理，获得信号模态数K=5、对应的初始中心频率ωF=[0，937，1 563，3 438，4 375]、惩罚因子αF=[5 000，2 365，1 913，1 139，888]，4种方法的重构信号包络频谱如图14所示。

图14 4种VMD分解方式的重构信号包络解调

结合图14进行分析可知：①图14(a)～图14(d)4种分解方法都能提取到故障的基频和倍频。但是，图14(a)FVMD方法提取到的基频和倍频幅值均大于图14(b)～图14(d)(α=2 000、ω=0或ω=呈线性或ω=ωF)的方法。②进一步计算重构信号的峭度指标峭度值，对比分析可知所提方法(FVMD)取得最大的峭度值，进一步验证FVMD方法的优越性。

5 结 论

本文通过仿真信号和CWRU、NASA滚动轴承数据的试验及对比分析，验证了基于FVMD的滚动轴承故障特征提取方法的有效性，可以得到以下结论：

(1)针对传统VMD方法需要预先给定分解模态数K和惩罚因子α的局限性，本文结合频谱趋势自适应分割方法，实现了VMD分解模态数K、惩罚因子α和模态初始中心频率ω的自适应选择。

(2)本文通过对模态初始中心频率ω和惩罚因子α进行动态优化降低了模态中心频率的迭代次数，提高了计算效率；同时，所提方法比传统VMD方法提取到了更丰富故障特征。

(3)针对如何从IMF分量中筛选有效分量的问题。结合峭度指标和相关系数指标建立了有效权重峭度指标筛选规则，实现了有效分量的选取。

(4)通过仿真信号和CWRU、NASA滚动轴承数据的对比分析，有效论证了FVMD滚动轴承故障特征提取方法的可行性和有效性。