AVMD-IMOMEDA在滚动轴承声学复合故障诊断的应用

周文杰, 周 俊, 柳小勤, 刘 韬

(1. 昆明理工大学 机电工程学院,昆明 650500;2. 云南省先进装备智能制造技术重点实验室,昆明 650500)

滚动轴承广泛用于工业制造,其健康状况直接决定了精度、可靠性、生产效率和设备寿命[1],且有资料统计,旋转机械失效案例中45%～55%是由滚动轴承失效引起的[2]。因此对滚动轴承进行状态检测和故障诊断是非常有必要的。事实上,滚动轴承的声信号中含有大量设备运行状态信息,通过对声信号进行分析处理,可以提取设备运行状态特征,对设备进行状态监测和故障诊断[3]。然而滚动轴承的声信号分析不可避免地受到其他部件和环境噪声的干扰[4],这使得难以恢复与故障相关的重复脉冲和提取故障信息。而且在现实生活应用中,滚动轴承存在不同故障同时存在的复合故障状态。因此研究一种可行的滚动轴承声信号复合故障特征提取方法具有重要意义。

Dragomiretskiy[5]提出的变分模式分解(variational mode decomposition,VMD)可将输入信号分解为一系列具有有限带宽和一定频率的固有模式函数。试验验证VMD克服了经验模态分解、集成经验模态分解(ensemble empirical mode decomposition,EEMD)等方法的缺点[6],并已广泛应用于滚动轴承故障诊断领域[7]。然而VMD的优异性能取决于适当的参数选择。模式数K的不适当预设值会导致过度分解或欠分解[8]。因此如何自适应地确定K是VMD需要解决的重要挑战。

Mcdonald等[9]提出了多点最优最小熵反褶积(multipoint optimal minimum entropy deconvolution adjusted,MOMEDA),可有效地增强故障脉冲的重复性,从而提取故障信息。与最小熵反褶积和最大相关峭度反褶积(maximum correlation kurtosis deconvolution,MCKD)相比,MOMEDA算法无需重采样过程。Cai等[10]通过MOMEDA对组合乘积函数进行降噪,以提取齿轮箱的故障特征。Zhao等[11]引入了改进的带自适应噪声的完整EEMD(complete EEMD with adaptive noise,CEEMDAN)和MOMOEDA,以减少强噪声干扰并从滚动轴承信号中提取缺陷特征。Li等[12]使用MOMEDA优化VMD处理的重构信号,从中提取故障周期脉冲分量。MOMEDA可进一步增强重构信号中与故障相关的周期脉冲,提高故障诊断的精度。然而MOMEDA也存在参数选择问题,如何确定参数周期T和滤波长度L是我们需要考虑的问题。

为了解决VMD和MOMEDA参数选择问题,以及分离提取复合故障相关特征,提出一种基于参数自适应变分模式分解(adaptive variational mode decomposition,AVMD)结合改进多点最优最小熵反褶积(improved multipoint optimal minimum entropy deconvolution adjusted,IMOMEDA)的复合故障声学诊断方法。首先采用综合指标(comprehensive index,CI)解决VMD参数选择问题;然后利用最大加权峭度(weighted kurtosis,WK)对AVMD分解的一系列模态进行有效模态选择及重构,以增强与故障特征相关的脉冲特征信息;最后利用IMOMEDA从重构信号中分离提取周期性脉冲信号,并通过包络解调获取故障特征频率,进而完成滚动轴承声信号复合故障的特征提取。仿真信号和实测声信号的试验结果证明了所提方法的有效性。并将所提出的方法与VMD-MCKD以及传统VMD和MOMEDA方法进行比较,结果表明该方法在故障特征分离提取方面比上述对比方法具有更好的性能。

1 VMD基础理论

VMD将输入信号分解为具有不同中心频率ωk和有限带宽的K个模式。为估计模式的带宽,约束变分问题可用式(1)描述。

(1)

通过使用拉格朗日乘子λ和二次罚项α,将上述问题转化为无约束变分问题。增广拉格朗日量L由式(2)表示。

(2)

(3)

(4)

(5)

当满足收敛条件等式(6)时,停止更新参数的迭代。

(6)

式中,ε为收敛误差。

模态数K的选择对分解结果的影响较大。因此,为了获得良好的去噪效果和计算效率,通常对模态数K进行优化。

2 MOMEDA算法

设传感器采集到的振动信号为x,其表达式为

x=h*y+e

(7)

式中:y为脉冲序列;h为传递函数;e为背景噪声。

McDonald等提出了多D范数(multi d-norm,MDN)对确定位置处的多个脉冲进行反卷积,其最大值为MOMEDA。该算法以非迭代方式选择最佳滤波器f,并基于振动信号x重置脉冲信号y。

(8)

式中,k=1,2,…,N-L。

滤波信号的MDN及其最大问题如下

(9)

(10)

式中:t为目标向量,用于确定目标冲击脉冲的位置和权重。式(10)的求解等价于求解式(11)。

(11)

式中:f=f1,f2,…,fL;t=t1,t2,…,tN-L。

由式(8)、式(10)、式(11)可求得

(12)

式中,k=1,2,…,N-L。

令X0=[M1,M2,…,Mk],则式(12)简化为

‖y‖-1X0t-‖y‖-3tTyX0y=0

(13)

整理得

(14)

(15)

取得特解作为一组最优滤波器,记为

(16)

3 基于AVMD-IMOMEDA的故障诊断方法

3.1 基于综合指标优化VMD参数

信息熵是一种优秀的评判信号稀疏特性指标,熵值大小反应了信号数值的不确定程度,信号不确定度 越大,其信息熵越大[13]。在此基础上,包络熵的概念相应被提出,即为先对信号进行解调运算,计算所得包络信号序列的信息熵,可更好的反映原信号的稀疏特性。平均包络熵(EME)是指在特定参数K下,信号经VMD分解后各模态分量的包络熵的平均值。具体计算公式如下

(17)

式中:K为VMD分解模态个数;Hen(i)为每一个模态分量uk的包络熵。

峭度作为一种无量纲指标,被广泛应用于适应度函数中,但其缺点是忽略了循环平稳性。因此,包络谱峭度被提出[14]。平均包络谱峭度(KMES)是指在特定参数K下,信号经VMD分解后各模态分量包络谱峭度的平均值。具体计算公式如下

(18)

式中:K为VMD分解模态个数;KES,K为第K个模态分量的包络谱峭度。

基于以上分析,结合两者优点,构建一种综合指标,其定义如下

IC=1/KNMES×ENME

(19)

式中:KNMES为归一化后的平均包络谱峭度;ENME为归一化后的平均包络熵。

VMD参数优化的具体步骤如下。

步骤1设置要优化的VMD参数的适当阈值范围。

步骤2采用CI作为目标函数进行循环迭代,以输出最小CI值的参数为最优解。

设定VMD参数K的寻优区间为[2,10],搜索步长为1。

3.2 分量选择

信号经VMD分解成多个子模态后,在区分有效成分、噪声及干扰成分方面,峭度指标和相关系数运用最为广泛[15]。然而,峭度指标只取决于冲击信号的分布密度,作为评判指标可能会忽略振幅较大,具有分散分布的成分。相关系数可以表示各子模态与原信号的相关性,但是易受到噪声的干扰。因此构造最大WK(KW)指标用于区分信号中有效成分、噪声及干扰成分,具体计算公式如下

KW=Kr·C

(20)

式中:Kr为各子模态的峭度指标;C为各子模态与原信号之间的相关系数。

选择具有最大加权峭度值的分量为有效模态分量重构信号。

3.3 基于多点峭度和变步长搜索优化MOMEDA参数

MOMEDA是检测轴承信号的有力工具,可通过精确的非迭代方法增强连续脉冲,然而滤波结果受到反褶积周期T和滤波器长度L的限制。为克服这些缺点,提出一种IMOMEDA方法。

首先为了确定反褶积周期T,引入了多点峭度,可描述为

(21)

当目标矢量t和故障产生的冲击信号y相同时,多点峭度归一化可得

进一步求解

经过标准化的多点峭度可定义为

(22)

多点峭度可使被控处的目标矢量转化为多个脉冲特征,从而更有效地识别出信号中与故障相关的周期性脉冲成分。滚动轴承在实际运行中每转动一圈,会激发出一系列的故障冲击信号,因此多点峭度的峰值通常出现在0.20T,0.50T,0.75T,T,1.50T,2.00T处,由此利用多点峭度确定反褶积周期T。

反褶积周期确定之后,本文以排列熵(permutation entropy,PE)作为目标函数,并使用变步长搜索方法搜索最佳滤波器长度L,其具体步骤如下。

步骤1设定初始搜索范围。

步骤2以大步长S1进行初步全局搜索,选择PE全局峰值对应的滤波器长度L1。

步骤3由步骤2的结果确定新的搜索范围为[L1-S1,L1+S1]。

步骤4以小步长S2进行局部搜索,得到PE峰值对应的滤波器长度L2,而L2即是最优滤波器长度。

基于文献[16]分析,在其基础上本文设定L值的初始搜索范围为[300,2 000],大步长S1为50,小步长S2为1。

3.4 AVMD-IMOMEDA算法步骤及流程

方法流程,如图1所示。其步骤如下。

图1 本文方法流程图Fig.1 The flowchart of the method in this paper

步骤1采集轴承的声信号作为数据输入。

步骤2通过AVMD对输入信号进行分解,得到一系列分量。

步骤3计算每个模式的WK值,选择最大WK值的分量重构信号。

步骤4利用IMOMEDA对重构信号进行处理。

步骤5对IMOMEDA滤波后的信号进行包络谱分析,得到轴承的故障特征频率。

4 仿真信号分析

为验证所提出方法的可行性,构建仿真信号进行验证。其中仿真内圈故障冲击数学模型表达式[17]如下

(23)

式中:s(t)为周期性冲击成分;n(t)为噪声;fr为转频;Ai为以1/fr为周期的调制幅值;CA为中心偏移量,CA=1;h(t)为指数衰减脉冲;τi为第i次冲击相对于平均周期T的微小波动;B为系统的衰减系数;φA和φω分别为初相位;fn为固有频率。

仿真外圈故障冲击信号数学模型为

(24)

式中:A1为外圈位移常数;t0为单周期采样时刻;g为阻尼系数。内、外圈故障仿真信号参数,如表1所示。

表1 故障仿真信号参数Tab.1 Parameters of fault simulation signals

信噪比为-10 dB的内外圈故障仿真信号时域图,如图2所示。内、外圈复合故障仿真信号的时域图和包络谱,如图3所示。由图3可知,由于背景噪声的强干扰,有效的故障特征信息很难识别。因此需要对信号进行处理以提取故障特征频率。

图2 内外圈故障仿真信号时域图Fig.2 Time domain diagram of inner and outer ring fault simulation signal

基于以上分析,利用本文提出的方法对复合故障仿真信号进行特征提取,首先依据以综合指标为CI目标函数的循环算法自适应VMD的参数K,确定参数K为3;再根据最大WK指标选择分量重构信号。重构信号的时域及多点峭度谱如图4所示。多点峭度谱显示,在周期T分别为96及其0.5倍和149及其0.5倍处的峰值明显突出,且符合故障周期的理论计算值,由此可确定复合故障仿真信号中两种不同故障的故障周期。确定故障周期后,使用变步长搜索法优化滤波器长度,最终得到参数组合分别为[96,394],[149,317],据此对重构信号进行MOMEDA分析。

图4 重构信号时域图及多点峭度谱Fig.4 Reconstruction of signal time domain diagram and multipoint kurtosis spectrum

内、外圈故障解卷积信号时域图,如图5所示。内、外圈故障解卷积信号包络谱分别为图6(a)和图6(b)所示;图5、图6中内、外圈故障特征频率及其倍频处明显突出,噪声干扰成分少,由此可知运用本文方法可有效提取复合故障仿真信号中的内外圈故障特征频率。

图5 内外圈故障解卷积信号时域图Fig.5 Time domain diagram of inner and outer circle fault deconvolution signal

图6 内外圈故障解卷积信号包络谱Fig.6 Envelope spectrum of inner and outer circle fault deconvolution signal

5 试验信号分析

试验数据来自实际环境下两个北京声望传声器拾取的QPZZ-Ⅱ旋转机械振动及故障模拟试验台运行时的轴承内外圈复合故障的声信号。故障轴承的缺陷为电火花切割处理,其故障尺寸为1.5 mm×0.5 mm×0.5 mm。故障轴承相关参数为:节圆直径D=39 mm,滚动体直径d=7.5 mm,滚动体数目Z=12,接触角β=0°。根据轴承的参数计算各个特征(故障)频率。当转速为800 r/min,即转频为13.33 Hz,内圈故障特征频率为95.38 Hz,外圈故障特征频率为64.61 Hz。

通过NI SignalExpress采集模块及NI-9234四通道采集卡进行信号采集,采样频率fs=8 192 Hz,采样点数N=8 192。频率间隔Δf=fs/N=1 Hz,将两个传声器相互垂直安装在距故障轴承试验台边沿0.5 m,高度0.5 m的位置,传声器1和故障轴承座侧面中线平齐,传声器2正对轴承主轴。试验台及传声器的布局如图7所示。本次试验分析来源于图7标识的传声器2的数据,滚动轴承故障类型为内外圈复合故障。

图7 故障模拟试验台及传声器布置图Fig.7 Layout of fault simulation test bed and microphone

轴承内外圈复合故障声信号的时域波形图和包络谱,如图8所示。由图8可知,滚动轴承内外圈复合故障声信号在时域中杂乱无章,冲击成分几乎被噪声掩盖。在其对应的包络谱中发现内外圈复合故障特征频率成分完全混在一起,且冲击不明显。因此需进行特征提取。

图8 轴承内外圈复合故障声信号时域图及包络谱Fig.8 Time domain diagram and envelope spectrum of compound fault acoustic signal of bearing inner and outer rings

利用本文方法对内外圈复合故障声信号进行特征提取。首先依据以综合指标CI为目标函数的循环算法自适应VMD的参数K,再根据最大WK选择最优分量重构信号。参数优化如图9所示。由图9可知,当参数K为3时,综合指标获最小值,确定此时参数K为最优值。最优分解各分量的加权峭度值,如表2所示。由表2可选择分量3作为重构信号。

表2 各分量加权峭度值Tab.2 Comparison of AVMD and VMD algorithms

图9 参数优化图Fig.9 Diagram of parameter optimization

内、外圈复合故障声信号经AVMD算法处理再选择具有最大WK值有效分量重构信号的时域图及多点峭度谱,如图10所示。由图10可知,在周期为126.7以及周期为86.5处的峰值明显突出,且符合故障周期的理论计算值,由此可确定复合故障声信号中两种不同故障的故障周期。

图10 重构声信号时域图及多点峭度谱Fig.10 Reconstruction of acoustic signal time domain diagram and multipoint kurtosis spectrum

确定故障周期后,以排列熵为目标函数,使用变步长搜索法优化滤波器长度,如图11所示。图11(a)和图11(b)分别为故障周期86.5和126.7对应的小区间滤波器长度优化示意图。由此结合3.3节和图11最终可得到优化后的两对参数组合分别为[86.5,1 599],[126.7,1 626],据此对重构信号进行MOMEDA分析。

内、外圈故障解卷积声信号时域图,如图12所示。内、外圈故障解卷积声信号包络谱,如图13所示。由图12、图13可知,内、外圈故障特征频率及其倍频处明显突出,噪声干扰成分少,由此可知运用本文方法可有效分离提取复合故障声信号中的内外圈故障特征频率。

图12 内外圈故障解卷积声信号时域图Fig.12 Time domain diagram of deconvolution acoustic signal for inner and outer circle faults

图13 内外圈故障解卷积声信号包络谱Fig.13 Envelope spectrum of deconvolution acoustic signal for inner and outer ring faults

试验信号经传统VMD算法(即任选参数)与AVMD算法处理后所得重构信号的峭度和相关系数,如表3所示。由表3可知,AVMD相比于VMD,其重构的信号峭度值和相关系数都更大,其恰好说明了优化VMD参数的必要性。

表3 AVMD与VMD算法对比Tab.3 Comparison of AVMD and VMD algorithms

传统MOMEDA算法(即任选参数T和L)对图10中的重构时域信号进行处理得到的包络谱,如图14所示。尽管在图14中可发现相应的内、外圈故障特征频率,但存在大量的干扰谱线,且整体来看峰值不明显以及倍频基本淹没在噪声中,其正说明了优化MOMEDA参数的必要性。

图14 基于MOMEDA重构信号包络谱Fig.14 Reconstruction of signal envelope spectrum based on MOMEDA

为验证本方法的优越性,利用VMD-MCKD对试验数据进行分析。试验数据经VMD-MCKD分析后的包络谱,如图15所示。虽然能够找到故障特征频率,但整体来看幅值明显小于本文方法,且背景噪声干扰较大。

图15 基于VMD-MCKD试验信号包络谱Fig.15 Based on VMD-MCKD test signal envelope spectrum

6 结 论

针对滚动轴承声信号中微弱复合故障难以提取的问题,提出了一种AVMD-IMOMEDA算法,经过构造仿真信号和试验信号对其进行验证,得出以下结论。

(1) 通过VMD与AVMD,MOMEDA与IMOMEDA算法的对比,说明了VMD和MOMEDA参数优化的必要性。

(2) 复合故障仿真信号和试验信号分析结果表明本文所提方法能有效准确地分离出滚动轴承内外圈故障,且包络谱中故障特征频率明显突出,噪声干扰成分少。

(3) 采用算法VMD-MCKD对试验信号分析,并将分析结果与所提方法进行了比较,说明了所提方法的优越性。