基于不完全Cholesky分解相关熵双谱的轴承故障诊断

李 辉, 郝如江

(1.天津职业技术师范大学 机械工程学院，天津 300222; 2.石家庄铁道大学 机械工程学院，石家庄 050043)

滚动轴承是机械传动中最常用的零件，滚动轴承的健康状态是保证机械传动正常运行的重要因素，因此，检测滚动轴承的健康状态至关重要，基于振动信号的滚动轴承故障检测与诊断是最常用的方法之一[1-2]。当滚动轴承内圈、外圈或滚动体产生表面点蚀、疲劳脱落、局部裂纹等故障时，由于故障部位对与其接触的轴承其它部件产生冲击作用，从而产生周期性的瞬态脉冲冲击现象，这种瞬态冲击会激励机电设备的固有振动，并与轴承的回转频率产生调制，使得振动信号产生复杂的幅值、相位调制现象，不仅具有非线性、非高斯信号特征，而且往往包含大量的背景噪声，将轴承微弱故障特征淹没，这些因素加大了轴承故障诊断的难度。轴承故障振动信号的非线性、非高斯特征使传统基于信号二阶统计量的方法，如相关分析、功率谱分析等，性能衰退，甚至失效[3-4]。因此，针对强背景噪声下滚动轴承故障特征难以有效提取的问题，许多学者进行了研究。王宏超等[5]对二阶循环统计量的谱相关进行了改进，提出一种新的时频分析方法。鄢小安等[6]提出了一种基于改进奇异谱分解的形态学解调方法，并用于轴承故障诊断。胥永刚等[7]提出了改进奇异谱分解方法，并成功应用于轴承故障诊断。胡超凡等[8]提出了基于张量分解的多维降噪技术，用于滚动轴承复合故障多维降噪以及特征提取。在众多的机械设备故障诊断方法中，基于信号高阶统计量的信号处理方法，由于在非线性、非高斯信号处理方面有着独特优势，得到了更广泛的研究和应用，如双谱[9-11]、切片双谱[12]、倒双谱[13]、循环双谱[14-16]、高阶谱[17]分析等，已广泛用于机械设备的故障诊断，取得了良好效果。

轴承故障振动信号呈周期性瞬态冲击和复杂调制的特点，具有明显的非线性、非平稳、非高斯特点。双谱是分析高斯信号的有效工具，理论上高斯噪声的双谱为零, 但非高斯噪声的双谱并不为零，传统的双谱不能有效消除非高斯噪声的干扰，因此会造成双谱分析性能衰退，甚至失效。相关熵是处理高斯、非高斯噪声的有效方法，已在雷达和通信信号检测、信号滤波、波达方向估计和时延估计等方面得到应用，取得了良好效果[18-20]。尽管相关熵方法在通信领域的应用已经展开，但在机电设备故障诊断领域的应用才刚刚起步[21-22]。本文针对传统双谱难以有效处理非高斯噪声干扰以及相关熵计算量大的问题，提出了一种基于相关熵和双谱的轴承故障诊断方法，该方法首先利用核函数和不完全Cholesky分解(incomplete Cholesky decomposition, ICD)算法计算核矩阵的低秩分解—下三角矩阵G；其次，利用Gini指数(Gini index，GI)选取下三角矩阵G的主分量，由下三角矩阵G的主分量计算信号的核矩阵；根据核矩阵计算信号的相关熵；最后，计算相关熵的双谱，根据相关熵的双谱特征识别轴承故障，并利用仿真信号和电机轴承故障试验信号验证了该方法的有效性和可靠性。

1 相关熵和ICD方法简介

1.1 相关熵

对于实信号x(t)，时变自相关熵可定义为

(1)

式中:Ε(·)为期望算子;τ为时延;κσ(·)通常采用高斯核函数。κσ(·)表达式为

(2)

(3)

式中，m为时延的采样点数，m=0,1,2,…,N-1。

1.2 不完全Cholesky矩阵分解

(4)

式(4)称为矩阵κ的不完全Cholesky分解。当d<

1.3 基于Gini指数的矩阵降维

GI最初用于判断收入分配公平程度的指标，是一个比例数值，在0和1之间，GI越小，表示信号的幅值越均衡，反之，表示信号的幅值变化越大[24]。因此，可利用GI有效提取轴承故障振动信号中周期性的脉冲冲击。

向量z的GI定义为

(5)

z(1)≤z(2)≤…≤z(N)

(6)

(7)

(8)

2 基于相关熵的双谱估计

设{x(n)}是零均值的平稳随机过程，其三阶累积量Rxx(τ1,τ2)可表示为

Rxx(τ1,τ2)=E[x(n)x(n+τ1)x(n+τ2)]

(9)

式中，τ1、τ2为时延的采样点数。

信号{x(n)}的双谱可表示为三阶累积量Rxx(τ1,τ2)的傅里叶变换，即

(10)

式中，ω为圆频率，且|ω1|≤π,|ω2|≤π, |ω1+ω2|≤π。

(11)

(12)

从式(12)可以看出基于相关熵的双谱是圆频率ω1和ω2的函数。

3 基于相关熵双谱的故障诊断步骤

在实际工程应用中，根据以下步骤进行基于相关熵双谱的轴承故障诊断：

(5) 根据双谱的频谱结构识别轴承故障。

基于相关熵双谱的轴承故障诊断流程图如图1所示。

图1 相关熵双谱故障诊断流程图Fig.1 Flowchart of correntropy bispectrum based fault diagnosis

4 相关熵双谱轴承故障诊断仿真

滚动轴承内圈、外圈或滚动体产生表面点蚀、疲劳脱落、局部裂纹等是滚动轴承的典型故障。滚动轴承点蚀故障振动响应信号模型可以表示为[25]

sin[2πfbi(t-ti)+θbi]

(13)

x(t)=x1(t)+n1(t)+n2(t)

(14)

式中:Ai为瞬态脉冲幅值;Ci为阻尼衰减因子;ti为冲击持续的时间;θbi为初始相位;fbi电动机系统的共振频率；n1(t)为零均值高斯噪声;n2(t)为脉冲噪声。

函数Θ(t-ti)用来指定冲击发生的时间，可用下式定义

(15)

利用信号仿真x(t)验证相关熵的降噪性能和本文提出方法的有效性。设轴承外圈故障特征频率为fouter=120 Hz，电动机系统的固有振动频率fb=3 000 Hz，瞬态冲击振幅Ai=4.5，采样频率fs=20 000 Hz，信号采样点数n=4 000。

图2(a)和图2(b)分别为仿真信号x1(t)的时域波形和FFT，可以清晰看到信号x1(t)在频域内围绕轴承系统共振频率fb形成边频带簇，边频带的间隔为轴承外圈故障特征频率120 Hz。

(a) 时域图

(b) 频域图图2 仿真信号x1(t)及其FFTFig.2 Simulative signal x1(t) and its FFT

当轴承外圈发生故障时，在外圈故障双谱平面内，会产生以系统共振频率为中心、轴承外圈故障特征频率为间隔的谱峰簇。图3为仿真信号x1(t)相关熵(σ=3)的双谱，从图3的双谱平面内可以清晰看到：围绕系统共振频率fb，在(±fb,0)、(0,±fb)、(fb,-fb)以及(-fb,fb)等位置附近存在显著的谱峰簇，相邻谱峰之间的间隔为轴承外圈故障特征频率fouter=120 Hz，这种频谱特征表明了轴承外圈的故障特征，与理论分析一致。

图3 仿真信号x1(t)的双谱Fig.3 Bi-spectrum of simulative signal x1(t)

在仿真信号x1(t)中加入零均值高斯噪声n1(t)，形成信噪比为SNR=-10 dB的染噪信号，之后再随机加入几个幅值不同的脉冲信号，以模拟非高斯脉冲噪声，图4(a)和图4(b)分别为仿真信号x2(t)的时域波形和FFT，由于仿真信号x2(t)完全被强高斯和非高斯噪声淹没，因此从图4(a)已完全看不出信号幅值的变化规律。在图4(b)中，除在系统共振频率3 000 Hz位置出现比较大的峰值外，已无法识别在共振频率两侧的边频带簇。

(a) 时域图

(b) 频域图图4 仿真信号x2(t)及其FFTFig.4 Simulative signal x2(t) and its FFT

图5 矩阵各列的Gini指数Fig.5 Gini index of each column for matrix

图6为信号x2(t)的相关熵 (σ=3)。从图6可以看出：仿真信号x2(t)相关熵的时域波形呈现出比较明显的周期性幅值调制特征，在图4(a)中，不但幅值较大的非高斯脉冲噪声已被完全消除，而且强高斯噪声也得到了很好地抑制，但相关熵仍然保留了信号中的周期瞬态冲击成分。

图6 仿真信号x2(t)的相关熵Fig.6 Correntropy of simulative signal x2(t)

图7和图8为仿真信号x2(t)相关熵的双谱。从图7和图8可以看出，在强高斯噪声和非高斯噪声的影响下，基于相关熵的双谱仍然能够准确识别轴承外圈的故障特征，在(±fb,0)、(0,±fb)、(fb,-fb)以及(-fb,fb)等位置附近存在显著的谱峰簇，这些独立的谱峰仍然清晰可辨，相邻谱峰之间的间隔为轴承外圈故障特征频率fouter=120 Hz，由此可见，在强高斯和非高斯噪声干扰下，基于相关熵的双谱仍然能准确提取轴承外圈的故障特征。

图7 仿真信号x2(t)相关熵的双谱Fig.7 Correntropy based bi-spectrum of signal x2(t)

图8 仿真信号x2(t)相关熵的双谱Fig.8 Correntropy based bi-spectrum of signal x2(t)

为彰显相关熵的降噪能力，将基于相关熵的双谱与传统双谱和文献[10]中基于小波阈值降噪双谱进行了对比。图9和图10给出了利用直接法计算的仿真信号x2(t)的双谱，图11和图12为采用小波阈值降噪后信号的双谱，采用5层sym8小波分解。从图9、图10和图11、图12可以看出，仿真信号x2(t)的传统双谱和小波阈值降噪双谱已完全被强高斯和非高斯噪声淹没，完全不能分辨轴承故障特征频率。对比图8和图10、图12，图8中基于相关熵的双谱噪声方差很小，而传统双谱(图10)和基于小波阈值降噪双谱(图12)的噪声方差仍然很大，轴承故障频率已完全被强噪声淹没，不能有效识别。由此可见，尽管双谱对高斯噪声具有一定的抑制能力，但当信号中含有很强的高斯噪声和非高斯噪声时，双谱也难以有效提取轴承外圈的故障特征。而小波阈值降噪方法，在处理强高斯噪声和非高斯噪声干扰时，也难以取得理想的效果。

图9 仿真信号x2(t)的双谱Fig.9 Bi-spectrum of simulative signal x2(t)

图10 仿真信号x2(t)的双谱Fig.10 Bi-spectrum of simulative signal x2(t)

图11 仿真信号x2(t)小波阈值降噪双谱Fig.11 Wavelet de-noising based bi-spectrum of x2(t)

图12 仿真信号x2(t)小波阈值降噪双谱Fig.12 Wavelet de-noising based bi-spectrum of x2(t)

通过上述仿真信号可知：由于相关熵能有效抑制高斯噪声和脉冲噪声，因此基于相关熵的双谱分析具有从强高斯噪声和非高斯噪声背景中提取轴承故障特征的能力，基于相关熵的双谱分析为高斯、非高斯噪声的处理提供了一种崭新的鲁棒性解决方法。

5 电机轴承故障诊断

采用美国凯斯西储大学(Case Western Reserve University，CWRU)轴承数据中心的网站上公布的试验数据[26]，CWRU轴承数据集采集系统如图13所示。试验轴承型号为：深沟球轴承6205-2RS JEM SKF，采样频率fs=12 000 Hz，电机负载为空载，电机转速1 797 r/min(fr=29.95 Hz)，轴承内圈故障振动信号数据记录号为105DE，轴承外圈故障振动信号数据记录号为130DE。计算得到滚动轴承内圈、外圈故障特征频率如下：

finner=162.185 Hz，fouter=107.365 Hz

图13 CWRU轴承数据采集系统Fig.13 Bearing data collection system of CWRU

5.1 轴承内圈故障诊断

图14为滚动轴承内圈故障振动信号及其FFT。在图14(a)中存在明显的幅值调制现象，但根据时域波形还不能识别轴承故障。在图14(b)中，频谱幅值最大值处的频率fb=3 585 Hz，反映了电机系统的共振频率，但在fb的周围没有出现清晰的以轴承内圈故障特征频率为间隔的边频带簇，因此根据频域图不能有效识别轴承的故障类型。

(a) 时域图

(b) 频域图图14 轴承内圈故障振动信号的时域图和频谱图Fig.14 Bearing inner race fault vibration signal and its FFT

图15 矩阵各列的Gini指数Fig.15 Gini index of each column for matrix

图16 轴承内圈故障振动信号的相关熵Fig.16 Correntropy of bearing inner race fault vibration signal

图17和图18为轴承内圈故障振动信号相关熵的双谱，从图17和图18可以看出，双谱图中的谱峰呈离散点状分布，几乎没有受到噪声的影响。在低频范围，沿f1=0和f2=0轴的两侧，约在2 720 Hz处存在明显的谱峰簇，谱峰的频率等于轴承内圈故障特征频率finner，谱峰之间的间隔为轴承内圈故障特征频率。在高频范围，在f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz、f1=±3 585 Hz等10条直线的交点附近，存在明显的谱峰簇且谱峰之间的间隔为轴承内圈故障特征频率，这种频谱结构表明了轴承的内圈故障特征，同时也说明电动机系统在1 310 Hz、2 720 Hz以及3 585 Hz频率处具有较大的振动能量，与图14(b)FFT频谱图的结果一致。

图17 轴承内圈故障振动信号相关熵的双谱Fig.17 Correntropy based bi-spectrum of bearing inner race fault vibration signal

图18 轴承内圈故障信号相关熵的双谱Fig.18 Correntropy based bi-spectrum of bearing inner race fault vibration signal

图19和图20给出了利用直接法计算的轴承内圈故障振动信号的双谱，图21和图22为采用小波阈值降噪后信号的双谱。在图19、图20中，因噪声的影响，可模糊看出f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz、f1=±3 585 Hz等10条频谱直线，但却不能识别轴承内圈故障特征频率。在图21、图22中，因小波滤波可视为带通滤波，经过小波阈值滤波后，去掉了信号的高频部分，只保留了信号中的中、低频部分，可模糊看出f1=±1 310 Hz、f2=±1 310 Hz、f1=±2 720 Hz、f2=±2 720 Hz等8条频谱直线，同样也不能识别轴承内圈故障特征频率。在上述频谱直线的交点处出现谱峰，表明轴承内圈故障振动信号具有很强的非线性耦合，但轴承内圈故障振动信号的双谱和小波阈值降噪后双谱已完全被强噪声淹没，因此，难以准确分辨轴承内圈故障特征频率及其高次谐波。

图19 轴承内圈故障信号的传统双谱Fig.19 Bi-spectrum of bearing inner race fault vibration signal

图20 轴承内圈故障信号的传统双谱Fig.20 Bi-spectrum of bearing inner race fault vibration signal

图21 轴承内圈故障信号小波阈值降噪后的双谱Fig.21 Wavelet de-noising based bi-spectrum of bearing inner race fault vibration signal

图22 轴承内圈故障信号小波阈值降噪后的双谱Fig.22 Wavelet de-noising based bi-spectrum of bearing inner race fault vibration signal

5.2 轴承外圈故障诊断

图23为滚动轴承外圈故障时采集的振动信号及其FFT。在图23(a)中存在明显的幅值调制现象，在图23(b)中的频率等于3 445 Hz处存在显著的峰值，但根据时域波形及其FFT还不能准确识别轴承故障类型。

(a) 时域图

(b) 频域图图23 轴承外圈故障振动信号的时域图和频谱图Fig.23 Bearing outer race fault vibration signal and its FFT

图24 矩阵各列的Gini指数Fig.24 Gini index of each column for matrix

图25 轴承外圈故障振动信号的相关熵Fig.25 Correntropy of bearing outer race fault vibration signal

图26和图27分别为轴承外圈故障振动信号相关熵(σ=3)的轮廓图和三维图。从图26和图27可以清楚地看到：在f1=±3 445 Hz、f2=±3 445 Hz等4条频谱直线的交点附近，存在6个明显的谱峰簇，且谱峰之间的间隔为轴承外圈故障特征频率，这种频谱结构表明了轴承的外圈故障特征，同时也说明电动机系统在3 445 Hz频率处具有较大的振动能量，与图23(b)FFT频谱图的结果一致。

图26 轴承外圈故障振动信号相关熵的双谱Fig.26 Correntropy based bi-spectrum of bearing outer race fault vibration signal

图27 轴承外圈故障信号相关熵的双谱Fig.27 Correntropy based bi-spectrum of bearing outer race fault vibration signal

为了与本文提出的方法进行对比，图28和图29给出了利用直接法计算的轴承外圈故障振动信号的双谱，图30和图31为采用小波阈值降噪后信号的双谱。在图28和图29中，因噪声的影响，可模糊看出f1=±3 445 Hz、f2=±3 445 Hz等频谱线，但却不能清晰识别以轴承外圈故障特征频率为间隔的边频带簇。在图30和图31中，小波阈值降噪后信号的双谱也难以正确分辨轴承外圈故障特征频率。

图28 轴承外圈故障信号的传统双谱Fig.28 Bi-spectrum of bearing outer race fault vibration signal

图29 轴承外圈故障信号的传统双谱Fig.29 Bi-spectrum of bearing outer race fault vibration signal

图30 轴承外圈故障信号小波阈值降噪后的双谱Fig.30 Wavelet de-noising based bi-spectrum of bearing outer race fault vibration signal

图31 轴承外圈故障信号小波阈值降噪后的双谱Fig.31 Wavelet de-noising based bi-spectrum of bearing outer race fault vibration signal

通过上述分析可以看出：基于相关熵的双谱具有很强的从强噪声背景中提取轴承故障特征的能力，在双谱图中能清晰地刻画轴承内圈、外圈故障特征，提高了轴承故障诊断的可靠性和准确性，为从强高斯和非高斯噪声环境中提取轴承故障特征的有效方法。

6 结 论

传统双谱具有较强的抑制高斯噪声的能力，但却不能有效抑制非高斯噪声的干扰。相关熵不仅能有效抑制高斯噪声，而且能有效抑制非高斯噪声，因此，相关熵为高斯、非高斯噪声的处理提供了一种崭新的鲁棒性方法。综合利用了相关熵和双谱的优点，提出了基于相关熵双谱分析的轴承故障诊断方法。首先利用不完全Cholesky分解算法，计算核矩阵的低秩下三角矩阵，再利用Gini指数选择下三角矩阵的主分量，不仅降低了运算量，提高了核矩阵的计算速度，而且突出了轴承故障的瞬态冲击特征，提高了轴承故障诊断的准确性和可靠性。仿真和试验结果表明：基于相关熵和双谱分析技术，能有效提取淹没在强高斯和非高斯噪声环境中的微弱信号，提高了信噪比，为轴承故障诊断的有效方法。