基于随机共振的轴承早期微弱故障诊断

赵团团，蒋甲丁，臧能义

(1.新疆工程学院工程技能实训学院，新疆乌鲁木齐 830011；2.新疆工程学院科研处，新疆乌鲁木齐 830011)

0 前言

滚动轴承作为现代工业生产中重要的机械组成部分，及时有效的早期故障诊断意义重大[1]。常用的故障诊断方法有小波变换[2]、经验模态分解[3-4]、盲源分离[5]等。对于早期微弱信号，传统滤波降噪方法在滤除噪声干扰的同时会将部分有用信号一同滤除，造成检测精度的降低，而随机共振(Stochastic Resonance，SR)的出现有效解决了这个问题。SR现象是BENZI等[6]在探讨气象相关内容时提出的，主要原理是借助噪声辅助增强微弱周期信号能量而非滤除噪声[7]，在生产尤其是微弱故障信号检测方面应用广泛[8]。

参数调节作为SR最常见的方法，应用广泛。文献[9]将峭度作为适应度函数优化由排列熵筛选出的振动信号；文献[10]以信噪比为目标函数将变分模态分解与SR结合，实现风机微弱故障信号的提取；文献[11]在进行SR时将步长作为参数自适应同步选取，更全面地考虑参数间相互关系；文献[12]引入级联SR有效消除了SR处理后信号的边缘脉冲，进而实现信号频率的准确检测；文献[13]通过对比二阶欠阻尼和一阶过阻尼的变尺度效果可知，二阶欠阻尼在轴承故障诊断中更具优越性。

根据以上文献，本文作者提出了一种新的基于SR的信号检测方法，首先充分利用排列熵的特性筛重构新的故障信号，其次将新的故障信号送入SR系统，并以信噪比为目标函数，采用麻雀搜索算法自适应优化SR参数，最后通过分析输出信号进行故障诊断。仿真结果表明：该方法可有效检测出早期微弱故障，进一步扩展了SR的应用，具有一定工程实用性。

1 理论基础

1.1 二阶双稳态系统

由噪声N(t)和微弱信号S(t)共同作用的二阶Duffing振子方程如下：

(1)

图1 系统势函数U(x) 图2 输入输出信噪比随噪声强度变化趋势

为进一步分析二阶双稳态共振系统的输出特性，将一个振幅a=0.3、频率f=0.01、采样频率为5 Hz的纯正弦信号送入共振系统，图2为k=0.5、a=0.65、b=0.5时的输入输出随噪声强度变化时信噪比的变化。

由图2可以看出：随着噪声强度的增加，输入信噪比逐渐减小，而输出信噪比呈现先增大后减小的趋势，在噪声强度为0.1～0.6时信噪比较大，说明共振系统中的噪声在一定情况下可以增强信号能量，且在一定条件下达到极值。

1.2 基于排列熵理论的筛选原则

由文献[14-15]可知：熵的计算简单、抗噪能力强，对复杂系统的动力学突变具有敏感性[16]，在机械故障诊断中应用广泛。

将长度为N的时间序列{x(i),i=1,2,…,N}进行相空间重构，得到如下公式：

(2)

式中：m为嵌入维数；τ为时延。

对X(i)={x(i),x(i+τ),…,x(i+(m-1)τ)}进行升序排列操作后，X(i)={x(i+(j1-1)τ)≤x(i+(j2-1)τ)≤…≤x(i+(jm-1)τ)},若存在x(i+(ji1-1)τ)=x(i+(ji2-1)τ)，则根据j值大小排序，即当ji1

S(l)=[j1,j2,…,jm]

(3)

(4)

其中：Hp(m)在Pl=1/m!时达到最大值ln(m!)，因此可将Hp(m)标准化为

Hp=Hp(m)/ln(m!)

(5)

其中：0≤Hp≤1，Hp数值大小代表时间序列的复杂度与随机度：Hp越小，时间序列越规则，含有周期脉冲信号的可能性越大；反之Hp越大，时间序列越杂乱，含有周期脉冲信号的可能性越小。

基于以上排列熵特性，本文作者构造了一种新的信号筛选判据——排列熵判据(Permutation Entropy Criterion，PEC)，即剔除复杂度高于原始信号50%的分量，将剩余的分量合并成新的待测信号，保留了大部分周期分量，其数学表达式如下：

fpe(i)

(6)

PEC=sum(fpe)

式中：fpe为排列熵值；i为第i个IMF分量

1.3 麻雀搜索算法

麻雀搜索算法(Sparrow Search Algorithm，SSA)[17]是薛凯健等于2020年提出的一种对麻雀觅食和反捕捉的新型仿生优化算法，其性能优于大多数智能算法[18]，SSA可简化为带侦察预警机制的发现者-加入者模型。设有N只麻雀在一个D维搜索空间中，其中第i只麻雀的位置为Xi=[xi1,…,xid,…,xiD]，其中i=1，2，…，N，xid为第i只麻雀在第d维的位置。发现者个数较少，一般为种群的10%～20%，其位置更新公式如下：

(7)

式中：L为大小为1×d、元素值都为1的矩阵；Q为服从标准正态分布的随机数；t为当前迭代次数；T为迭代的最大次数；α为(0，1]之间的均匀随机数；R2∈[0,1]、fST∈[0.5,1]，分别表示预警值和安全值。

当R2fST时，侦查麻雀发现捕食者并作出报警信号，种群立刻调整搜索策略，迅速向安全区域靠拢。

除了发现者，加入者位置更新公式如下：

(8)

式中：Xworst为当前全局最差位置；Xp为发现者目前的全局最优位置；A为大小为1×d、元素值为1或-1的矩阵，且满足A+=AT(AAT)-1。根据公式，当i>n/2时，由于适应度值较低的第i个加入者没有获得食物，此时需飞向其他地方寻找食物，以获得更多的能量。

1.4 基于SSA的风机轴承未知微弱故障信号检测

将信噪比作为适应度函数进行自适应参数寻优，信噪比最大时的参数条件下SR输出性能最优。信噪比表达式如下所示：

(9)

式中：n为故障特征频率附近的噪声功率谱点数；P0为故障特征处对应的功率谱；Pi为各计算点处对应的功率谱。

对于公式(1)的数值计算，主要采用四阶龙格-库塔法求解，数学表达式如下：

(10)

结合公式(1)与公式(10)可知，信号输出除了受到系统参数k、a、b的影响外，步长h也会影响输出。而目前大多SR方法只考虑参数而忽视了步长，本文作者充分考虑系统参数及步长h的影响，采用SSA进行自适应参数寻优。(k*，a*，b*，h*)=argmax CKZCR (k，a，b，h)，其中(k*，a*，b*，h*)为输出最优参数。基于SSA的轴承故障诊断流程如图3所示。

图3 基于SSA的风机轴承故障诊断流程

2 工程应用

为了验证该方法是否有效，选取凯斯西储大学(CWRU)滚动轴承的实验数据集进行仿真，试验台如图4所示，表1为轴承的相关参数。

图4 西储大学轴承试验台

表1 6205-2RS轴承相关参数

公式(11)为轴承故障特征频率计算公式：

fBPI=rn/120(1+d/Dcosα)

fBPO=rn/120(1-d/Dcosα)

(11)

式中:n为滚珠的个数；d为滚珠的直径；r为转速；D为节径。

已知采样频率fs=12 kHz，采样点数N=8 192。内圈故障频率159.93 Hz，尺度变换系数m取2 400，同时为体现所提方法优越性，向信号中加入强度为1的高斯白噪声。原始信号时域频谱图及加噪后的频谱如图5所示。

图5 内圈信号时域频谱

可以看出，采集的原始数据与信号包含了大量轴承信息，但无法直接识别出故障特征的频率，加噪后信号被完全淹没。应用文中所提出方法对信号进行处理，首先对带强噪的信号进行EEMD经验模态分解(见图6)并采用SEC判据，计算出各IMF分量的排列熵值(见图7)，并将值小于原始信号50%的信号滤除。

图6 内圈故障信号EEMD分解后IMF分量

图7 各IMF分量排列熵值

图8(a)为重构信号时域图，已经可以看到明显的尖峰。由图8(b)可知：信号集中在0～400 Hz区间，故障频率159.7 Hz与内圈故障频率一致，但周围干扰噪声使得故障特征难以辨别。将筛选后的SEC信号输入SR系统，寻优得(k，a，b，h)=(0.747 4，1.593 3，0.156 3， 0.2)，经自适应SR后信号如图8(c)与8(d)所示，可以看到信号经SR参数调整后，周期特征明显，且频谱可看到明显的故障特征频率。

图8 SR处理前后重构信号时域频谱

将所提方法与经典滤波算法——快速谱峭度[19]做对比，得到内圈故障谱峭度图如图9所示，最大的频段位于Level6.6，中心频率为2 812.5 Hz。将筛选出的频带进行自适应SR参数优化后，其包络如图10所示。

图9 内圈故障信号快速谱峭度图

图10 快速谱峭度经SR后包络图

由图10可以看出：经过快速谱峭度滤波后的信号经SR后故障特征频率为161.8 Hz，与内圈故障理论值159.93 Hz基本一致，但周围的干扰使得故障特征不明显，可见，快速谱峭度在强噪声背景下的检测性能会受到干扰。

3 结论

(1)针对早期微弱信号难以提取，利用排列熵对周期信号的敏感特性提出一种新的早期微弱故障信号的滤波方法。

(2)为仿真传统滤波方法噪声早期微弱故障信号丢失这一问题，采用随机共振噪声辅助增强微弱信号。

(3)仿真结果表明：参数调节后的SR系统可有效放大微弱周期信号。文中提出的方法能准确诊断出早期轴承的未知微弱故障，具有一定的实用价值。