本发明属于机械设备信号处理领域，具体涉及一种滚动轴承早期微弱故障的诊断方法。
背景技术：
滚动轴承是旋转机械中应用最广泛最关键的零件之一，它的运行状态往往决定了整机的性能，任何轻微故障都会对设备运行稳定性产生很大影响，甚至会导致安全隐患并造成重大经济损失。若能在故障早期阶段提取出轴承微弱的故障信号，对信号进行分析处理并及时给出准确的诊断结果，可以使得维修人员能够针对故障制定有效且合理的维修计划，从而延长设备寿命，极大地降低故障所带来的危害。因此，针对滚动轴承早期微弱故障的诊断研究不仅具有重要的理论意义，更有很好的实际应用价值。
滚动轴承出现局部损伤时，每当滚动体和滚道在损伤部位接触，便会产生一个冲击并激起系统结构自身的高频振动，这个周期性的冲击就是故障诊断的重要特征。现有的滚动轴承故障诊断方法主要利用离散小波或小波包、经验模态分解和希尔伯特-黄变换等方法，将采集的振动信号分解为从高频到低频的不同频带分量，在每个分量上提取如相关系数、峭度、能量或熵等多种可以表征周期性故障冲击存在的特征，将不同分量上提取的特征整理成特征向量或特征矩阵输入bp(backpropagation，反向传播)神经网络、svm(supportvectormachine)支撑向量机等网络模型中进行训练，利用训练好的模型对故障进行诊断。
然而在轴承的早期故障期间，特别是微弱的故障诊断中，轴承的局部缺陷和损伤很小，所引起的冲击振动非常微弱，在不同频带分量上提取到的故障特征容易被转频和设备运行的各种噪声掩盖，而且离散小波或小波包变换以及经验模态分解等频带划分方法离散间隔太大过于粗糙，也会影响故障特征的提取，在早期的微弱故障诊断中难以取得理想的效果。
技术实现要素：
发明目的：为了克服现有技术的不足，本发明提出一种基于小波变换和深度学习的滚动轴承微弱故障诊断方法，能够对早期轴承故障中微弱的周期性冲击成分进行有效提取并且进行快速准确的故障识别分类。
技术方案：为了实现以上目的，本发明所述的基于小波变换和深度学习的滚动轴承微弱故障诊断方法包括以下步骤：
1)获取滚动轴承在不同状态下的振动信号，包括正常状态、外圈故障状态、内圈故障状态、滚珠体故障状态。
2)对采集的振动信号进行连续小波分解、去噪和解调处理，得到改进的小波时频图，具体包括：
21)利用连续morlet小波变换对采集的信号样本进行分解得到连续小波变换时频图；
22)对时频图上每个频率对应的小波系数进行自相关运算，得到小波系数的自相关函数，通过自相关运算滤除噪声干扰并提取出周期性的故障冲击成分；
23)对小波系数的自相关函数进行hilbert变换，求得自相关函数的包络，对包络再进行傅里叶变换求得包络的功率谱，得到改进后的小波时频图。
3)将改进的小波时频图作为分类模型的特征输入，训练分类模型。其中分类模型采用基于深度学习的神经网络模型，具体为卷积神经网络(convolutionneuralnetwork，cnn)模型。并且将改进的小波时频图分为训练集样本和验证集样本，使用训练集对cnn模型进行训练，使用验证集对网络模型进行调参，得到训练好的cnn模型。
4)使用训练好的分类模型对滚动轴承的故障类别进行判断。将待分类的振动信号按照步骤2的处理得到改进的小波时频图，输入到训练好的cnn模型中，cnn模型给出对应的故障类别。
有益效果：与现有技术相比，本发明具有以下有益效果：
1、本发明采用连续小波变换来分解信号，与离散小波和小波包变换相比，连续小波变换具有很好的时频分辨率和瞬态检测能力，在时域和频域同时具有良好的局部化性质，能够对信号的细节进行精确刻画，非常适合于对早期轴承故障中微弱的周期性冲击成分进行提取。利用连续小波变换分解信号得到的小波时频图能够准确地反映出滚动轴承早期故障中微弱的周期性冲击成分发生的时刻。
2、在连续小波时频图的基础上，本发明进一步通过自相关运算滤除了大量的非周期噪声信号，将对故障类别判断起决定性作用的微弱周期性故障冲击成分提取了出来，降低了环境噪声对故障检测的影响，使得深度网络模型在微弱故障上的检测精度和速度有了大幅度的提升。
3、本发明采用cnn模型，它是一种含有多个隐层的深度学习网络模型，能够通过训练实现自动对滚动轴承的微弱故障特征进行学习，不需要人为对故障特征进行选择；并且与bp神经网络和svm支撑向量机等浅层网络相比，cnn模型对复杂特征的学习表达能力更强，运算速度更快，对非线性模式识别有很强的分类能力，其自动化程度、准确性以及泛化能力都要更高。
4、本发明结合改进的小波时频图作为输入，使cnn模型能够在滚动轴承产生微弱故障的早期对故障类别进行非常准确的判断，在不同的轴承上使用时也能给出正确的判断结果，训练的速度也很快，具有非常好的实际应用价值。在测试集上的正确率达到了97％以上，具有很好的泛化性，使得微弱故障检测在实际生产中的应用成为了可能。
附图说明
图1为本发明方法的流程图；
图2为根据本发明实施例的正常、内圈故障、外圈故障、滚珠体故障的时域振动信号示意图；
图3为根据本发明实施例的训练集、验证集、测试集的构建示意图；
图4为根据本发明实施例的外圈故障的morlet连续小波时频图；
图5为根据本发明实施例的外圈故障小波时频图中4khz和500hz频率沿着时间轴对应的小波系数示意图；
图6为根据图5的小波系数经自相关运算后的得到的结果图；
图7为滚动轴承在不同状态下的改进小波时频图；
图8为本发明使用的cnn模型的结构图；
图9为三种不同算法的故障诊断正确率对比图。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
本发明提供了一种基于小波变换和深度学习的滚动轴承微弱故障诊断方法，如图1所示，该方法包括以下步骤：
步骤1、数据获取
准备200个型号为nsk-6304的滚动轴承，其中50个为正常轴承，50个为内圈故障轴承，50个为外圈故障轴承，50个为滚珠体故障轴承。三种故障轴承是使用雕刻机分别在相应轴承的内圈、外圈和滚动体上以人为方式制造轻微的划痕得到，以模拟滚动轴承的早期故障。利用振动加速度传感器以16khz的采样频率分别采集各种滚动轴承在600r/min转速下的振动加速度信号，采集的信号在时域上的波形如图2所示。从图中可以看出，故障早期，振动信号中反映故障特征的周期性冲击成分非常微弱，被转频和噪声掩盖，在时域中无法直接观测到故障特征的存在，需要对振动信号进一步的处理。
步骤2、样本划分
在处理之前，为了适应后续分类模型的训练，先将信号数据进行划分，分为训练集、测试集和验证集。参照图3，将每种类别对应的50个轴承中，选取30个作为训练轴承，余下的20个中9个作为验证集，11个作为测试集，每个轴承采集总长为1min的振动信号，将这1min的信号划分为60个1s的等长段落x(t)，每一段作为一个样本，每种类别的样本数和对应标签如表1所示。
表1不同故障类别的样本数及对应标签
步骤3、特征提取
分别对训练集样本、验证集样本和测试集样本进行特征提取，得到改进后的小波时频图作为分类模型的特征输入，主要包括如下步骤：
步骤31)采用morlet连续小波变换对训练集、验证集和测试集中采集的振动信号x(t)进行分解，得到连续小波变换时频图cwtx(a,b)，计算公式如下：
其中，a为尺度因子，表示与频率相关的伸缩；b为平移因子；ψ(t)为基小波；
其中，σ为形状系数；f0为中心频率；i为虚部；t为时间。
经过上述处理后得到了连续小波的时频图，竖轴对应的是不同的频率，采样率为16khz时，信号的频率在[0,8khz]之间，被等分成了512个频段，频率的分辨率为15.625hz，横轴为时间，外圈故障的小波时频图如图4所示，在小波时频图上可以清晰的看到周期性故障冲击发生的时刻。
步骤32)对时频图上每个频率对应的小波系数进行自相关运算，滤除噪声干扰并提取出周期性的故障冲击成分。样本集中的振动信号x(t)经过morlet连续小波变换之后得到了512×16000大小的时频图coefs(f,t)，竖轴代表频率，横轴代表时间，在程序中是以矩阵形式存在的，矩阵的每一行就是每个频率所对应的小波系数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，经过自相关运算消除了噪声的干扰，提取出了振动信号中的周期故障成分。
图5给出了图4的外圈故障小波时频图中4khz和500hz频率沿着时间轴对应的小波系数c(t)的值，从图中可以看出，4khz频率上，morlet小波提取出了振动信号中周期性的故障冲击成分，而500hz左右的低频部分则主要包含轴承的转频以及各种非周期的噪声，与低频部分的噪声相比，周期性故障冲击的能量太小，直接将小波时频图作为分类模型的特征输入，振动信号中有助于故障判断的故障特征成分会被噪声所掩盖，使得分类模型不能够对早期故障的类别进行准确判断。4khz和500hz频率对应的小波系数经自相关运算后的结果如图6所示，自相关运算主要是对周期信号的周期进行提取，与原信号本身能量大小无关，从图中可以看出，虽然低频部分的非周期噪声信号能量远大于高频部分周期性的故障冲击成分，但经自相关运算之后，噪声的自相关函数值变得很低，远小于故障冲击的自相关函数，在morlet小波时频图的基础上，对每个频率对应小波系数进行自相关运算，大量的非周期噪声被滤除，微弱的周期性故障冲击成分被提取出来，信噪比得到了提升。
步骤33)采用hilbert变换求得自相关函数rcc(τ)的包络，对包络再进行傅里叶变换求得包络的功率谱，得到改进的小波时频图。经过了自相关运算之后，小波时频图的每个频率对应的是小波系数的自相关函数rcc(τ)，对rcc(τ)进行hilbert变换求取包络波形，计算公式如下：
其中，h[·]为hilbert变换运算符。经过hilbert变换，所有频率成分被相移90o，得到了新的时间信号，由此构造的新的解析信号r(τ)为：
其中j表示虚部，解析信号r(τ)的幅值就是rcc(τ)信号的包络，计算公式如下：
对包络再进行快速傅里叶变换求得包络的功率谱，就可以得到改进的小波时频图。
自相关函数提取出了故障冲击发生的周期，经过hilbert和快速傅里叶变换之后，在包络的功率谱上可以得到故障冲击的频率，轴承工作在不同状态下故障冲击的频率是不同的，通过故障频率计算公式可知，本实施例使用的轴承转速为600r/min时，故障频率及其谐波分量主要集中在150hz以内，对每个频率的小波系数经过自相关和hilbert变换后得到的包络功率谱保留150hz以内的频段，则最终的改进小波时频图的大小为512×150。
图7为滚动轴承在不同状态下的改进小波时频图，从图中可以看出，经过自相关运算后，低频部分的噪声被滤除，周期性的故障冲击成分被提取出来，再进行hilbert包络解调，可以在时频图上清晰的看到故障特征频率，四种类别之间的差别非常明显。
步骤4、网络训练
将训练集样本和验证集样本提取的改进小波时频图作为分类模型的输入，使用训练集对模型进行训练，使用验证集对网络模型进行调参，得到训练好的分类模型。
分类模型选用卷积神经网络(cnn)模型，cnn是为处理二维特征图像识别而设计的一类包含多个隐藏层的深度学习网络模型，可直接将二维图像作为输入，所需预处理工作少，学习能力强，在手写字符识别、人脸识别等多个领域得到了广泛应用，但cnn在滚动轴承故障检测中的研究还比较少。本实施例中将cnn模型引入轴承早期微弱故障诊断中，将经过改进的小波时频图作为cnn模型的特征输入，轴承的故障类别作为模型输出对应的标签，训练cnn模型对故障类别进行判断。
本实施例使用的cnn模型的结构如图8所示。cnn网络由输入层、隐藏层、全连接层和输出层组成，其中，隐藏层由两个卷积层c1、c2和两个采样层s1、s2交替组成，卷积层c1和c2的卷积核个数分别为16和32，均为5×5的大小，激活函数为relu，可以使用0对边界进行填充，所以经过卷积运算卷积层图像大小不变；采样层采用maxpooling方式，区域大小为2×3，且区域不重叠；全连接层包含64个神经元和输出层一起构成softmax分类器，输出层为4个神经元，输出4种类别的判断结果。
cnn模型的训练主要包含数据的正向传播和误差的反向传播两个部分。首先，设定好模型训练的各个参数，初始化网络各层的权重w和偏置b；接着，进行模型训练的正向传播过程，将改进的小波时频图作为cnn模型的特征输入，依次经过卷积层、采样层和全连接层的处理后送入输出层，每一层的输出都是下一层的输入；在得到了模型的输出之后，进行反向传播运算，将模型的输出与期望输出之间进行比较，得到两者的误差，通过bp反向传播算法将误差分配到各层，对模型的权重和偏置进行调整，直至满足收敛条件，结束cnn模型的训练，得到训练好的cnn网络模型。本实施例中，cnn的正向传播过程如下：
(1)输入。改进的小波时频图，大小为512×150；
(2)卷积层c1。卷积层主要进行特征提取，c1一共包含16个5×5大小的卷积
核，不同卷积核用于提取不同特征，依次将这16个卷积核与输入的小波时频图进行卷积，得到的结果并不直接作为c1的输出，而是先经过一个激活函数再加上偏置项再作为卷积层c1提取的特征输出，激活函数一般取relu函数，计算过程如下：
y＝relu(wx+b)
其中，w为5×5大小的卷积核，可以看为一个滑动窗口，在512×150的特征图上按照从左到右，从上到下的顺序依次滑动，x为卷积核从特征图上取出的5×5大小的图像块，对于图像块，采用卷积核w进行卷积，偏置项为b，输出的y为卷积运算的结果，由于使用了0对边界进行填充，所以经过卷积运算后c1输出的特征图像大小不变，仍然为512×150的矩阵，一共有16个卷积核，所以最终c1的输出为16×512×150的矩阵；
(3)采样层s1。s1层也称为特征映射层，负责对c1层获得的特征进行子采样，降低特征图大小的同时增强对噪声的鲁棒性，本实施例中采样层采用maxpooling方式，区域大小为2×3，且区域不重叠，即在c1输出的512×150的特征图上，在每个2×3大小的区域上只选取最大的数作为这个区域的特征进行保留，经过了s1采样之后，特征图的大小变为了16×171×75的矩阵；
(3)余下的卷积层c2与采样层s2。这些层的原理与上述一致，将s1输出的特征图依次经过卷积层c2和采样层s2的作用，最终得到32×57×38大小的特征矩阵；
(4)全连接层f1。全连接层f1包含64个神经元，与采样层s2进行全连接，s2的输出包含32×57×38＝69312个特征，在与f1层相连之前，需要先把s2层的3维特征矩阵转化为1维的特征向量，每个特征都与f1层的一个神经元相连；
(5)输出层。输出层与f1进行全连接，输出层为4个神经元，输出4种类别的判断结果，(1,0,0,0)t对应正常状态，(0,1,0,0)t对应滚珠体故障，(0,0,1,0)t对应内圈故障，(0,0,0,1)t对应外圈故障。
本实施例的实验平台配置如下：计算机为ubutu16.04系统，gpu为nvidagtx-980ti，cnn网络模型使用tensorflow1.0和python3.5搭建，网络搭建的具体步骤如下：
(41)在tensorflow中使用truncated_normal函数和constant函数来分别初始化第一层卷积层的权重w和偏置b，调用格式为w＝tf.truncated_normal(shape,stddev)，truncated_normal函数的作用为生成一个shape大小，方差为stddev的截断正态分布，本实施例中shape＝[5,5,1,16]，stddev＝0.1；b＝tf.constant(0.1,shape＝shape)，constant函数的作用为生成一个值为0.1，大小为shape的常数，本实施例中shape＝[16]，其中tf为tensorflow的简称。
(42)调用tf.nn.conv2d函数和tf.nn.relu函数计算第一层卷积层的输出结果，首先使用tf.nn.conv2d函数计算输入图像经过第一层卷积后的值，调用格式为conv＝tf.nn.conv2d(x,w,strides,padding)，其中x为输入的特征图像，本实施例中为改进的小波时频图，w为卷积层的权重，strides指定了卷积核在每一维度滑动的步长，padding可以选择在卷积运算时是否进行填充，conv为输入经过卷积运算后的值，本实施例中，strides＝[1,1,1,1]，padding＝’same’可以使用0对边界进行填充；然后，调用tf.nn.relu函数计算conv加上偏置b后经过激励函数relu运算的值，调用格式为h_conv＝tf.nn.relu(conv+b)，h_conv就为第一个卷积层的输出。
(43)使用tf.nn.max_pool函数对卷积层的输出进行采样，调用格式为h_pool＝tf.nn.max_pool(h_conv,ksize,strides,padding)，其中ksize指定了采样范围，本实施例中设定ksize＝[1,2,3,1]，strides＝[1,2,3,1],padding＝’same’，h_pool就为第一个采样层的输出，第二个卷积层c2和采样层s2的处理与上述一致。
(44)在进行全连接运算之前，需要将2维的特征图转化为1维的，调用tf.reshape函数改变特征的维度，调用格式为h_flat＝tf.reshape(h_pool2,[-1,38*57*32])，其中h_pool2为第二个采样层的输出，经过两次卷积和采样之后，每个特征图的大小为38×57，第二个卷积层有32个卷积核，所以转化的1维特征向量h_flat的大小为38×57×32。
(45)进行第一个全连接运算，使用tf.truncated_normal函数初始化全连接的权重w_fc，shape＝[38*57*32,64]，输入层为1维特征向量的大小，隐藏层为64个神经元，调用tf.matmul函数计算输入特征和隐层神经元的积，调用格式为fc＝tf.matmul(h_flat,w_fc)，然后调用tf.nn.relu函数计算fc加上偏置b_c后经过激励函数relu运算的值，调用格式为h_fc＝tf.nn.relu(fc+b_c)，h_fc就为全连接层的输出。
(46)最后，计算输出层的输出，使用tf.truncated_normal函数初始化输出层的权重w_out，shape＝[64,4]，一共有4个输出，和全连接层一样计算出输入特征和输出层权重w_out的积，调用格式为y＝tf.matmul(h_fc,w_out)，在输出层使用softmax函数来计算输出的结果，调用格式为y_out＝tf.nn.softmax(y,b_out)，y_out就为cnn卷积网络的输出。
cnn模型的网络结构确定后，使用训练集和验证集对模型进行训练，训练的参数设定如下：损失函数采用tensorflow中的softmax_cross_entropy_with_logits函数，调用格式为loss＝tf.nn.softmax_cross_entropy_with_logits(labels,logits)，其中labels是cnn网络输出对应的标签，logits是cnn网络的输出y_out，loss为计算出的损失函数，每次训练的batch＝100，调用reduce_mean计算出整个batch的误差，调用格式为loss＝tf.reduce_mean(loss)，这里的loss就是网络训练所要减小的误差，训练迭代总次数epoch＝2500，学习率lr＝1e-4。
步骤5、故障诊断
将测试集样本提取的改进小波时频图输入到训练好的cnn模型中，cnn模型会给出对应的故障类别。
本实施例通过试验对几种故障检测方法的准确性进行对比。在信号时频分析领域常用的有两种方法：短时傅里叶变换(short-timefouriertransform，stft)和小波变换，通过短时傅里叶变换后得到stft时频图；对小波变换，本发明进行了多种方法的研究和尝试，最后选用了两种不同于现有技术的小波时频图，一是在小波变换的基础上直接进行hilbert包络解调得到小波时频图，二是如上所述小波变换后进行自相关运算再进行hilbert解调得到改进的小波时频图。将上述stft时频图和本发明提供的两种不同于现有技术的小波时频图进行对比试验，将时频谱作为特征输入，利用训练集在cnn模型上进行训练后，使用测试集进行验证。训练集上三种方法的结果如表2所示。
表2cnn模型在训练集上的正确率
stft使用的窗长为1024个采样点，窗移64个采样点，morlet连续小波的最大尺度a设为512，相同采样率下，两者的频率分辨率是一致的。从表2中可以看出，训练集上，stft时频图的正确率只有80％左右，与stft相比，使用小波变换进行故障诊断的正确率有所提升，直接使用hilbert对小波系数解调得到的时频图准确率能够达到92％，但与改进的小波时频图相比，仍然存在一定程度的误判，改进的小波时频图的正确率达到了99％以上，4种类别被很好的区分开来。
实际的工业应用中，训练好的模型需要在不同的轴承上进行检测，为了验证模型的泛化性，本发明挑选了训练集之外的轴承组成了测试集，其判断结果如表3所示。
表3cnn模型在测试集上的正确率
从表3中可以看出，在测试集上，stft时频图的正确率只有70％左右，小波+hilbert时频图的正确率下降到90％以下，而改进的小波时频图的正确率依旧保持在97％以上，说明使用改进的小波时频图训练的cnn模型在不同的轴承上仍然能够给出准确的判断，特别是实际工程中，小波+hilbert和stft训练出的模型在训练集之外的轴承上进行检测时可能会造成很多的误判，而改进小波时频图则可以给出非常准确的结果，具有很好的实用价值。
模型训练的速度也是实际应用时需要考虑的一个重要因素，图9比较了stft时频谱、小波+hilbert和改进小波时频谱作为特征输入时，在训练集上cnn模型同样进行1500次迭代时的准确率。从图中可以看出，采用改进小波时频图作为特征输入，cnn模型的训练速度非常快，300次迭代，模型的准确率已经达到80％左右，600次迭代之后准确率在90％以上，1000次迭代以内准确率已经接近100％。而300次迭代时，小波+hilbert和stft时频图的准确率均低于60％；训练结束时小波+hilbert变换的准确率最多能达到92％左右，而stft时频图的准确率只有88％左右。使用改进的小波时频图训练的cnn模型能够在滚动轴承产生微弱故障的早期对故障类别进行非常准确的判断，在不同的轴承也能给出正确的判断结果，而且训练的速度也很快，具有很好的实际应用价值。