基于增强最大二阶循环平稳盲解卷积的滚动轴承弱故障特征提取
时间：2023-06-25 08:55:04 来源：雅意学习网 本文已影响 人 
陈兆峰，高伟鹏，朱丹宸
(1.海军研究院，北京 100161；
2.海军士官学校，安徽 蚌埠 233012)
滚动轴承安装在各类旋转机械中，在工业生产中得到广泛使用，并发挥着重要作用。然而，滚动轴承常工作在高温、高压和高转速的环境下，极易产生故障，给机械设备的安全稳定运行带来较大的不确定性，因此急需对滚动轴承的故障特点、故障诊断方法进行更深入的研究。目前，振动监测方法是滚动轴承故障诊断中的常用方法，但也存在一些问题：1)设备振源多，导致滚动轴承故障信号常常淹没在强背景噪声中；
2)为了不破坏设备结构，测点往往只能选取在远离故障轴承的机壳外侧，故障点与测试点之间存在复杂的传递路径，进一步加剧噪声干扰。
事实上，滚动轴承故障特征提取所要解决的就是在强背景噪声干扰下提取弱故障特征。为解决这个难题，专家学者们提出了很多有效的方法，如共振解调[1-2]、经验模式分解[3]、变分模态分解[4-5]、形态学滤波[6-7]等。上述方法都能够对滚动轴承的故障状态进行判断，但也存在一定问题，如强背景干扰下共振频带的准确选取问题，分解参数设置对变分模态分解的影响问题，形态学算子以及结构参数的确定问题等。考虑传递路径的影响，为更有效地提取故障源信号，盲解卷积算法近些年得到了快速发展，常见的包括最小熵卷积(Minimum Entropy Deconvolution, MED)[8]、最大相关峭度解卷积(Maximum Correlated Kurtosis Deconvolution, MCKD)[9-10]、自适应多点最优最小熵解卷积(Multipoint Optimal Minimum Entropy Deconvolution Adjusted, MOMEDA)[11]等。这些方法都是以特征指标最大为目标进行解卷积，去除无关成分的干扰，进而获取故障特征更突出的信号；
但这些算法也受到众多参数的制约，如MCKD中移位数M的选取，MOMEDA中解卷积周期的确定，以及广泛存在的滤波器长度选取问题。针对这些问题，专家学者们对算法进行了众多改进，提升了实际使用效果[12-13]。近来，最大二阶循环平稳盲解卷积(Maximum Second-Order Cyclostationarity Blind Deconvolution, CYCBD)算法[14]的提出进一步增强了盲解卷积算法的有效性，它将最大二阶循环平稳指标(Second-Order Cyclostationarity, ICS2)作为分析目标，提取包含滚动轴承周期性故障冲击特征的源信号；
然而，基于CYCBD的特点，为保证分析效果，准确预估滤波器长度和特征频率值尤为重要，但如何更为准确地选取参数仍然值得进行深入探讨[15-17]。
综上所述，基于CYCBD的特征提取能力，本文提出一种增强最大二阶循环平稳盲解卷积(Enhanced CYCBD, ECYCBD)的方法以实现强背景噪声干扰下的滚动轴承弱故障特征提取。ECYCBD以重加权峭度值作为适应度函数，预先确定滤波器长度和特征频率的搜索范围；
利用鲸鱼优化算法获得最优解，确保CYCBD的分析效果；
利用Teager能量算子进一步突出故障成分，实现滚动轴承故障特征的准确提取。
1.1 最大二阶循环平稳盲解卷积
假设观测信号为x，目标源信号为s0，则盲解卷积算法的计算过程可以表示为
s=x*h≈s0，
(1)
式中：s，h分别为通过解卷积得到的源信号和滤波器；
*为卷积算子。为了方便计算，可以用矩阵形式进行表示
s=Xh，
(2)
(3)
式中：L和N分别为s和h的长度。
基于此，最大二阶循环平稳指标可表示为[14]
(4)
(5)
(6)
E=[e1,…,ek,…,eK]，
(7)
(8)
式中：H为矩阵的共轭转置；
RXWX，RXX和W分别为加权相关矩阵、相关矩阵和加权矩阵；
diag为构建对角矩阵；
k为样本数；
Ts为故障周期的长度。
基于循环频率与信号能量波动相关的特性，对于轴承、齿轮故障信号，循环频率值与故障现象、频率相关，循环频率值可设定为k/Ts[14]。此时，向量ek任意一项的上标中，k(N-1)/Ts可改写为kfstN-1/Ts，其中tN-1为第N-1个数据对应的时间，fs为采样频率，则fs/Ts即为滚动轴承的故障特征频率。
求解下列方程
RXWXh=RXXhλ，
(9)
最大的特征值λ即为最大的ICS2值[14]。
由于算法运行过程中h的初始值并不确定，需要预先设定以获得初始的加权矩阵W。受到初始化过程随机性的影响，需要进行如下迭代过程以获得最大ICS2值：1)利用尤尔—沃克方程计算自回归模型的系数，并作为白化滤波器对h进行初始化；
2)借助信号x和h计算得到矩阵W；
3)通过(9)式计算得到最大特征值λ及其特征向量h；
4)返回第2步，将特征向量h代入重新进行计算直至收敛。
本文的收敛条件定义为：1)迭代次数达到100次；
2)|λi-λi-1|/|λi-1|<0.001，其中λi，λi-1分别为本次和上次迭代中获得的特征值。
上述计算过程的输入信号和输出信号均只有一个，为提高分析效果，CYCBD可用于多个信号输入的情况，从而综合利用多信号特征。假设有Q个观测信号xq，通过Q个滤波器hq同时对相应的观测信号进行解卷积计算，则最终的源信号s是多个信号的解卷积结果之和。
(10)
此时，整个算法的分析过程仍可利用上述步骤，通过(9)式获得最大的ICS2值。此时获得的滤波器h是多个与输出相关的滤波器系数的组合[14]，即
h=[h1，…，hq，…，hQ]T，
(11)
第q个输出对原信号的贡献可表示为
sq=Xqhq，
(12)
矩阵Xq的构造方法与(3)式中一致且与信号xq相关，最终的滤波信号为
(13)
1.2 鲸鱼优化算法
鲸鱼优化算法(Whale Optimization Algorithm,WOA)[18]模拟了鲸鱼围捕猎物的行为，主要包括鲸鱼在觅食过程中的3种行为，分别是包围猎物、搜寻猎物和气泡网环形游动驱赶猎物并进行位置更新，这几种行为可通过数学表达式表示。
1)包围猎物。第t+1代中备选解的位置X(t+1)可以通过(14)式进行更新
(14)
(15)
式中：X*(t)为第t代中的最优位置；
a为动态变量，随迭代次数的增加从2至0线性递减；
r为[0,1]范围内的随机数。
2)搜寻猎物。搜索猎物的过程与包围猎物的过程非常类似，不同的是，搜寻猎物的过程不再使用当次迭代的最优解X*(t)，而是随机选取了一个备选位置Xrandom(t)
(16)
3)螺旋式位置更新。鲸鱼在围捕猎物的同时会喷出气泡网用以驱赶猎物，同时不断更新自身位置。使用气泡网时，鲸鱼位置的更新公式为
(17)
式中：D′为第t代的最优位置与任意备选位置的距离；
b为常数；
l为[-1,1]范围内的随机数。
鲸鱼优化算法的总过程如下：
1)初始化鲸鱼种群位置Xi(i=1,2,…,N)，计算各鲸鱼所在位置对应的适应度值，X*(t)即为最优适应度值对应的位置。
2)设定参数a，A，C，l，p的值，本文计算时参数b设定为1。
3)若p<0.5，当|A|<1时，利用(14)式对鲸鱼当前的位置进行更新，否则利用(16)式对鲸鱼当前的位置进行更新；
若0.5≤p<1，则用(17)式更新鲸鱼的当前位置。
4)利用新位置更新a，A和C，并重新计算更新后鲸鱼位置的适应度值，重新确定最优位置X*(t)。
5)重复步骤3和4直到达到最大迭代数。
1.3 重加权峭度值
本文使用一种重加权峭度指标[19]解决峭度指标容易受到少量冲击干扰成分影响的问题，该指标的计算过程如下：
1)将采集到的信号x共划分n段，得到的子信号分别记为x1,x2,…,xn。
2)计算各子信号的峭度值k1,k2,…,kn；
对n个峭度值进行升序排序，得到峭度序列K。
3)计算每个峭度值所占比
(18)
4)重加权峭度值为
(19)
1.4 Teager能量算子
Teager能量算子本质上是对信号动态成分总能量的估计，包括信号的瞬时值和微分形式，由于其计算简便，得到了广泛使用。
对于连续信号x(t)，可将其Teager能量算子表示为
(20)
对于离散时间信号x(n)，其Teager能量算子为
ψ[x(n)]=x2(n)-x(n-1)x(n+1)。
(21)
由于Teager能量算子可以对瞬态冲击成分进行增强，因此适用于突出滚动轴承的故障特征成分。
1.5 故障特征提取流程
基于上述理论，本文提出基于ECYCBD的滚动轴承弱故障特征提取方法，适用于强背景噪声干扰下的滚动轴承故障特征提取问题，该方法的流程如图1所示，具体过程为：
图1 基于ECYCBD的滚动轴承弱故障特征提取方法流程
1)获取滚动轴承故障信号并利用信号的最大幅值(取绝对值)进行标准化处理(文中仿真信号和试验信号均进行了量纲一化处理)，计算该滚动轴承故障特征频率的理论值f。
2)基于故障特征频率的理论值确定滤波器长度和特征频率值的搜索范围，为保证CYCBD的分析效果和计算效率，本文将故障周期的长度记为L，进而将滤波器长度的选取范围确定为[0.75L,3L]，特征频率的选取范围确定为[f-5,f+5]。
3)以最大重加权峭度值为目标，利用鲸鱼优化算法优化选取CYCBD的参数。
4)借助最优参数组合，获得包含明显故障特征的解卷积信号。
5)利用Teager能量算子和快速傅里叶变换得到该解卷积信号的Teager能量谱。
6)通过识别Teager能量谱中的特征判断滚动轴承故障类型。
首先利用仿真信号对本文方法的有效性进行验证，仿真信号构建时模拟了滚动轴承内圈故障的特点。为模拟较为复杂的背景环境，该信号包含幅值调制、谐波信号和白噪声成分，可以表示为
(22)
式中：Ai为幅值调制成分；
B(t)为谐波成分；
s(t)为自由衰减信号；
n(t)为随机白噪声；
fr为转频，取42 Hz；
fm为谐波信号频率，取100 Hz；
fn为系统的固有频率，取3 200 Hz；
r为阻尼系数，取0.05；
T=1/185 s，即该信号的故障特征频率为185 Hz；
ti为第i个周期内由于滚动体滑移引起的延迟，ti=0.01T～0.02T。
仿真信号构造时，为突出背景成分的干扰，设定信号的信噪比为-10 dB，同时信号的采样频率和时长分别为32 768 Hz和2 s。该仿真信号的时域波形和包络谱如图2所示，时域波形中未能识别出周期性的轴承故障冲击成分，包络谱中仅能勉强识别故障特征频率及其2倍频成分，谱图中存在较强的干扰，说明对于该仿真信号，包络分析效果不够理想。本文的包络谱均是对信号直接进行希尔伯特变换和快速傅里叶变换获得。
(a) 时域波形
(b) 包络谱
下面利用本文提出的方法处理该仿真信号，在优化过程中，若WOA的种群数量和迭代次数设定过大，将导致整个算法的计算效率较低，若参数设置过小，会导致最后的寻优结果不理想。为兼顾算法的有效性和运行时间，通过多次尝试，本文将WOA中的参数设定为：种群数量20，最大迭代次数30，其余参数都选取默认值；
滤波器长度和特征频率的搜索范围确定为132～531和180～190 Hz。分析得到的迭代收敛结果如图3a所示，经过16次迭代，适应度函数取得最大值4.815，此时对应的鲸鱼位置为[525,185]。将此位置带入CYCBD进行计算，得到的最优解如图3b所示，该信号中的噪声成分得到明显去除，谐波成分的干扰得到有效抑制，且能够体现出一些冲击成分。该信号的Teager能量谱如图3c所示， 图中故障特征频率185 Hz及其倍频成分得到了准确提取，无关频率成分的干扰得到明显抑制。通过上述特征可以判断该轴承存在内圈故障，表明本文提出的方法能够解决强背景噪声带来的干扰，准确有效提取滚动轴承弱故障特征。
(a) 迭代收敛曲线
(b) 最优解卷积信号
(c) Teager能量谱
为进一步突出本文方法的分析效果，利用CYCBD和MCKD进行对比。
使用CYCBD对仿真信号进行分析，对参数进行预先确定，将滤波器长度和特征频率值分别设定为312和185 Hz，分析结果如图4所示。解卷积信号中噪声干扰仍较为明显，Teager能量谱提取出滚动轴承故障特征成分但幅值较低，尤其185和555 Hz处谱线幅值不突出，无关频率成分仍产生较明显干扰，分析效果欠佳。上述结果说明在分析过程中，随机确定的CYCBD参数并不准确。
(a) 解卷积信号
(b) Teager能量谱
利用经典的MCKD方法分析该仿真信号，并计算Teager能量谱。MCKD的参数依据本文优化过程得到的最优解[525,185]进行设置，其中滤波器长度取为525，冲击信号的周期T=172。MCKD方法对仿真信号的分析结果如图5所示，解卷积信号中故障冲击成分难以得到识别且噪声干扰较强，Teager能量谱虽然在特征频率附近出现了幅值增大现象，但未能提取出明显的故障特征频率成分。
(a) 解卷积信号
上述对比结果进一步突出本文算法在确定最优参数，准确提取故障特征上的有效性。
为进一步凸显本文算法在工程实际中的价值，对某型试验台获得的滚动轴承故障信号进行分析，试验台结构如图6所示，可用于模拟滚动轴承在实际工作中的常见故障类型和运行工况。为模拟实际设备结构中的复杂传递路径，确定测点时不考虑靠近故障轴承的轴承座测点，而是选取加载轴承端(测点1)和支承轴承端(测点2)作为本次试验的测点，测得的信号分别记为信号1和信号2，振动测量方向为径向。此时，轴承故障信号在向外传递时需要经过多个连接界面，还要受到主轴、正常轴承等结构的干扰，信号较为复杂，给故障特征提取带来一定难度。
图6 试验台
试验轴承型号为6010，内、外径分别为50，80 mm，球组节圆直径Dpw为65 mm，球径Dw为9 mm，球数Z为13，接触角α为0°。利用激光切割法在内圈沟道上割出一道宽0.2 mm、深0.2 mm的窄缝模拟内圈故障。试验过程中，试验台转速为3 000 r/min(转频fr=50 Hz)，径向和轴向加载分别为1和2 kN，采样频率为32 768 Hz，轴承内圈故障特征频率fi为370 Hz。
为更好提取故障源信号，本文在实测信号分析时采用多输入单输出的CYCBD方法，同时使用信号1和信号2作为CYCBD的输入进行解卷积计算。实测信号及其包络谱如图7所示：
由于复杂传递路径及噪声成分的干扰，时域波形中故障导致的冲击成分未能得到识别；
包络谱中无关干扰频率成分较为明显，故障特征未能得到有效提取：上述结果说明仅利用包络解调方法难以在强噪声干扰下获得理想效果。
(a) 信号1
(b) 信号2
利用本文方法分析该实测信号，结果如图8所示。分析时，滤波器长度的搜索范围为66～265，特征频率的搜索范围为365～375 Hz。通过WOA得到的迭代曲线如图8a所示，经过7次迭代可获得最大适应度值为4.742，此时对应的最优位置为[254，366]，此组参数也就是CYCBD计算时的最优参数组合。最优特征频率值366 Hz与理论值不一致可能是由于轴承尺寸存在制造误差以及试验过程中转速控制不准确。将优化结果代入CYCBD进行计算，得到的最优解如图8b所示，能够识别出周期性的冲击成分，噪声干扰得到明显抑制。Teager能量谱如图8c所示，图中存在转频成分、故障特征频率及其倍频成分，围绕故障特征频率成分的调制边频带都得到了准确提取， 故障特征较为丰富。通过上述分析结果，能够确定该轴承存在内圈故障，验证了本文方法在弱特征信号提取中的有效性。
(a) 迭代收敛曲线
(b) 最优解卷积信号
(c) Teager能量谱
与仿真分析类似，同样使用2种方法进行对比分析。
将CYCBD分析时的滤波器长度和特征频率值分别随机设定为321和370 Hz，将信号1和2共同作为CYCBD的输入， 分析结果如图9所示。解
(a) 解卷积信号
(b) Teager能量谱
积信号中难以识别故障冲击成分且噪声干扰成分的抑制效果不明显，Teager能量谱中也未能识别出任何有关故障的特征，说明随机确定CYCBD的分析参数难以获得较好的效果。
利用最优参数组合[254,366]确定MCKD的参数，MCKD的滤波器长度和冲击信号的周期分别取254和90，对2个信号进行分析。信号1的MCKD分析结果如图10所示，解卷积信号中未能识别出明显的冲击成分，且Teager能量谱仅能勉强确定故障特征频率的2倍频成分，故障特征并不丰富，噪声干扰较为明显，分析效果欠佳；
信号2的MCKD分析结果如图11所示，信号2的采集位置离故障轴承更远， 从其时域波形或Teager能量谱中均不能识别出有关故障的任何信息。
(a) 解卷积信号
(b) Teager能量谱
(a) 解卷积信号
(b) Teager能量谱
针对滚动轴承故障信号中噪声干扰成分强，故障冲击成分弱的问题，提出了基于ECYCBD的故障特征提取方法，其具有以下优点：
1)利用增强的CYCBD方法可以抑制无关成分干扰，有效提取故障特征，借助Teager能量谱可以进一步突出故障成分，提高解调效果；
2)以最大重加权峭度值为目标，借助WOA的参数寻优，保证了CYCBD参数选取的准确性，提高算法的分析效果。仿真和实测信号的对比分析结果凸显了本文方法在噪声干扰成分抑制和故障特征成分提取上的有效性。
