本发明属于机械设备信号处理领域，具体涉及一种滚动轴承微弱故障信号的特征提取方法及设备。
背景技术：
滚动轴承是旋转机械中应用最广泛最关键的零件之一，它的运行状态往往决定了整机的性能，任何轻微故障都会对机器运行稳定性产生很大影响。而滚动轴承在运转过程中又容易损坏，若能在故障早期阶段提取出轴承微弱的故障信号，对信号进行分析处理并及时给出准确的诊断结果，可以使得维修人员能够针对故障制定有效且合理的维修计划，从而延长机器寿命，极大地降低故障所带来的危害。
滚动轴承出现局部损伤时，每当滚动体和滚道在损伤部位接触，便会产生一个冲击并激起系统结构自身的高频振动，这个周期性的冲击就是故障诊断的重要特征。现有的分析方法主要利用离散小波或小波包、经验模态分解和希尔伯特-黄变换等方法，分解为从高频到低频的不同分量，利用相关系数、峭度或熵等多种特征对分量进行筛选，在选取分量上提取故障特征从而确定故障类别。然而在轴承的早期故障期间，特别是微弱的故障诊断中，轴承的局部缺陷和损伤很小，所引起的冲击振动非常微弱，在不同频带分量上提取到的故障特征容易被转频和机器运行的各种噪声掩盖，而且离散小波或小波包变换以及经验模态分解等频带划分方法离散间隔太大过于粗糙，也会影响故障特征的提取，在早期的微弱故障诊断中难以取得理想的效果。
技术实现要素：
发明目的：为了克服现有技术的不足，本发明提出一种基于改进的连续小波时频图的滚动轴承微弱故障信号特征提取方法，能够对早期轴承故障中微弱的周期性冲击成分进行有效提取，从而有助于快速准确的故障识别分类。
本发明的另一目的在于提供一种相应的滚动轴承微弱故障信号特征提取设备。
技术方案：为了实现以上目的，本发明所述的滚动轴承微弱故障信号特征提取方法包括以下步骤：
(1)获取滚动轴承的振动信号；
(2)对振动信号进行连续小波分解，得到连续小波时频图；
(3)对小波系数进行自相关运算，滤除噪声干扰；
(4)提取小波系数经过自相关运算得到的自相关函数的包络特征并进行包络谱分析，得到故障特征频率。
优选的，上述各步骤具体为：
(1)通过传感器和/或麦克风采集/获取滚动轴承的振动信号；
(2)利用Morlet连续小波变换对获取的滚动轴承振动信号进行分解，得到连续小波时频图；
(3)对时频图上每个频率对应的小波系数进行自相关运算，滤除噪声干扰并提取出周期性的故障冲击成分；
(4)在每个频率上对步骤(3)得到的自相关函数进行Hilbert变换，求得自相关函数的包络，对包络再进行傅里叶变换求得包络的功率谱，得到改进的小波时频图。
其中，上述步骤(2)中连续小波时频图的计算公式为：
其中，CWTx(a,b)表示连续小波变换时频图；x(t)为采集的振动信号；a为尺度因子，表示与频率相关的伸缩；b为平移因子；ψ(t)表示基小波，是基小波经过位移和伸缩产生的族函数，称为小波基函数，t表示时间。基小波选用Morlet小波，其形式为：
其中，σ为形状系数；f0为中心频率；i表示虚部。
步骤(3)中自相关运算的计算公式为：
其中，Rcc(τ)为小波系数c(t)的自相关函数，τ为延迟时间，T为信号观测时长。在小波时频图中，小波系数c(t)由周期性的故障冲击信号s(t)和非周期的噪声信号n(t)组成，将c(t)＝s(t)+n(t)带入公式(3)，计算可得：
将公式(4)展开可得：
Rcc(τ)＝Rss(τ)+Rsn(τ)+Rnn(τ) (5)
其中，Rss(τ)、Rsn(τ)和Rnn(τ)分别为周期故障冲击信号的自相关函数、周期信号与噪声信号的互相关函数、噪声的自相关函数。
步骤(4)中Hilbert变换的公式为：
其中，H[·]为Hilbert变换运算符，经过Hilbert变换，所有频率成分被相移90o，得到了新的时间信号，由此构造的新的解析信号R(τ)为：
其中j表示虚部，解析信号R(τ)的幅值就是Rcc(τ)信号的包络，计算公式如下：
对包络再进行快速傅里叶变换求得包络的功率谱，得到故障的特征频率。
一种滚动轴承微弱故障信号特征提取设备，包括信号采集装置和信号处理装置，其中，信号采集装置采用振动加速度传感器、麦克风、超声传感器中的一种或多种，用于获取滚动轴承的振动信号；信号处理装置包括存储器和处理器，用于对信号进行分析和处理，存储器存储有计算机程序，所述程序被处理器执行时能够实现以下步骤：
1)对振动信号进行连续小波分解，得到连续小波时频图；
2)对小波系数进行自相关运算，滤除噪声干扰；
3)提取小波系数经过自相关运算得到的自相关函数的包络特征并进行包络谱分析，得到故障特征频率。
进一步地，连续小波分解采用Morlet小波，其形式为：
其中，σ为形状系数；f0为中心频率；i表示虚部。
自相关运算的计算公式为：
其中，Rcc(τ)为小波系数c(t)的自相关函数，τ为延迟时间，T为信号观测时长；
自相关函数包络特征的提取采用Hilbert变换，其计算公式为：
其中，H[·]为Hilbert变换运算符。
有益效果：本发明提供的振动信号提取方法，能很好地凸显早期故障状态下振动信号中微弱的周期性故障冲击，尤其适用于滚动轴承早期微弱故障下含噪振动信号中调制特征的分析和特征提取，进而帮助机组开展故障诊断工作。具体而言，本发明的突出优势表现在：
1、本发明采用连续小波变换来分解信号，与离散小波和小波包变换相比，连续小波变换具有很好的时频分辨率和瞬态检测能力，在时域和频域同时具有良好的局部化性质，能够对信号的细节进行精确刻画，非常适合于对早期轴承故障中微弱的周期性冲击成分进行提取。利用连续小波变换分解信号得到的小波时频图能够准确地反映出滚动轴承早期故障中微弱的周期性冲击成分发生的时刻。
2、在连续小波时频图的基础上，本发明进一步通过自相关运算滤除了大量的非周期噪声信号，将对故障类别判断起决定性作用的微弱周期性故障冲击成分提取了出来，大大降低了环境噪声的干扰和影响，从而大幅提高故障检测的精度和速度。
附图说明
图1为本发明的故障信号特征提取方法的流程图；
图2为本发明的故障信号特征提取设备结构框图；
图3为根据本发明实施例的四种轴承的时域振动信号波形图；
图4为根据本发明实施例的外圈故障的Morlet连续小波时频图；
图5为根据本发明实施例的外圈故障小波时频图中4kHz和500Hz频率沿着时间轴对应的小波系数示意图；
图6为基于图5的小波系数经自相关运算后的得到的结果图；
图7为最终得到的滚动轴承在不同状态下的改进小波时频图；
图8为三种不同算法的故障诊断正确率对比图。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
参照图1和图2，本发明提供了一种基于改进的小波时频图的滚动轴承微弱故障信号特征提取方法，该方法包括以下步骤：步骤10、获取滚动轴承的振动信号；步骤20、对振动信号进行连续小波分解，得到连续小波时频图；步骤30、对小波系数进行自相关运算，滤除噪声干扰；步骤40、提取小波系数经过自相关运算得到的自相关函数的包络特征并进行包络谱分析，得到故障特征频率。不同的步骤由提取设备中相应的装置来执行，其中，步骤10由信号采集装置完成，步骤20-40由信号处理装置完成。本实施例中信号处理装置包括存储器和处理器，存储器存储有计算机程序，所述程序被处理器执行时能够实现步骤20-40所述的方法。以下详细描述具体过程。
步骤10、获取振动信号
本实施例中，利用振动加速度传感器以16kHz的采样频率采集滚动轴承在600r/min转速下的振动加速度信号，每个轴承的采集时长为1min。为了便于验证本发明方法的有效性，在本实施例中采集了四种类型的滚动轴承振动信号，分别是正常轴承、内圈故障轴承、外圈故障轴承、滚珠体故障轴承，三种故障轴承是使用雕刻机分别在相应轴承的内圈、外圈和滚动体上以人为方式制造轻微的划痕得到，以模拟滚动轴承的早期故障。对四种类型轴承所采集的信号在时域上的波形如图3所示。从图中可以看出，故障早期，振动信号中反映故障特征的周期性冲击成分非常微弱，被转频和噪声掩盖，在时域中无法直接观测到故障特征的存在，需要对振动信号进一步的处理。除了振动加速度传感器之外，还可以使用麦克风、超声传感器等来获取振动信号，获取的振动信号传输至信号处理装置进行分析和处理。
步骤20、连续小波变换
将每个轴承1min的信号划分为60个1s的等长段落x(t)，每一段作为一个信号样本。采用Morlet连续小波变换对采集的振动信号x(t)进行分解，得到连续小波变换时频图CWTx(a,b)，计算公式如下：
其中，a为尺度因子，表示与频率相关的伸缩；b为平移因子；ψ(t)为基小波；是基小波经过位移和伸缩产生的族函数，称为小波基函数；Morlet小波的定义为：
其中，σ为形状系数；f0为中心频率；i为虚部；t为时间。
本实施例利用MATLAB程序完成上述过程，使用MATLAB中的连续小波分解函数cwt，其调用格式为coefs＝cwt(x,scales,'wname')，x为输入的信号，wname为所使用的连续小波的名称，scales为尺度因子a的取值，直接指定scales的值计算得到的coefs是小波的尺度图而不是小波的时频图，尺度图的竖轴对应的是尺度因子a的不同取值，横轴是时间。为了得到小波的时频图，需要将小波的尺度因子a转化为对应的实际频率fa，尺度因子a和实际频率fa之间的对应关系如下：
fa＝fc*fs/a (3)
其中，fc为所用连续小波的中心频率；fs为信号的采样频率，为了使小波尺度图的频率范围在[0,fs/2]之间，尺度a的范围应该为[2*fc,inf]，其中inf表示无穷大。由公式3可以看出，为使转换后的频率序列为等差序列，尺度a的取值序列必须取为如下形式：
其中，totalscal是对信号进行小波变换时所用尺度序列的长度，c为一常数，从公式4中可以看出，c/totalscal对应的实际频率为fs/2，则c的计算公式如下：
c＝2*fc*totalscal (5)
将公式5带入公式4中就可以得到所需的尺度序列。本实施例在MATLAB中先使用函数centfrq得到Morlet小波的中心频率fc，其调用格式为wcf＝centfrq(‘morl’)，wcf为Morlet小波的中心频率，指定小波的尺度序列长度totalscal为512，连续小波变换函数cwt中的尺度序列scales可以由公式4计算得到，最后调用MATLAB中的scal2frq计算出小波尺度a对应的实际频率fa，调用格式为f＝scal2frq(scales,'wname',1/fs)，其中wname为所用连续小波的名称，fs为信号采样频率。
经过上述处理后得到了连续小波的时频图，横轴为时间，竖轴对应的是不同的频率，采样率为16kHz时，信号的频率在[0,8kHz]之间，被等分成了512个频段，频率的分辨率为15.625Hz。外圈故障的小波时频图如图4所示，在小波时频图上可以清晰的看到周期性故障冲击发生的时刻。
步骤30、自相关运算
然而，从图4中还可以看出，故障发生的早期，虽然Morlet小波从振动信号中提取出了周期性的故障冲击，但与低频部分的转频和机器运行的各种噪声相比，能量非常微弱。图5给出了图4的小波时频图中4kHz和500Hz频率沿着时间轴对应的小波系数。从图5中可以看出，4kHz频率上，Morlet小波提取出了振动信号中周期性的故障冲击成分，而500Hz左右的低频部分则主要包含轴承的转频以及各种非周期的噪声，与低频部分的噪声相比，周期性故障冲击的能量太小，振动信号中有助于故障判断的故障特征成分会被噪声所掩盖，从而影响故障类别的准确判断。为了解决这一问题，需要通过进一步的处理对微弱的故障特征进行加强。
本实施例中通过对时频图上每个频率对应的小波系数进行自相关运算，滤除噪声干扰并提取出周期性的故障冲击成分。振动信号x(t)经过Morlet连续小波变换之后得到了512×16000大小的时频图coefs(f,t)，竖轴代表频率，横轴代表时间，在MATLAB中是以矩阵形式存在的，矩阵的每一行就是每个频率所对应的小波系数coefs(fi,t)＝c(t)，fi代表第i个频率，小波系数c(t)的自相关函数计算公式如下：
其中，Rcc(τ)为信号c(t)的自相关函数，τ为延迟时间，T为信号的观测时长。在小波时频图中，当某个频率存在周期性的故障冲击时，沿着时间轴，其对应的小波系数c(t)应当由周期性的故障冲击信号s(t)和非周期的噪声信号n(t)组成，将c(t)＝s(t)+n(t)带入公式6，计算可得：
将公式7展开可得：
Rcc(τ)＝Rss(τ)+Rsn(τ)+Rnn(τ) (8)
其中，Rss(τ)、Rsn(τ)和Rnn(τ)分别为周期故障冲击信号的自相关函数、周期信号与噪声信号的互相关函数、噪声的自相关函数。由相关函数性质可知，经过自相关运算后，Rss(τ)的周期性与原信号s(t)一致，但s(t)与噪声信号n(t)的互相关函数Rsn(τ)和噪声信号的自相关函数Rnn(τ)均趋近于0，经过自相关运算消除了噪声的干扰，提取出了振动信号中的周期故障成分。
4kHz和500Hz频率对应的小波系数经自相关运算后的结果如图6所示，自相关运算主要是对周期信号的周期进行提取，与原信号本身能量大小无关，从图中可以看出，虽然低频部分的非周期噪声信号能量远大于高频部分周期性的故障冲击成分，但经自相关运算之后，噪声的自相关函数值变得很低，远小于故障冲击的自相关函数，在Morlet小波时频图的基础上，对每个频率对应小波系数进行自相关运算，大量的非周期噪声被滤除，微弱的周期性故障冲击成分被提取出来，信噪比得到了提升。
本实施例利用MATLAB程序完成上述过程，使用循环语句每次取出时频图矩阵的一行，先对这一行数据进行零均值化处理，调用MATLAB中的求均值函数mean，调用格式为ct＝ct-mean(ct)，其中ct为某个频率对应的小波系数，然后调用MATLAB中的自相关函数autocorr，调用格式为ac＝autocorr(ct,numLags)，其中numLags是自相关函数的计算长度，本实施例中取numLags的大小为信号ct的长度减去1。
步骤40、包络谱分析
自相关函数提取出了故障冲击发生的周期，采用Hilbert变换求得自相关函数Rcc(τ)的包络，对包络再进行傅里叶变换求得包络的功率谱，在包络的功率谱上可以得到故障冲击的频率。经过了自相关运算之后，小波时频图的每个频率对应的是小波系数的自相关函数Rcc(τ)，对Rcc(τ)进行Hilbert变换求取包络波形，计算公式如下：
其中，H[·]为Hilbert变换运算符。经过Hilbert变换，所有频率成分被相移90o，得到了新的时间信号，由此构造的新的解析信号R(τ)为：
其中j表示虚部，解析信号R(τ)的幅值就是Rcc(τ)信号的包络，计算公式如下：
对包络再进行快速傅里叶变换求得包络的功率谱，就可以得到改进的小波时频图。
本实施例利用MATLAB程序完成上述过程。调用MATLAB中的hilbert函数，调用格式为y＝hilbert(ac)，其中ac为小波系数的自相关函数，y为hilbert变换构造的解析信号，解析信号y的幅值就是ac的包络；调用MATLAB中求模的函数abs计算解析信号y的幅值，调用格式为yabs＝abs(y)；对幅值yabs进行零均值化处理，调用格式为yabs＝yabs-mean(yabs)；最后调用MATLAB中的功率密度函数psd计算yabs的功率谱，调用格式为ypsd＝psd(ypsd,nfft)，其中nfft是快速傅里叶变换的点数，本实施例取nfft＝16384。
图7为最终得到的滚动轴承不同状态下的改进小波时频图，从图中可以看出，故障轴承与正常轴承之间的差异得到很好的凸显，因为经过自相关运算后，低频部分的噪声被滤除，周期性的故障冲击成分被提取出来，再进行Hilbert包络解调，可以在时频图上清晰的看到故障特征频率，从而可以对早期的微弱故障进行非常准确的判断。
最后，本实施例通过试验对几种不同的故障信号特征提取方法在检测轴承微弱故障上的准确性进行比较。在信号时频分析领域常用的有两种方法：短时傅里叶变换(short-time fourier transform，STFT)和小波变换，通过短时傅里叶变换后得到STFT时频图；对小波变换，本发明进行了多种方法的研究和尝试，最后选用了两种不同于现有技术的小波时频图，一是在小波变换的基础上直接进行Hilbert包络解调得到小波时频图，二是如上所述小波变换后进行自相关运算再进行Hilbert解调得到改进的小波时频图。将上述STFT时频图和本发明提供的两种不同于现有技术的小波时频图进行对比试验，将时频谱作为特征输入，利用训练集在分类模型上进行训练后，使用测试集进行验证。图8示出了STFT时频谱、小波+Hilbert以及改进小波时频谱作为特征输入时，在训练集上同样进行1500次迭代时的准确率。从图中可以看出，采用改进小波时频图作为特征输入，分类模型的训练速度非常快且准确率更高。300次迭代，模型的准确率已经达到80％左右，600次迭代之后准确率在90％以上，1000次迭代以内准确率已经接近100％。而300次迭代时，小波+Hilbert和STFT时频图的准确率均低于60％；训练结束时小波+Hilbert变换的准确率最多能达到92％左右，而STFT时频图的准确率只有88％左右。使用改进的小波时频图故障信号特征提取方法能够在滚动轴承产生微弱故障的早期对故障特征进行准确的提取，使得滚动轴承微弱故障检测的准确率有了很大的提升，具有很好的实际应用价值。