基于自适应MOMEDA与VMD的滚动轴承早期故障特征提取

刘 岩，伍 星，刘 韬，陈 庆

(昆明理工大学 云南省高校振动与噪声重点实验室，昆明 650500)

轴承在机械设备中的应用十分广泛，与此同时轴承也是设备中较易损坏的零件。一旦轴承失效会发生机械故障并可能出现严重的事故，因此对轴承进行早期故障特征提取具有重要意义。早期故障信号冲击成分微弱，?但实际工况下又有广泛的背景噪声存在，并且轴承信号在轴承座中传递会被衰减，传感器得到的振动信号是经过系统衰减后的信号，所以采集到的信号信噪比低，特征提取困难。因此，早期故障的特征提取一直是轴承故障诊断领域的热点。

针对此问题，近年来学者们进行了深入的研究。例如王宏超等[1]对轴承微弱故障信号利用MED降噪，并将降噪之后的信号使用匹配追踪(Matching Pursuit，MP)的方法进行了稀疏重构，对重构信号进行了特征提取，取得了良好的效果。唐贵基等[2]提出使用最大相关峭度解卷积(Maximum Correlation Kurtosis Deconvolution，MCKD)降噪并结合1.5维谱的滚动轴承早期故障提取方法，并取得了一定的效果。刘韬等[3]对轴承微弱故障使用了短时傅里叶变换计算了各个频带的频带熵，并依据“最小频带熵处存在轴承的共振频率”设计了带通滤波器，提取滤波之后的信号的故障特征，取得了较好的结果。苏文胜等[4]使用EMD分解了滚动轴承早期故障信号，并依照峭度和相关系数将分量信号叠加获得了重构信号。对重构信号计算了谱峭度，构建了最优带通滤波器进行滤波，并提取滤波之后的信号的故障特征，取得了一定的效果。毕果等[5]利用循环平稳理论对轴承信号进行了分析，并针对轴承微弱故障利用二阶循环理论对信号的谱相关密度函数做切片分析取得了较好的效果。Huang等[6]使用盲分离理论并结合谱峭度对滚动轴承早期故障进行了诊断，取得了一定的效果。唐先广等[7]针对激振器干扰情况下的滚动轴承故障信号使用了ICA对EMD分出的IMF进行降噪，取得了较好的效果。Rubini等[8]提出使用小波变换与包络分析结合的方法对滚动轴承早期故障进行特征提取，取得了一定的效果。

最近，Mc Donald等[9]在优化最小熵解卷(Optimal Minimum Entropy Deconvolution，OMED)方法的基础上提出了一种新的MED算法——多点优化最小熵解卷积(Multipoint Optimal Minimum Entropy Deconvolution，MOMED)。该方法以D范数最大为目标，并根据预先估计的故障频率构建了目标向量(Ttarget Vector)，不使用迭代的方法求出了解卷积的最优滤波器的系数。本文尝试将VMD与MOMED方法结合，针对轴承早期故障，利用MOMED降噪之后的信号进行VMD分解，并将分解之后的信号依照全频带谱峭度大小进行重构，再对重构信号进行故障特征提取。

1 基础理论

1.1 MED

MED目的是使解卷积之后的结果中部分脉冲被突出，该特性非常适用于滚动轴承低信噪比故障信号的降噪，降噪后的轴承冲击被突出。假设轴承发生故障时的信号如式(1)所示

xn=h*y+e

(1)

式中:xn为发生故障时采集到的信号，y为系统卷积前的故障冲击信号，h为系统冲击响应函数，e为噪声。在不考虑噪声影响下，最小熵解卷积问题等价于寻找找一个长度为L逆滤波器f使得y能被滤波器f还原出来，也就是如式(2)所示

y=f*xn

(2)

为了达到这个目的，MED算法需要利用信号的熵求极值来构建滤波器，并对原信号滤波从而达到解卷积的作用。但是，由于信号的熵求极值是采用迭代法，因此最终只能求出局部最优解而无法求出全局最优解。严格上说这种方法只是找到了一个“相对较好”的滤波器而不一定是“最好的”滤波器。

此外，如文献[9]中所述MED算法在初始采样点处还依赖于一个假设即：初始采样点之前的信号值为0，如式(3)所示

xn=0,n<1

(3)

当条件不满足时可能会导致滤波后的时域信号中出现虚假冲击，这些冲击的幅值会干扰后续的故障特征提取和辨识。

2.2 MOMEDA算法

为解决MED初始采样点处的假设与实际工况不符的问题，有学者提出了改进卷积定义来避免端点假设的方法[9]。传统的卷积定义如式(4)所示

(4)

yn为滤波后的实际信号，X0定义如式(5)所示

(5)

改进卷积定义中的X0如式(6)所示

(6)

将式(6)代入式(4)，可得到改进最小熵解卷积(Minimum Entropy Deconvolution Adjusted,MEDA)法。

为了解决不使用迭代求出“最优”MED滤波器的问题，Cabrelli等[10]提出了不使用熵而使用D-范数作为指标来求出“最优”的MED滤波器，也就是优化最小熵解卷积OMED(Optimal Minimum Entropy Deconvolution)方法。D-范数定义如式(7)所示

(7)

式(7)的优点就是对滤波器参数求极值就可得到D-范数最大的解析式，避免了迭代而直接求出了“最优”MED滤波器。但OMED滤波后的信号会有多个输出信号，每个输出中只含有一个冲击脉冲，这就限制了OMED的实际工程应用。

针对OMED的不足，Mc Donald提出多点优化最小熵解卷积MOMED，该方法定义了一个t由0和1组成目标向量(Target Vector)，其中1代表冲击数。利用目标向量和D-范数重新定义了新的指标多点D-范数(MultiD-Norm)如式(8)所示：

(8)

将最小熵解卷问题变为：

(9)

通过对式(9)求解，就可以得到MOMED方法计算出的滤波信号：

(10)

当式(10)中X0取式(6)时MOMED变为改进多点优化最小熵解卷(Multipoint Optimal Minimum Entropy Deconvolution Adjusted,MOMEDA)。

1.3 基于VMD信号重构

Dragomiretski等提出的分解方法变分模态分解(Variational Mode Decomposition，VMD)将信号分解转化非递归、变分模态分解模式，其实质是多个自适应的维纳滤波器组，因此比EMD和LMD等方法具有坚实的理论基础[12]。其次，VMD使用交替方向乘子法，将各模态解调到各自中心频率附近的频带上，并且保证各频带长度最小，避免了各模态在频域上的混叠现象。

设输入信号Y的每阶模态函数为ui共有k个模态，VMD就是在各阶模态相加等于输入信号为约束条件下模态的长度求和最小的优化问题。设一个模态的预估中心频率为ωi，于是VMD问题变为：

(11)

使用惩罚函数法求解上式。设惩罚因子为C，拉格朗日乘法算子为λ(t)。构造拉格朗日函数为：

y,λ(t)>

(12)

(13)

(14)

(15)

在进行了VMD分解后需要使用一些标准选择几个“较好”的IMF分量相加获得重构信号。常用的标准主要有：峭度准则和相关系数准则。

相关系数准则是通过计算各IMF与原信号之间的相关系数来判断各IMF与原信号之间的相关性，并剔除与原信号相关性较小IMF。这种方法主要用于EMD分解时剔除伪分量。但VMD在计算IMF时已经保证了各IMF频带不混叠，因此这种方法对VMD不适用。

峭度衡量了波形的尖峰程度，如果一个IMF分量中峭度越大，则这个分量可能含有更多的冲击成分。选择峭度较高的分量进行重构，故障特征就会越明显。

峭度对于噪声干扰较小的信号有良好的效果，但作为一个全局指标在高噪声干扰的情况下不适合反映特点分量的变化情况。早期故障的冲击信号常常有极低的信噪比，虽然经过了MOMEDA降噪但仍然存在相当的噪声干扰。为了克服上面的缺点，本文使用了谱峭度作为指标来重构信号。

谱峭度(Spectral Kurtosis,SK)最早由Dwyer提出[13]。它是归一化的四阶累积量，它可以看做是一个频带内信号概率分布函数的极大值。由于谱峭度是高阶统计量，因此它具有对噪声的“盲性”以及对瞬态成分(冲击)的敏感性。由于严格意义上的谱峭度难以计算，通常需要进行估计。本文采用了Antoniz在文献[14]中提出的估计方法，对于第i个分量ui取其包络为ci，分量的谱峭度定义为：

(16)

由于本文只将谱峭度作为一个重构标准，并不选取共振解调频带，因此只计算每个分量在全频带的谱峭度，并选取谱峭度最大的两个分量进行重构。

1.4 自适应MOMEDA与VMD的滚动轴承早期故障特征提取

MOMEDA算法虽然有诸多优点，但其最终的效果还和滤波器长度宽度L的设定有关。本文基于峭度指标使用进退法自适应寻找最优滤波器长度宽度，但直接进退法导致的计算量过大，先在确定的寻优范围[a,b]内使用峭度指标寻找最佳起始点，然后使用进退法计算最佳滤波器长度。具体流程如图1所示。

其中自适应算法步骤如表1所示。

通过进退法以峭度为指标搜索出了较好的滤波器长度从而使MOMEDA方法有了自适应性。这种方法适合低信噪比具有周期性的滚动轴承早期故障信号的降噪。

2 滚动轴承微弱故障仿真验证

使用文献[1]所示的滚动轴承故障模型对轴承外圈点蚀故障进行仿真。设置采样频率为fs=51 200 Hz，外圈故障频率为fo=100 Hz，系统固有频率为fn=3 000 Hz，并在信号加入-19 db的白噪声模拟早期故障。仿真信号波形与包络谱，如图2所示。

图1 算法流程

表1 进退法选最优滤波器算法

观察图2(a)发现此时冲击成分极为微弱，且观察图2(b)仿真信号的包络谱发现外圈故障频率fo=100 Hz被噪声干扰，几乎被淹没。仿真信号使用MOMEDA进行滤波得到滤波后信号的包络谱如图3(b)所示，对比图3(b)与图2(b)，发现此时故障特征频率被提取出来，但仍受到噪声干扰。

(a) 仿真信号时域波形

(b) 仿真信号包络谱

(a) 仿真信号MOMEDA滤波后时域波形

(b) 仿真信号MOMEDA

再对MOMEDA滤波后信号使用VMD分解重构，设置分解个数为K=4。根据实际情况本文选择α=200，分解出4阶IMF分量如图4所示。

利用谱峭度准则优选分量得到重构信号及其包络谱，如图5所示。

图5(b)中出现了非常明显的外圈故障频率及其倍频，对比图3(b)和图5(b)可以发现VMD重构信号的效果比仅使用MOMEDA滤波后的效果更好。

图4 VMD法分解出的4阶IMF时域波形图

(a) 重构信号时域波形

(b) 重构信号包络谱

3 实验验证

实验使用ABLT-1A全寿命试验台早期故障数据和加工轴承外圈故障实验数据并添加白噪声模拟轴承早期外圈故障，对本文提出方法进行验证，并使用全寿命数据将本文提出方法与其它方法进行了比较。

3.1 全寿命轴承数据采集说明

本文所用的早期故障信号来自于轴承全寿命实验台ABLT-1A，ABLT-1A实验台系统如图6所示。

实验时转速3 000 r/min，转频fr=50 Hz，加载向径载荷FR=11.8 N，实验轴承型号为GB6205，轴承滚珠数为z=9，滚珠直径为d=7.5 mm，滚动体节圆直径为D=39 mm,接触角为β=0°。轴承故障特征频率计算公式如下

(17)

(18)

(19)

(20)

计算得内圈故障频率为fi=268.269 Hz，计算得外圈故障频率为fo=181.731 Hz，计算得保持架故障频率为fc=20.192 Hz，计算得滚子故障频率为fbs=125.192 Hz。

轴承全寿命实验数据采集系统利用NI 9234采集卡采集振动信号，采样频率为fs=51 200 Hz，每隔7分钟采集10 s的信号，振动传感器型号为DH112。整个实验持续约17天。通过NI 9234采集卡采集振动信号，采集程序基于NI LabVIEW编写。

图6 ABLT-1A实验台

由图7的轴承全寿命RMS和峭度变化可以看到，轴承的全寿命数据大致可分为4个阶段：正常、早期故障、严重退化和失效。轴承全寿命周期RMS变化如图7(a)所示。观察图7(a)发现在B时刻之前RMS没有明显的波动，在B时刻之后轴承的RMS值快速增加。因此B时刻之后轴承进入严重退化阶段，但从图7(a)中无法确定早期故障的时间段。

轴承全寿命周期峭度变化如图7(b)所示。观察图7(b)发现轴承峭度在A时刻之前没有明显的波动，在A～B之间有了微小的波动但波动幅度不大，判定段时间内的信号作为早期故障。观察图9(a)的正常信号(在15 050 min时采集)与9(b)，发现与正常信号相比早期故障幅值和冲击成分都略微增加。如图7(b)所示，B～C之间峭度波动增加，且B～C的峭度平均值相比A～B上升了约1.45；观察图9(c)发现此时冲击成分与冲击幅度都明显增加但轴承还可以工作，再结合B～C段峭度变化判定B～C段为严重退化阶段。观察图7(b)发现C点之后峭度陡然上升，再观察图9(d)发现冲击成分大幅增加，综合峭度与时域波形判定此阶段为失效阶段。

(a) 轴承全寿命RMS变化趋势

(b) 轴承全寿命峭度变化趋势

图8 轴承失效后的内圈损伤

(a) 正常信号时域

(b) 早期故障信号时域

(c) 严重退化信号时域

(d) 失效信号时域

根据上面分析结合A、B、C三点时间，确定轴承四个阶段时间为：正常(0～16 021 min)、早期故障(16 021～23 400 min)、严重退化(23 400～23 930 min)和失效(23 930～24 745 min)。

当实验结束后拆下轴承在其外圈处可以清楚的看到内圈的裂纹如图8所示。取第24 740 min时的数据(失效阶段)时域信号及其包络谱如图10所示。

从图10(a)中可以清楚的看到时域冲击及其调制现象；从图10(b)中可以清晰的看出轴承的转频、故障频率、故障频率的边带，通过上面的计算可以确定轴承此时出现了内圈故障，并结合图5发现实际情况相符判定内圈故障使轴承失效。

(a) 失效阶段时域波形

(b) 失效阶段包络谱

取第18 991 min时的数据(早期故障)时域信号及其包络谱如图11所示，观察图11(b)，在内圈故障频率与外圈故障频率处只有微弱的峰值，且内圈故障边带不可见，因此早期故障直接做包络谱无法确定其故障类型。

3.2 滤波器长度对MOMEDA方法滤波效果的影响

为了确定4.1节中较优的滤波器长度范围[a,b]，分别取全寿命实验四段早期故障的数据(四段数据采集时间为：18 550 min、18 921 min、19 264 min以及19 425 min)用MOMEDA方法在不同滤波器长度下滤波。每组数据使用MOMEDA方法从滤波器长度为10开始滤波，同时设置步长增量为10直至滤波器长度为1 000终止滤波并计算每个长度下的滤波效果，滤波效果评估指标为：滤波后信号的包络谱上故障频率处峰值与最大峰值的比值，比值越大则滤波效果越好，滤波器长度对MOMEDA滤波效果影响，如图12所示。

(a) 早期故障时域波形

(b) 早期故障包络谱

图12 滤波器长度对MOMEDA滤波效果影响

图12较粗的实线为四组数据在不同长度下滤波效果曲线的平均曲线。从图上可以看出在滤波器长度区间为[30,130]以及[600,700]内会有一个局部的较大值，寻优区间[a,b]应该是这两个区间中的一个。观察平均曲线发现区间[30,130]的平均效果比区间[600,700]更好。当信号长度相同时理论上使用MOMEDA方法的计算时间与长度成反比。在不同的滤波器长度下使用MOMEDA对18 550 min时长度为2 s的数据进行滤波，得出MOMEDA所需时间如图9所示。观察图13可以发现滤波器长度越宽所需计算时间越长，这与理论分析的结果一致。综合图12与图13，选定区间[30,130]为滤波器寻优范围。

3.3 全寿命轴承实验数据分析

3.3.1 自适应MOMEDA与其它方法及比较

取18 991 min的信号进行分析，并截取信号长度为1 s，采样频率为fs=51 200 Hz。作为比较，除自适应MOMEDA方法外，还对MED和MCKD以及小波包分解重构信号结果进行比较。原始信号的时域波形及包络谱如图7所示。通过4.1节分析，早期故障信号直接做包络谱无法确定其故障类型。

图13 滤波器长度对MOMEDA的计算时间影响

首先对这段信号使用N=10的Daubechies小波进行4层小波包分解。再选取重构分量信号峭度最大的作为重构信号，小波包重构信号及其包络谱如图14所示。观察图14(b)，发现外圈故障频率较为明显，但内圈故障频率非常微弱，且边带被噪声淹没。

(a) 小波包重构信号

(b) 小波包重构信号包络谱

使用MED法对这段信号进行滤波，并设置滤波器长度为41。滤波后效果如图15所示。图15(b)中可以观察到转频和内圈故障频率及其边带，且在外圈故障频率处出现峰值。

使用MCKD法对信号进行滤波，并设置滤波器长度为41。滤波后效果如图16所示。

观察图16(b)发现内圈故障频率相较于原始信号没有较大提升且看不出内圈故障的右边带，并出现了较大的外圈故障频率峰值。

(a) 早期故障使用MED滤波后信号

(b) 早期故障使用MED滤波后包络谱

(a) 早期故障使用MCKD滤波后信号

(b) 早期故障使用MCKD滤波后包络谱

使用自适应MOMEDA法对信号进行滤波，通过进退法在[30,130]区间上确定了MOMEDA的最优滤波器长度为41。滤波后效果如图12所示。

从图17(b)可以清楚的识别出转频fr、转频的二倍频以及内圈故障频率fi及其边带fi±fr此时内圈故障频率峰值为包络谱上的最大值，从图17(b)中可以诊断出此时轴承内圈存在故障，这与实际情况相符。三种方法都出现了外圈故障的特征，但最终实验结束后拆解轴承时发现外圈并没有明显的裂纹，可以判断轴承在早期故障阶段外圈也应该有轻微的故障，但外圈故障最终并未使轴承失效。

(a) 早期故障使用MOMEDA滤波后信号

(b) 早期故障使用MOMEDA滤波后包络谱

图18中列出了三种MED类方法滤波信号与小波包分解重构信号的包络谱中故障特征峰值与0～350 Hz谱总值之比，比值越大则对应故障特征约明显。从图18可以看出自适应MOMEDA方法对内圈故障及其边带提取的效果最好。MCKD与小波包分解重构虽对外圈故障提取较好但对内圈故障提取较差。综上面分析可以说明自适应MOMEDA方法优于其它两种方法。

图18 5种方法效果比较

3.3.2 滤波后基于VMD分解的重构信号特征提取

4.3.2节中的分析依然使用4.3.1节中的信号，所以自适应MOMEDA降噪后的时域图与包络谱依然如图17所示。在使用VMD分解时为了避免过多的计算，在保证精度的情况下将滤波后信号的采样频率下降为fs=25 600 Hz，并截取了0.5 s的信号进行分析。

VMD分解的参数选择对分解的效果有影响。VMD分解的分量个数设置需要适中，分量个数太少会导致每个分量有较宽的频带，这容易在分量上引入过多的冗余噪声。但分解的分量个数也不易过多，过多的分解会导致大量的计算从而降低算法的运行效率；并且由于分量过多，分量的频带过窄，分量可能无法包含足够的诊断信息。设置分解个数为K=4。惩罚因子α的取值对分解效果也有影响。当K确定时如果α越大，则分量的长度就越窄，根据实际情况本文选择α=200。分解得到的4个分量，如图19所示。

图19 VMD法分解出的4阶IMF时域波形图

使用谱峭度为指标计算各分量信号的全频带谱峭度，并选取谱峭度最大的两个分量相加得到重构信号。重构信号与重构信号的包络谱，如图20所示。

(a) 重构信号

(b) 重构信号包络谱

从图18对比MOMEDA与重构信号可以看出重构信号增强了内圈、外圈和左边带频率处的故障特征，且其它特征也与MOMEDA方法效果持平。因此，重构信号的故障特征更加明显。

3.4 加工轴承外圈故障数据采集说明

为说明MOMEDA与VMD结合的必要性，使用加工轴承外圈故障实验数据并添加-16 db的白噪声模拟轴承早期外圈故障。

加工故障轴承的实验，使用PQZZ-II轴承故障模拟实验台对外圈故障轴承进行实验，如图21所示。轴承故障为电火花加工故障。

图21 PQZZ-II实验台

实验转频fr=15 Hz，实验轴承型号为NU205EM，轴承滚柱数为z=12，滚柱直径为d=7.938 mm，节圆直径为D=38.5 mm,接触角为β=0°。根据式(18)，计算得外圈故障频率为fo=71.443 Hz。

(a) 原始信号时域波形

(b) 原始信号包络谱

试验台故障轴承座的水平位置与垂直位置分别安装有型号为PCB603C01加速度传感器。数采设备为NI PXI硬件和Labview软件平台。采样频率为25 600 Hz，采样时间为4 s。原始时域波形与包络谱如图22(a)、(b)所示。观察图22(b)，可以确定轴承存在外圈故障，这与已知故障类型一致。

3.5 加工轴承外圈故障数据分析

加入-16 db的白噪声后时域波形与包络谱如图23(a)、(b)所示。

(a) 加噪信号时域波形

(b) 加噪信号包络谱

(a) MOMEDA滤波后时域波形

(b) MOMEDA滤波后包络谱

观察图23(b)声几乎将特征频率淹没。使用MOMEDA对加噪信号进行滤波，滤波后的时域波形与包络谱如图24(a)、(b)所示。观察图24(b)，MOMEDA对加噪信号特征频率有所增强，但特征频率附近仍有大量噪声干扰。

使用VMD对滤波后的信号进行分解重构，VMD参数设置与4.3.2节中一致，VMD分解重构后的时域波形与包络谱如图25(a)、(b)所示。

(a) 重构信号

(b) 重构信号包络谱

观察图25(b)，特征频率峰值极为突出。分别计算提出方法各步特征频率与谱总值(0～350 Hz)比值，如图26所示。

图26 3种方法效果比较

观察图26，相比原始信号MOMEDA滤波后特征频率与谱总值比值提升了161%；而重构信号相比原始信号特征频率与谱总值比值提升了400%。可以说明在信噪比极低情况下MOMEDA与VMD分解重构结合对信噪比提升有重要意义。

4 结 论

本文提出了自适应MOMEDA方法并与VMD方法结合对滚动轴承早期故障进行了特征提取。对比传统的MED类方法，自适应MOMEDA避免了MED类方法可能会产生的虚假冲击和局部最优，从而提高了对滚动轴承早期故障特征提取的准确性。相比文献[9]提出的方法，本文提出的方法可以自适应确定最优的滤波器长度从而获得较好的滤波效果。对滤波后的信号进行VMD分解，再根据谱峭度准则得到重构信号，并对重构信号进行包络分析提取了故障特征。实验结果验证了本方法的有效性。