基于平方包络谱负熵准则的轴承早期复合故障特征提取方法

陈 鹏， 赵小强,2,3

(1.兰州理工大学 电气工程与信息工程学院, 兰州 730050；2.甘肃省工业过程先进控制重点实验室，兰州 730050；3.兰州理工大学 国家级电气与控制工程实验教学中心，兰州 730050)

在石油化工、有色冶金等领域的生产装置中，鼓风机、压缩机、风机和泵等旋转设备是整个生产装置的“心脏”，也被称为“关键机组”。这些关键机组具有转速高、功率大、无备机和维修周期长的特点，如一旦发生故障而停机检修，则会对整个企业生产造成牵一发而动全身的重大影响，而轴承是关键机组在长期运行过程中最易发生故障的零部件，如果不对其早期故障进行及时监测和诊断，则会引发严重故障而导致停机检查和维修，从而造成严重经济损失。而在早期故障中，早期复合故障由于具有多信号耦合、信号微弱和强噪声干扰的特点，相对早期单故障和成熟复合故障而言更加难以辨识[1]。因此，研究能够提取轴承早期复合故障特征的有效方法一直是故障诊断领域的难点[2]。

目前，以时频分析为主的信号分解方法在轴承的降噪和故障诊断中得到了广泛应用[3]。如Huang等[4]提出了经验模态分解(empirical mode decomposition,EMD)方法，但是EMD方法存在对噪声敏感、模态混叠和端点效应的问题，因此有学者们针对EMD的问题，提出了集总经验模态分解(ensemble empirical mode decomposition，EEMD)[5]、互补集合经验模态分解(coplementary ensemble empirical mode decomposition,CEEMD)[6]。同时，学者们提出了其他一些自适应分解方法，如优化集合局部均值分解(ensemble local mean decomposition,ELMD)[7]、局部特征尺度分解(local characteristic-scale decomposition,LCD)[8]和固有时间尺度分解(intrinsic time scale decomposition,ITD)[9]等。但这些方法都是基于递归分解的原理，不可避免地存在一定的模态混叠和端点效应[10]。而Dragomiretskiy等[11]于2014年提出了一种新的信号自适应分解方法—变分模态分解(variational mode decomposition，VMD)，近年来被广泛应用于旋转机械的振动信号处理。研究结果表明，VMD可以在一定程度上改善模态混叠和端点效应[12]，但是该方法在早期复合故障诊断领域的效果欠佳[13]。对此，Apostolidis等[14]于2016年在群滤波(swarm filter, SWF)的基础上，提出了一种新的信号智能分解方法—群分解(swarm decomposition, SWD)，它通过迭代滤波，可将多个频率相近的复合信号自适应地分解为若干单一模态的振荡分量，具有较高的频率区分能力，在解决多模式混合信号的分解方面具有突出的优势。李娟等[15]将SWD分解方法与平均差值算子方法相结合实现了齿轮箱的复合故障诊断。但是，未对阈值参数值PSth和Tth设置，导致很难达到预期效果，故建立能够提取复合故障特征的准则，获得到合理阈值参数，从而使 SWD将不同故障特征分解在不同的分量下，对于提高SWD的分解能力具有重要意义。

基于以上分析，本文构建基于平方包络谱负熵的优化准则，采用在全局、局部搜索能力和稳定性方面较优的改进蝗虫优化算法(improved grasshopper optimization algorithm, IGOA)对SWD中的阈值参数进行优化，实现对轴承早期复合故障特征的提取。

1 SWD方法

SWD作为一种智能分解方法，在参数化SWF的基础上，可以从复合信号中迭代提取出主振荡分量(oscillatory component ,OC)[14]。SWF的基本原理是将输入信号x[n]作为群体的捕食路径，输出群体轨迹。在群滤波的模型中，猎物在第n步的位置表示为pprey[n]，M为群的数量，群的第i个个体在第n步的位置和速度分别为pi[n]和vi[n]。群个体在猎物过程受两种作用力的控制，第一种作用力是由自身的驱动力产生，第二种作用力是由相互之间的吸引或排斥而产生。分别定义为式(1)和式(2)

(1)

(2)

式中，f()函数定义为

(3)

式中：sgn()为符号函数；ln()自然对数；d|为两个群个体之间距离的绝对值；dcr为2个群个体之间相互不影响的临界距离，即f(dcr)=0。一般来说，dcr可以用来控制群个体的分布，其大小为输入信号的均方根值。

为跟踪猎物，群个体必须通过改变自己的位置来更新其状态，其第i个个体在第n步的速度和位置可表示为

vi[n]=vi[n-1]+δ(Fdr-i[n]+Fco-i[n])

(4)

pi[n]=pi[n-1]+δvi[n]

(5)

式中，δ为控制群个体的灵活性。

SWF的输出为群的轨迹，可以通过式(6)表示

(6)

式中：β为尺度因子，一般取β=0.005；群参数β和M为控制群行为的关键决定因素。对于不同的信号，通过式(7)来选择δ和M的值

(7)

式中：Yδ,M(k)和S(k)为序列Yδ,M(n)和S(n)傅里叶变换的幅值；Yδ,M(n)为群滤波在参数为δ和M时的输出；S(n)为一个预设的非平稳信号，该信号可以是包含多个分量的复合信号。参数δ和M与各M频率之间的关系为

(8)

(9)

通过对SWF的迭代执行，不断分解最新的振荡分量，直到当输入信号中不包含任何足够能量的振荡模式时算法终止，即当连续两次迭代的偏差小于阈值Tth时，SWD迭代终止。在每一次迭代中，选取能量谱密度幅值较大的频带作为提取新模式的候选分量。

由于SG(savitzky-golay)滤波器具有平滑能量谱的作用，为了提高计算效率，在找到最高峰值之前对分解信号进行SG滤波具有重要意义。同时，为了减少搜索区域，预先设置一个合适的峰值阈值PSth，其最优频率ωm可以通过下式估计

(10)

(11)

式中，S为使用SG滤波xit[n]信号后的傅里叶变换。

通过对SWD分解原理的分析，SWD分解算法的具体步骤总结如下：

步骤1输入信号x[n]并初始化阈值PSth和Tth；

步骤2离散化xit[n]←x[n]，xit[n]为离散数据序列，it为序列号。首先，设置it=0，然后指定y0[n]←xit[n]；

步骤3平滑y0[n]的能量谱，估计最优频带；

步骤4初始化群参数和δ；

步骤5当i=1时，通过SWF进行群滤波，yi[n]←FSWF(xit[n],M,δ)。

在SWD的分解过程中，阈值PSth和Tth对于SWF的分解过程起着关键作用，其大小决定了是否能够合理分解。如果设置不合理，则会导致分解时间增加或分解分量增多。

2 提出的早期复合故障诊断方法

2.1 建立平方包络谱负熵优化准则

当滚动轴承发生故障时产生循环瞬态冲击，并引起振动信号瞬时能量的波动，而平方包络方差可以用来检测信号时域内的瞬时能量波动[16]。对于长度为N的信号x[n]，其平方包络SSE(x[n])可表示为

(12)

但是平方包络方差不能检测周期性冲击，而故障产生的冲击往往是周期性的，对于周期性冲击可以通过平方包络谱来检测，平方包络谱定义为[17]

SSES(f)=TDFT(SSE(x[n]))

(13)

信息熵由于可以衡量信号的不规则性而得到广泛应用，其也可以作为用来测量信号周期性冲击强度的指标。因此，根据信息熵定义平方包络谱熵为

(14)

式中，〈〉为取均值。由于平方包络谱熵随着信号的周期性冲击越强，其值越小。为了使其变化性同周期性冲击强度变化一致，定义平方包络谱负熵为[18]

(15)

在智能算法优化参数的过程中，建立一个简单、鲁棒的优化准则对于提高SWD的分解能力非常重要[19]。滚动轴承振动信号的平方包络谱在故障特征频率及其多次谐频处均有峰值，即故障特征频率在平方包络谱中是周期分布的。因此本文根据以上分析，提出以振动信号的平方包络谱负熵为准则，进行自适应SWD分解，将分解信号分量的平均平方包络负熵最大对应的阈值参数作为自适应SWD的分解参数，即

(16)

式中，k为分解后模态分量的个数。

2.2 蝗虫优化算法及改进

Saremi等[20]提出了蝗虫优化算法(grasshopper optimization algorithm, GOA) 。该算法由于原理简单、收敛速度较快以及平衡全局和局部搜索能力强而得到一定应用，其数学模型为

Xi=r1Si+r2Gi+r3Ai

(17)

式中：Xi为第i只蝗虫的位置；Si为在群体中第i只蝗虫与其他蝗虫的相互影响；Gi为第i只蝗虫的重力引力；Ai为第i只蝗虫所受风力影响；r1，r2和r3为[0,1]的随机数，表示Si,Gi和Ai受不同环境影响的系数。

(18)

式中：dij=xj-xi|表示第i只蝗虫和第j只蝗虫之间的距离；dij=xj-xi|/dij为第i只蝗虫和第j只蝗虫之间的单位向量；s(r)群体中蝗虫相互影响力，其函数表达式为

(19)

式中：f为吸引力常数；l为吸引力长度尺度参数。在蝗虫群体觅食的过程中，形成了舒适区、排斥区和吸引区三种互区域，三种区域的范围通过s(r)的值来界定：当s(r)>0时，蝗虫间相互吸引，称r的取值范围为吸引区；当s(r)<0时，蝗虫间相互排斥，称r的取值范围为排斥区；当s(r)=0时，蝗虫之间既不相互吸引也不相互排斥，称r的取值为舒适距离。但当r的取值过大时，s(r)的取值同样接近于0，而此时r的取值并非舒适距离。此外，f和l二者的取值情况可以控制吸引区、排斥区及舒适区的分布情况，且二者取值通常较小，本文中取f=1.5，l=0.5，式(17)中G分量的表达式为

Gi=-geg

(20)

式中：g为重力常数；eg为朝向地心的单位向量。A分量的计算公式为

Ai=-uew

(21)

式中：u为漂移常数；ew为风向的单位向量。在式(17)中代入S，G和A的值，得到

(22)

式中：ubd为d维度变量的上限；lbd为d维度变量的下限；Td为目标中d维度变量的值，最外侧的c用于平衡全局与局部寻优能力，内侧的c用于控制舒适域、排斥域及吸引域的伸缩，c的值按线性递减，其表达式为

(23)

式中：cmax为最大值；cmin为最小值；l为当前迭代次数；L为最大迭代次数。

在此，通过改进的圆混沌映射来代替式(23)来更新参数c的变化，从而提高GOA算法的全局搜索能力和稳定性，其表达式为

(24)

2.3 早期复合故障诊断实现步骤

通过平方包络谱负熵准则，建立IGOA优化SWD的复合早期故障诊断方法如图1所示。其具体步骤如下：

图1 提出的早期复合故障诊断流程Fig.1 The flow of early composite fault diagnosis

步骤1输入采集的轴承振动信号和初始化GOA的参数，蝗虫的种群N=20、最大迭代次数L=10和阈值参数ubd=0.99，lbd=0.01，根据设定参数产生初始化蝗虫群体；

步骤2计算输入信号的能量谱，估计最优频带，初始化群参数M和δ；

步骤3进行SWD迭代分解，并计算分解模态分量的适应度值；

步骤4如果l≥L，迭代准止，否则转入步骤3，且l=l+1，更新当前最优的个体；

步骤5保存输出最优阈值参数，构建优化群分解(optimized swarm decomposition，OSWD)方法；

步骤6利用OSWD方法分解振动信号，获得所有模态分量；

步骤7对分解后的分量进行包络谱分析，得到故障特征频率。

3 试验验证

3.1 仿真试验

为了验证所提方法的准确性和有效性，构建一组周期性脉冲序列信号并添加高斯白噪声，其仿真信号模型为

y(t)=s(t)+n(t)

(25)

式中：s(t)为模拟故障冲击的周期性冲击分量；固有频率fn=4 000 Hz；采样点数N=4 096；采样频率为20 kHz；n(t)为高斯白噪声。其中：外圈位移常数A0=3；阻尼系数g=0.05；故障特征频率为100 Hz；产生的模拟外圈故障振动信号如图2所示。内圈位移常数A0=2.5，阻尼系数g=0.1，故障特征频率为125 Hz，产生内圈模拟故障振动信号如图3所示。通过内外圈模拟信号并添加一组随机扰动噪声复合而成的故障信号，如图4所示。与内外圈单故障信号相比，复合故障信号中反映内外圈故障的周期性冲击基本被淹没。对图4所示的复合故障仿真信号直接进行包络分析得到如图5所示的结果。从图5可知，外圈故障及倍频较为突出，内圈故障特征较微弱，容易被忽略，可见直接进行包络谱分析实现不了内外圈的复合故障诊断。

图2 内圈故障时域波形Fig.2 Time domain waveform of inner ring fault

图3 外圈故障时域波形Fig.3 Time domain waveform of outer ring fault

图4 复合故障信号时域波形Fig.4 Time domain waveform of composite fault signal

图5 复合故障信号包络谱Fig.5 Envelope spectrum of composite fault signal

在此，设置蝗虫的种群数N=20、最大迭代次数L=10以及阈值参数ubd=0.99，lbd=0.01。通过提出的OSWD方法对内外圈复合故障仿真信号进行故障特征提取，首先经过OSWD分解后得到的分量如图6所示。从图6可知，故障信号被分解为2个分量，其次对分解后的分量进行包络谱分析如图7所示。由图7可知，在第1个分量中可以清晰地得到内外圈故障特征频率的一倍频和二倍频；在第2个分量中同样可以观察到内外圈故障特征频率的一倍频和二倍频。与理论故障频率进行对比，易判断出内外圈发生了故障。因此，OSWD对内外圈复合故障特征的分离是有效的。

图6 OSWD分解后OC分量Fig.6 OSWD decomposition component

图7 OSWD分解后OC分量包络谱Fig.7 Envelope spectrum of OSWD decomposition component

为了验证本文所提方法的有效性，设置SWD中2个阈值参数为默认值0.2，进行内外圈复合故障信号的分解和特征提取。首先，对复合故障仿真信号进行分解，分解后的分量和包络谱如图8和图9所示。由图可知，复合故障仿真信号被分解为5个分量，在第2个分量中得到了外圈故障特征频率的一倍频和内圈故障特征频率的一倍频与二倍频；在第3个分量中得到了内外圈故障特征频率的一倍频和二倍频。因此，通过SWD方法也可以判断出内圈产生故障，但是SWD的阈值参数是手动选择而不具备自适应性，同时分解得到了3个冗余分解分量，这也必然导致分解过程中计算时间的增加。

图8 SWD分解后OC分量Fig.8 SWD decomposition component

图9 SWD分解后OC分量包络谱Fig.9 Envelope spectrum of SWD decomposition component

为了进一步验证OSWD的优越性，本文继续采用当前在模态混叠和端点效应问题方面较有优势的VMD分解方法对复合故障仿真信号进行分解，设置参数k=6，α=8 000，其中k一般认为值为5～6时可以实现较好的分解，α为VMD中的默认值。分解得到的分量和分量的包络谱如图10和11所示。由图可知，在第2个和第3个分量的包络谱中得到了内外圈故障特征频率的一倍频和二倍频。但是选取不同的[k,α]组合参数，在VMD分解中得到的结果不同，此结果为本文经过反复选择参数测试得到的VMD分解提取效果较优的结果。

图10 VMD分解后分量Fig.10 VMD decomposition component

图11 VMD分解后分量包络谱Fig.11 Envelope spectrum of VMD decomposition component

通过上述三种方法对仿真信号进行复合故障特征提取分析，可以得知OSWD在复合故障诊断方面具有一定的优越性，能够实现复合故障的诊断，有必要进行进一步的应用研究。

3.2 工程验证

将所提方法应用于某石油化工企业生产装置中关键机组风机轴承的早期复合故障诊断，验证该方法的有效性和优越性。在风机运行过程中，其型号为SKF6030C3的风机驱动端轴承内外圈发生早期磨损，如图12所示。从图12可知，其故障程度还比较轻微，属于早期故障，相比内圈，外圈发生了程度更轻微。振动信号通过安装在驱动端的水平和垂直加速度传感器获得，其中电机转速为2 985 r/min，信号的采样频率为25.6 kHz，数据采样点为16 384。电机驱动风机运行的系统结构图如图13所示。该轴承不同部件的理论故障特征频率，如表1所示。

图12 内外圈早期复合故障Fig.12 Early composite failure of inner and outer rings

图13 关键机组风机运行系统结构图Fig.13 The system structure of key unit operation

表1 轴承故障特征频率

风机振动信号时域波形如图14所示。从图14可知，不太规律的周期性冲击，但仅通过时域波形很难判断出内外圈故障。对其直接进行包络谱分析如图15所示。从图15可知，仅有内圈故障特征频率的一倍频比较突出而外圈故障可能是因为内圈故障发生后才产生而比较微弱，因此故障特征频率难以发现。

图14 振动信号时域波形Fig.14 Time domain waveform of vibration signal

图15 振动信号包络谱Fig.15 Envelope spectrum of vibration signal

设置IGOA的初始化参数与仿真案例中的参数相同，通过IGOA优化SWD得到的最佳参数为[0.950 55,0.202 66]。应用提出的OSWD方法对风机轴承的振动信号进行分解得到的分量，如图16所示。从图16可知，振动信号被分解为2个分量，然后对2个分量进行包络谱分析得到的结果，如图17所示。在分解后的第1个分量的包络谱中较为突出的得到外圈故障特征频率。同时，在第2个分量的包络谱中可以明显的看到转频、转频倍频和外圈故障特征频率。由此可见，通过本文提出的方法可以实现关键机组轴承内外圈早期复合故障的特征提取和诊断。

图16 OSWD分解后OC分量Fig.16 OSWD decomposition component

图17 OSWD分解后OC分量包络谱Fig.17 Envelope spectrum of OSWD decomposition component

为验证本文所提方法的有效性和优越性，采用SWD分解方法对关键机组轴承振动信号进行分析，按原始文献设置阈值参数为默认值0.2，振动信号被SWD分解为4个分量如图18所示。对其进行包络谱分析如图19所示。由图19可知，在第1个分量的包络谱中获得了外圈故障特征频率，在第2、第3、第4个分量中都得到了内圈故障特征频率和转频的倍频，实质上在第2个分量已经得到了内圈故障特征频率，但是因为未通过建立的平方包络负熵准则优化来实现自适应分解导致分解后产生冗余分量，增加了分解时间。

图18 SWD分解后OC分量Fig.18 SWD decomposition component

图19 SWD分解后OC分量包络谱Fig.19 Envelope spectrum of SWD decomposition component

为了进一步对比分析，设置VMD参数k=6，α=8 000，继续采用VMD分解方法对采集的风机振动信号进行分析，得到分解后的分量和分量包络谱，分别如图20和图21所示。由图20可知，在分解后第2个分量中得到了外圈的故障特征频率，而内圈故障特征频率没有在任何分量中有效提取，因此通过VMD方法无法判断出内圈故障，故无法实现复合故障特征的有效提取。

图20 VMD分解后分量Fig.20 VWD decomposition component

图21 VMD分解后包络谱Fig.21 Envelope spectrum of VMD decomposition component

通过以上试验方法的对比分析可以得到，通过本文提出的OSWD方法对于关键机组轴承早期复合故障振动信号进行自适应分解，能够实现早期复合故障特征的有效提取，解决了早期复合故障信号由于微弱、多信号相互耦合以及强噪声干扰导致的故障特征难以提取的问题，实现了滚动轴承早期复合故障的准确诊断。

4 结 论

本文提出了基于平方包络谱负熵准则的轴承早期复合故障特征提取方法，主要结论如下：

(1)针对轴承早期复合故障诊断中故障特征难以提取的问题，提出了基于OSWD的复合故障特征提取方法。

(2)建立了能够检测周期性故障冲击能量的平方包络谱负熵。

(3)就阈值参数在SWD在分解过程中对SWD分解效果影响较大的问题，通过建立的平方包络谱负熵准则，应用IGOA优化SWD实现了阈值参数组合优化。

(4)通过与SWD、VMD对比分析表明，该方法在复合故障诊断方面的性能更加优越，具备较好的应用性。