基于固有时间尺度分解与最大相关峭度解卷积的轴承故障特征提取

王银玲 尹显明 李华聪

(1. 西南科技大学工程技术中心 四川绵阳 621010； 2. 西北工业大学动力与能源学院 西安 710072)

滚动轴承广泛应用在旋转机械中，受负荷、摩擦、冲击等非线性因素影响会出现损伤。通常由于滚动轴承工作环境复杂，传感器采集到的信号含有大量的背景噪声，造成损伤信号的非线性、非平稳、冲击性特征被淹没。因此，如何通过对损伤过程中产生的非线性、非平稳随机信号进行分析，准确检测出设备的故障信号是滚动轴承故障检测研究的热点及难点[1-3]。

在轴承故障诊断中，损伤信号通常被看作是周期性冲击信号，其与系统响应的卷积即为传感器测量的输出信号。若不考虑噪声的影响，可以寻找一个逆滤波器，通过传感器输出的观测信号，利用解卷积恢复出原始的故障信号。Wiggin[4]结合熵值最小卷积滤波思想，提出最小熵解卷积方法(Minimum Entropy Deconvolution，MED)，在解卷积过程中使得峭度最大化，该方法又称为最大峭度解卷积。但MED只能处理单一冲击信号，并不适合旋转机械的周期性冲击信号的处理。

故障信号通过EMD分解[5]过程中存在端点效应以及模态混叠现象。为了解决这些问题，学者们提出了不同的解决方案。陈哲等[6]提出一种集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)[7]与多尺度排列熵(Multiscale Permutation Entropy，MPE)的舰船辐射噪声复杂度特征提取方法，解决复杂海洋环境中舰船辐射噪声的特征提取问题。耿读艳等[8]为了有效识别心冲击(Ballistocardiogram)信号，利用噪声自适应完备总体平均经验模态分解(Complete EEMD with Adaptive Noise，CEEMDAN)的自适应降噪优势，提出一种基于自适应噪声的CEEMDAN联合排列熵(PE)的BCG降噪方法，降噪后信号的幅频特性得到明显改善且信噪比较传统方法有明显提高。赵荣珍[9]提出一种基于经验小波变换、多尺度排列熵、GG聚类算法相结合的故障诊断方法。刘兴教[10]通过EEMD和最大相关峭度解卷积(Maximum Correlated Kurtosis Deconvolution，MCKD)提取柔性薄壁轴承故障特征。

由于以上方法不能产生适当的旋转分量的余量，且在信号分解过程中为了降低噪声及模态混叠，在模态分解过程引入了随机噪声，但是引入随机噪声后不能从根本上消除模态混叠，并可能会产生虚假分量。为了解决噪声及模态混叠对故障信号的不利影响，本文提出一种改进的固有时间尺度分解(Intrinsic Time-scale Decomposition, ITD)与最大相关峭度解卷积(MCKD)[11]相结合的故障信号特征提取方法。首先通过ITD将原始信号分解为一系列固有旋转分量和一个剩余的单调趋势分量，然后选择与原信号相似度高的旋转分量作为滤波信号的主要成分，利用粒子群算法优化MCKD的参数，使得峭度值最大，最后对信号进行包络解调，提取滚动轴承信号的故障特征。

1 改进的固有时间尺度分解

固有时间尺度分解(ITD)是由Frei等[12]提出的一种非平稳信号的时频分析方法。ITD可以实时提取信号的瞬时特征，将信号分解为一系列适当的旋转分量和一个趋势分量。对于t>0的待分解信号Xt，ITD的分解过程是定义一个运算算子L，从Xt提取出基线信号Lt，再通过Xt-Lt得到旋转分量Ht，用公式表示为：

Xt=LXt+(1-L)Xt=Lt+Ht

(1)

ITD分解[13]首先找出Xt上所有的局部极值点所对应的时间集合{τk,k=1,2,3…}，时间集合包含信号起始时间点。为了表示方便，用Xk和Lk代替X(τk)和L(τk)。假设Ht和Lt在时间区间[0,τk]上有定义，Xt中时间t∈[0,τk+2]，这样就可以在时间(τk,τk+1]上定义一个分段线性的基提取算子L：

t∈(τk,τk+1]

(2)

其中：

(1-α)Xk+1

(3)

式中，α∈(0,1)，通常取值为0.5。以这种方式构造的基线信号Lt，可以使其在极值之间保持单调性，并在极值点的包络内。由式(2)和式(3)得到基线点信号后，通过对其进行线性插值，得到基信号Lt。用待分解信号减去基信号Lt，就得到一个适当的旋转分量Ht。将提取出的基信号Lt作为原始信号，继续进行ITD分解，直到所分解的基信号小于某设定值时，ITD分解结束。最终ITD分解可以表示为：

H1+H2+…HP+LP

(4)

式中，Hi表示第i(i=1,2,…p)个适当的旋转分量，而Lp表示P次后的残余基线分量。

由于ITD分解过程中采用线性插值，线性插值只顾及其附近两点的影响，函数上两点之间的近似随着所近似函数的二阶导数增大而逐渐变差，为了解决这个问题，本文采用Akima插值方法。Akima插值方法是一种5点求导分段三次多项式插值算法，其5点求导法利用所求点与前后连续两点组成5个点，然后利用点之间延长线得到交点，最后得到所求点的导数[14]。

为了验证基于Akima插值ITD分解的有效性，给出一个合成的数字信号，其数学表达式为：

y(t)=y1(t)+y2(t)+0.07·random

(5)

其中random为随机信号：

y1(t)=[1+0.5cos(8πt)]·cos[200πt+

2cos(10πt)]

y2(t)=0.8sin(πt)sin(30πt)

图1为原始合成信号及其分量波形，图2为通过EMD分解得到前3个主要IMF分量的时域波形，图3为通过ITD分解后得到的特征波形。可以看到ITD分解方法要优于EMD分解方法，验证了该方法的有效性。

图1 原始合成信号及其分量Fig.1 Original composite signal and its components

图2 合成信号EMD分解后前3个IMF分量Fig.2 The first three IMF components ofcomposite signal after EMD decomposition

图3 合成信号ITD分解后前3个分量Fig.3 The first three components of the composite signal after ITD decomposition

2 最大相关峭度解卷积

传感器采集到的滚动轴承故障信号x(n)可以表示为：

x(n)=h(n)·y(n)+e(n)

(6)

式中，h(n)为系统的冲击响应，y(n)为系统故障信号，e(n)为噪声。在不考虑噪声的情况下，寻找最优滤波器f=[f1,f2,…fL]T，对系统采集的信号x(n)进行解卷积处理，可还原出故障信号y(n)，即

y(n)=f(n)·x(n)

(7)

熵表示系统内信息的混乱程度，熵值越大则系统越难预测。MED方法将熵值最小与卷积滤波思想进行了结合，对传感器采集到的原始信号进行解卷积，以反卷积信号x(n)峭度作为目标函数求解最优滤波器，以突出信号的冲击成分。但是滤波器选择不当将会提取出虚假成分，因此出现了以相关峭度(Correlation Kurtosis，CK)为最优指标的MCKD，相关峭度的计算式为：

(8)

(9)

式中L为滤波器长度，而f=[f1,f2,…fL]T。对目标函数求导，使得：

(10)

最终得到滤波器为：

(11)

其中：

(12)

r=[0,T,2T,…,mT]

(13)

(14)

从以上分析可知，通过MCKD将传感器采集信号x(n)进行解卷积滤波时，需要设定周期T、偏移量M、滤波器长度L以及迭代次数N。

3 试验验证

实验数据来源于某变速器轴承的外圈故障，变速器转速为n=923 r/min，轴承节径D=57.50 mm，滚动体直径d=14.22 mm，压力角α=0°，滚动体个数Z=7。外圈故障的特征频率满足以下关系：

(15)

通过上式可以计算出轴承外圈故障频率fr=17.43 Hz。

采集变速器轴承外圈故障信号，采样频率为40 kHz，采集数据长度为84×103。首先对信号做去均值处理，处理后信号的时域与频域特征如图4所示。

图4 原始故障信号时域及频域特征Fig.4 Time domain and frequency domain characteristics of the original fault signal

图5为原始信号EMD分解后一阶IMF分量的包络谱，由图5可以看出，故障信号的特征频率在35.4 Hz，105 Hz处都有较大的峰值，但是由于噪声的影响，不能判断是否含有冲击成分，也就是故障信号的特征包络谱。所以通过对原始信号的包络谱分析，不能有效提取故障特征频率。

图5 信号EMD分解后一阶IMF信号包络谱Fig.5 First-order IMF signal envelope spectrum after signal EMD decomposition

接下来需要去除系统的噪声，突出损伤冲击信号特征。对信号进行基于Akima插值ITD分解，分解后前4个分量的时域波形和其频谱如图6所示。

图6 ITD分解后各分量波形及频域特征Fig.6 Waveform and frequency domain characteristics of each component after ITD decomposition

ITD分解后各分量的能量比、峭度、峭度因子和相似性的统计特征如图7所示。由图7可知，第一阶次的能量比最大，是信号的主要成分；第一、二阶次的峭度最大，证明其变化最为剧烈；通过相似性分析可知第一阶次与原始信号最接近，故将第一阶分量作为提取信号。

图7 ITD分解后各分量的统计特征Fig.7 Statistical characteristics of each component after ITD decomposition

通过ITD分解得到信号的主要成分后，再对信号进行MCKD处理以得到信号的冲击成分。通过MCKD对信号进行解卷积时，需要设定周期T、偏移量M、滤波器长度L及迭代次数N。其中：周期T计算值并不一定最优，与实际还有偏差；而滤波器长度L直接影响着整个计算存储数据大小；偏移量M一般取1～7之间；迭代次数太小有可能不是最优解，太大影响计算效率，一般选取30。以上参数的选取直接影响着峭度，为了获得最优参数，本文以最大峭度为目标函数，采用粒子群算法(PSO)对上述参数进行优化，选取最优组合。最终选择滤波器尺寸L=100、最终迭代次数N=30、解卷积周期为T=30、偏移量M=5时取得最优解。图8比较了MCKD处理前后的信号谱峭度，图9为MCKD处理后信号的包络谱。图9清晰显示出故障频率为17.09，外圈故障频率及其倍频上出现峰值，能够清晰地诊断出设备的故障特征。

图8 MCKD处理前后信号谱峭度Fig.8 Signal spectral kurtosis before and after MCKD processing

图9 MCKD处理后信号的包络谱Fig.9 Envelope spectrum of signal processed by MCKD

4 结束语

对于非平稳性、非线性的轴承故障信号，本文提出了基于固有时间尺度分解和相关峭度解卷积相结合的特征提取方法。借助固有尺度分解，提取滚动轴承故障信号的主要特征，为了使得到的波形更加光滑平稳，在尺度分解中采用Akima插值方法代替线性插值。借助MCKD提取信号的冲击成分过程中，为了选择合适的参数，通过粒子群算法对参数进行优化，以最大峭度为目标函数，得到最优参数，最后对MCKD后信号进行包络谱分解，可以清晰地提取出故障特征频率。实验结果与理论一致，有利于对噪声环境下轴承的状态进行特征提取。