CN113405799B - 基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法 - Google Patents

Abstract

本发明公开了一种基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法：首先，采用短时傅里叶变换STFT将原始振动数据转换为二维时频图，充分提取健康状态信息，对得到的时频图进行灰度处理以便后续模型输入；其次，利用高斯可视单元构造了新的高斯矩阵变量卷积深度信念网络GMVCDBN模型，对转换后的时频信息自动构造健康状态指标；最后，拟合健康状态指标概率分布函数，自学习预警控制极限。将本发明应用于滚动轴承实验振动信号的分析，能够实现滚动轴承的早期故障检测，降低误报率和漏报率，具有较好的泛化性，能更有效地跟踪故障传播趋势和估计轴承的动态退化行为。

Description

基于健康状态指标构建和故障预警限自学习的轴承早期故障 检测方法

技术领域

本发明涉及一种基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，属于机械装备的故障预测与健康管理领域。

背景技术

预测性维护是智能制造领域潜在爆发点。设备早期故障预测作为其关键一环，成为近年来工程领域研究的热点，发展出一系列极具应用前景的故障预测方法。然而机械系统采集到的振动信号通常是复杂的、非平稳的、背景噪声大的信号，健康状态演化趋势通常是非单调或者甚至是突变的，因此有效地从振动信号里学习能表征设备健康表征指标，以提前检测、告警设备即将发生的轻微或不正常故障征兆，指导设备预测性维修，是一个很大的挑战。

与许多可用的故障诊断方法不同，关于机械零部件的早期故障预测处于起步阶段，相关研究也在迅速发展。目前，早期故障预测模型可分为物理模型和数据驱动模型。基于模型的方法需要根据设备特性，搭建精确的物理模型对设备健康状况进行预测。其局限性在于对研究对象的先验知识要求非常高、泛化能力差。而且对于安装在结构复杂、运行环境恶劣以及部件之间相互影响较大的复杂机械设备建立完备的物理失效预测模型非常困难，使得该方法的使用和推广存在很大限制。

数据驱动的预测方法通过分析设备的传感器监测数据，得到设备潜在的退化规律，从而预测设备的健康状态演化趋势。它仅仅需要收集足够的性能退化数据，并将收集到的数据转换为相关健康状态特征或性能退化模型，以对未来数据进行预测。数据驱动的预测方法对设备相关先验知识的要求比较低，能够很好地解决在设备失效机理复杂以及难以建立精确模型时的预测问题，效率高且模型泛化能力好，获得越来越多的关注和研究。数据驱动的预测方法又可分为三类：第一类是统计模型预测方法，如灰色模型(Grey Model，GM)、粒子滤波(Particle Filter，PF)和贝叶斯预测模型(Bayesian prediction model，BPM)等；第二类是浅层机器学习预测方法，如人工神经网络(Artificial Neural Network,ANN)和支持向量机(Support Vector Machine,SVM)等；第三类是深度学习预测方法，如深度置信网络(Deep Belief Networks，DBN)和循环神经网络(Recurrent Neural Network,RNN)等。

虽然前两类数据驱动模型在健康状态预测中得到了广泛的应用，但它们仍然存在以下固有不足：SVM和ANN被认为具有浅层结构，很难有效地学习到健康状态演化过程中的高度非线性关系。通常情况下，由于运行条件的变化，轴承健康状态演化是一个非线性、非平稳的时间序列，更糟的是，在轴承故障的早期阶段，健康状态演化信息往往被噪声淹没。在这种情况下，采用浅层架构进行准确的故障预测是一个很大的挑战。因此，对轴承故障预测问题中复杂的非线性关系进行深入的研究是十分必要的。

深度学习模型是一种拥有多个非线性映射层级的深度神经网络模型，它们比浅层学习模型更能有效、灵活地表达任何高度变化的非线性函数，能够对输入信号逐层抽象并提取特征。换句话说，深度学习方法可以摆脱对手工信号处理技术和领域专家的依赖。深度置信网络(DBN)、循环神经网络(RNN)及其变体作为流行的深度学习预测方法，逐渐被应用于机械健康状态预测中。DBN 是一种生成式神经网络，具有很强的无监督特征学习能力。

即使以往故障预测取得一定成果，但大多数关注的是预测算法的改进，而基于深度学习的故障预测最重要的前提就是从振动信号中提取能够全面表征健康状态的特征指标(health state indicator，HI)，故现有方法存在以下不足： (1)现有的健康状态评估指标多是通过从原始振动信号中手动提取的部分时频域统计特征或熵特征并进行数据融合构造的，仅依靠某个域或某些域的特征量作为轴承健康状态评估指标，不能满足对其退化信息的全面描述，而且，每个原始特征只对某个退化阶段的某个故障敏感，也较难准确地捕捉到轴承退化过程的故障信息。(2)在振动特征信号中，有一些关键特征相互独立，对轴承的某种特定故障信息比较敏感，可以很好地描述轴承的退化状态，应该充分利用。然而，还有一些冗余或者不相关的特征，它们对轴承的退化状态不够敏感，会降低故障状态趋势预测的计算效率，甚至干扰评估和预测结果。所以在进行特征信号集输入到健康状态预测模型之前，需要利用特征评估方法将敏感特征从特征集中筛选出来，剔除特征集中不相关的冗余信息，而现有的特征选择和评估在很大程度上都取决于专家的工程经验。(3)健康状态指标一般是针对特定退化过程构建的，因此不能很好地推广到其他退化过程。

卷积神经网络(Convolutional neural network，CNN)因其具有表征学习能力和位移不变特征，以及通过卷积连接的权值共享策略使得训练出来的模型的泛化能力更强等优点。CNN一般是监督训练式的判别模型。而DBN则是一种生成式神经网络，具有很强的无监督特征学习能力，但要求输入数据具有平移不变性，使用小核的离散卷积来替换矩阵乘法是解决具有空间平移不变性或时间结构的输入问题的标准方式。虽然无监督学习模式的卷积深度置信网络 (Convolutional Deep belief network，CDBN)方法在机械状态监测和故障诊断方面取得了良好的效果，然而，将以往的研究CDBN应用于轴承健康状态指标构建时仍存在一个主要问题:隐藏层节点取0或1值，由此模型得到的隐藏层输出在一定程度上损失了原始数据的信息，不能全面表征健康状态。因此，CDBN模型的改进和健康状态指标的构建已成为一项紧迫的任务。

发明内容

本发明所要解决的技术问题是：提供一种健康状态指标自动构建与预警控制限自学习的滚动轴承早期故障预测方法，克服上述现有问题之不足，能有效解决现有的滚动轴承故障预测研究中存在的较难完备捕捉到轴承退化过程的故障特征、预测模型鲁棒性不强、漏报警和误报警率较高等问题。

为解决上述技术问题，本发明所采用的技术方案为：基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，包括以下步骤：

S1：时频域故障特征提取

S11：获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值；

S12：对得到的时频图进行灰度处理以便后续模型输入，将三通道图像处理为单通道灰度图像，调用MATLAB中的rgb2gray(RGB)函数，消除色调和饱和度信息，同时保留亮度，将RGB图像转换为灰度图；通过计算R、G和B 分量的加权和：0.2989*R+0.5870*G+0.1140*B，将RGB值转换为灰度值，图像像素处理为32*32；

S2：基于MVGCDBN的健康状态表征指标构建

采用矩阵变量高斯卷积深度置信网络MVGCDBN用于健康状态指标的构建：

模型采用由一个可视层V和一个隐藏层H构成的两层结构；可视层是一个 N_V×N_V二进制单元矩阵，且包含K个N_W×N_W(N_W＝N_V-N_H+1)的卷积核Wk；隐藏层是由K个N_H×N_H二进制单元矩阵的特征映射面组成，标记其中第k个特征映射面H^k的卷积核为W^k，偏置为b_k，所有的可视层共享偏置a；

假设可视层和隐藏层每个单元分别服从方差为σ₁ ²和γ₁ ²的高斯分布，则模型的能量函数如式(1)：

式中v和h分别代表所有可视层和隐藏层的状态向量；v_i,j是矩阵v第i行第j 列的元素，

是第k个隐藏单元第i行第j列的元素，

是第k个滤波器第r 行第s列的元素；

是

的翻转矩阵，

*表示卷积运算，·表示元素乘积后求和，即A·B＝trA^TB；

隐藏层和可视层单元条件激活概率分别为：

模型训练的目的是确定最优参数集，采用基于Gibbs采样的对比散度CD 算法更新模型的参数；根据计算得到的隐藏层条件概率对隐藏层节点进行采样，对隐藏层节点按均值为

方差为σ₁ ²的高斯分布进行采样，不再是将计算得到的概率值与0到1直接的随机数进行比较，得到0 或1值，此方法得到的采样数据不再是0或1，而是在均值附近采样的实数；可视层节点的采样与隐藏层采样方法相同；

输入：训练样本即可视层向量v，隐藏层中滤波器个数K，学习率、正则化常数ζ、稀疏度系数ρ、最大迭代次数Q；

输出：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k；

随机初始化：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k为随机的较小数值；

S3：自学习故障预警控制限

S31：分布拟合

(1)将“正常”阶段的HI数据做归一化处理；

(2)采用最大似然估计计算“正常”阶段的HI数据的正态分布及伽马分布的形状参数；

若随机变量x(HI)服从一个位置参数为μ、尺度参数为σ的概率分布，其概率密度函数为：

则称随机变量x服从参数为μ，σ²的正态分布，记作：x～N(μ,σ²)；

若随机变量x服从一个形状参数为α、逆尺度参数为β的概率分布，其概率密度函数为：

令x～Γ(α,β)，且令

即

其中α＞0，β＞0；则称随机变量X服从参数为α，β的伽马分布，记作：x～G(α,β)；

S32：拟合优度检验：对拟合的分布进行检验，确定HI所服从的分布；

S33：计算HI上侧分位数对应的阈值来确定归一化的控制限，再对其进行反归一化求得故障预警控制限，实现对故障点的早期预测。

作为一种优选的方案，所述步骤S11采用短时傅里叶变换STFT获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值；设g_τν(t)是一个以(τ,ν)为中心的局部化标准单窗函数，即低通滤波，在时域上集中于区间[τ-Δg,τ+Δg]，其中e^iνt是频移因子；

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R (a)

式中T_x(τ,ν)表示对窗口内的信号进行傅里叶变换，x(t)是需要被转换的原始信号，

看作是用一个以τ为中心，宽度为2Δg的窗口截取信号x(t)的一部分，当τ移动时，信号的不同部分将被依次取入窗口；STFT解释为：先将信号进行频移x(t)e^-iνt，然后再通过g_τν(t)进行低通滤波，从而得到在频率ν处的时域表示；其中，采用宽度为25ms的汉明窗函数，帧的位移大小为10ms。

作为另一种优选的方案，所述步骤S11对原始信号x(t)进行小波变换获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值，具体如下：

式中，a₀是尺度因子，a₀＞0；b₀是时移因子；

是小波基函数；小波变换的实质就是用小波基函数

通过内积运算*将信号x(t)分解为不同频带的子信号，得到其频谱。

作为一种优选的方案，所述步骤S32拟合优度检验采用K-S检验对拟合的分布进行检验，确定HI所服从的分布；此处令显著性水平α₁＝0.05，则置信度为1-α₁＝95％；

(1)令f(x)为频率分布，g(x)理论分布，假设H₀：数据符合理论分布；

(2)计算D＝max|f(x)-g(x)|，当实际观测值D＞D(n，α₁)则拒绝H₀，否则则接受H₀假设，即当

接受原假设。

作为另一种优选的方案，所述步骤S32拟合优度检验采用χ²检验对拟合的分布进行检验，确定HI所服从的分布；此处令显著性水平α₁＝0.05，则置信度为1-α₁＝95％，

式中，f₀表示第i个组的样本观测值的实际观测频数；f_e为理论期望的频数； nbins是组数；由公式(d)可以看出，如果观测结果与期望结果一致则χ²＝0，说明拟合度处于完美情况；当χ²大到一定程度，即累计贡献度超出设定的置信度范围，拒绝原假设即数据符合理论分布，χ²越大越不接受原假设。

本发明的有益效果是：

本发明提出了一种基于健康状态指标构建和故障预警控制限自学习的轴承早期故障预测方法，构建能够在边缘计算设备“端”运行的早期故障预测模型，仅需要正常运行状态的一段历史数据，不需要外部专家干预和人工分析即可准确获取早期故障信息，在信息物理系统环境下实现预测性维护提供技术支撑；以解决现有研究存在的所构建的健康状态评估指标不能满足对其退化信息的全面描述，较难准确地捕捉到轴承退化过程的故障信息的问题。

本发明相较于传统的固定阈值限报警方式能大幅提前故障预测时间，降低了误报警和漏报警率，并且具有较好的泛化性，这对于实现预测性维修、备件供应储备计划制定、避免非计划停机可能造成的事故损失具有重要的价值。

附图说明

图1是5个轴承全生命周期的水平方向振动加速度数据、峭度和均方根变化趋势图

图2是另外5个轴承全生命周期的水平方向振动加速度数据、峭度和均方根变化趋势图

图3是另外5个轴承全生命周期的水平方向振动加速度数据、峭度和均方根变化趋势图

图4是另外2个轴承全生命周期的水平方向振动加速度数据、峭度和均方根变化趋势图

图5是实施例的具体流程图

图6是基于MVGCDBN模型构建轴承的HI值

图7是轴承B1-3的故障预测结果(a)伽马分布(b)正态分布

图8是轴承B2-3的故障预测结果(a)伽马分布(b)正态分布

图9是轴承B3-3的故障预测结果(a)伽马分布(b)正态分布

图10是不同故障预测方法的对比结果

具体实施方式

下面结合附图1至10、表1至表4和具体实施例对本发明进行详细说明。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

采用PRONOSTIA试验台实测的滚动轴承加速寿命实验标准数据集验证所发明的早期故障预测方法。标准数据集收集三种工况下的数据，采样频率为 25.6kHz，每隔10秒采集一次数据，每次采样时长持续0.1s，即每次采样可以得到2560个加速度数据点。该数据集包含了17个滚动轴承在加速退化实验中所获取的运行至失效的振动监测信号，工况信息如表1所示，轴承样本描述如表2所示：

表1轴承工况表

表2轴承样本描述表

在滚动轴承加速寿命实验中，每个轴承性能自然退化的早期故障时间并没有明确被说明。然而，敏感的监测指标仍然可以在很大程度上检测出动态行为的明显变化。峭度和均方根(Kurtosis and root mean square，RMS)是两种最常用的量化轴承从健康到失效的统计参数，常被用作轴承早期故障预测指标。计算每个轴承的峭度和均方根，绘制其时域参数曲线，图1-4为17个轴承全生命周期的水平方向振动加速度数据、峭度和均方根变化趋势图。从图1-4可以看出：

B1-1的峭度在生命周期内波动很大，没有出现明显的增长趋势和变化；RMS 则在全生命周期内呈递增趋势，但没有突增点。因此峭度和RMS均无法判断轴承B1-1的早期故障点。B1-2的峭度在生命周期内波动很大，没有出现明显的增长趋势和变化；RMS则在827点以前处于较平稳阶段，827点后显示出明显增长趋势和变化。因此峭度无法判断轴承B1-2的早期故障，RMS认为轴承 B1-2在827点发生早期故障。B1-3的峭度在生命周期内波动很大，没有出现明显的增长趋势和变化；RMS则在全生命周期内呈递增趋势，但没有突增点，无法判断B1-3早期故障。B1-4的峭度在1083点以前处于较平稳阶段，1083 点后显示出明显增长趋势和变化；RMS在1090点以前处于较平稳阶段，1090 点后显示出明显增长趋势和变化。因此峭度认为轴承B1-4在1083点发生早期故障，RMS认为轴承B1-4在1090点发生早期故障。B1-5的峭度在2412点以前处于较平稳阶段，2412点后显示出明显增长趋势和变化；RMS在2450点以前处于较平稳阶段，2450点后显示出明显增长趋势和变化。因此峭度和RMS 分别预测轴承B1-4在2412点和2450点发生早期故障。B1-6的峭度和RMS 在1631点以前处于较平稳阶段，1631点后显示出明显增长趋势和变化。因此峭度和RMS均认为轴承B1-6在1631点发生早期故障。B1-7的峭度和RMS 在2210点以前处于较平稳阶段，2210点后显示出明显增长趋势和变化。因此峭度和RMS均认为轴承B1-7在2210点发生早期故障。

B2-1的峭度在生命周期内波动很大，没有出现明显的增长趋势和变化；RMS 则在前中后期均有趋势性波动，在875点后显示出明显增长趋势和变化。因此峭度无法判断轴承B2-1的早期故障，RMS认为轴承B2-1在875点发生早期故障。B2-2的峭度除了早期运行的不平稳数据外，752点以前处于较平稳阶段， 752点后显示出明显增长趋势和变化；RMS在752点后显示出明显增长趋势和变化。因此峭度和RMS都认为轴承B2-2在752点发生早期故障。B2-3的峭度在全生命周期内波动很大，没有出现明显的增长趋势和变化；RMS除早期运行的不平稳数据外，在1947点以前处于较平稳阶段，1947点后显示出明显增长趋势和变化，因此，峭度无法判断轴承B2-3的早期故障点，而RMS认为轴承B2-3在1947点发生早期故障。B2-4的峭度在262点以前处于较平稳阶段， 262至410点之间波动很大，410点后又处于较平稳阶段，在603点至失效阶段出现多次波动，但无明显增长趋势和变化；RMS除早期运行的不平稳数据外，在744点以前处于较平稳阶段，744点后显示出明显增长趋势和变化。因此峭度无法判断轴承B2-4的早期故障点，RMS认为在744点发生早期故障。 B2-5的峭度除早期运行的不平稳数据外，在2269点以前处于较平稳阶段，2269 点后显示出明显增长趋势和变化；RMS在530点之前波动很大，之后较平稳，但无显示出明显增长趋势和变化。因此峭度认为轴承B2-5在2269点发生早期故障，RMS无法判断轴承B12的早期故障点。B2-6的峭度除早期运行的不平稳数据外，在688点以前处于较平稳阶段，688点后显示出明显增长趋势和变化；RMS在688点以前处于较平稳阶段，688点后显示出明显增长趋势和变化。因此峭度和RMS均认为轴承B2-6在688点发生早期故障。B2-7的峭度在生命周期内有多处突增点，但无明显增长趋势和变化；RMS则在226点以前处于较平稳阶段，226点后显示出明显增长趋势和变化。因此峭度无法判断轴承B2-7的早期故障点，RMS认为轴承B2-7在226点发生早期故障。

B3-1的峭度和RMS在491点以前处于较平稳阶段，491点后显示出明显增长趋势和变化。因此峭度和RMS均认为轴承B3-1在491点发生早期故障。 B3-2的峭度出现多处波动，在1433点有明显增长趋势；RMS则在1614点以前处于较平稳阶段，1614点后显示出明显增长趋势和变化。因此，峭度认为轴承B3-2在1433点发生早期故障，RMS则认为轴承B3-2在1614点发生早期故障。B3-3的峭度在生命周期内有多处突增点，但没有出现明显的增长趋势和变化；RMS则在310点之前处于较平稳阶段，310点后显示出明显增长趋势和变化。因此，峭度无法判断轴承B3-3的早期故障点，RMS则认为轴承B3-3 在310点发生早期故障。

由于退化过程的随机性，多个轴承的全寿命振动信号存在一定差异。为了避免结果的偶然性和特殊性，选择三种不同工况下前2个轴承(即轴承编号为 1-1、1-2、2-1、2-2、3-1、3-2)作为训练集，其余轴承作为测试集。

本发明提出的一种早期故障预测方法实施例的具体流程图如附图5所示，主要包括以下步骤：

S1：时频域故障特征提取

S11：采用STFT获取17个滚动轴承的非平稳原始振动信号中的时频域信息(即，信号频率随时间序列变化的情况，各个时刻的瞬时频率及其幅值)。设g_τν(t)是一个以(τ,ν)为中心的局部化标准单窗函数，即低通滤波，在时域上集中于区间[τ-Δg,τ+Δg]，其中e^iνt是频移因子。

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R (1)

式中T_x(τ,ν)表示对窗口内的信号进行傅里叶变换，x(t)是需要被转换的原始信号，

可以看作是用一个以τ为中心，宽度为2Δg的窗口截取信号x(t)的一部分，当τ移动时，信号的不同部分将被依次取入窗口。STFT可以解释为：先将信号进行频移x(t)e^-iνt，然后再通过g_τν(t)进行低通滤波，从而得到在频率ν处的时域表示。其中，采用宽度为25ms的汉明窗函数，帧的位移大小为10ms。

S12：对得到的时频图进行灰度处理以便后续模型输入，将三通道图像处理为单通道灰度图像，调用MATLAB中的rgb2gray(RGB)函数，消除色调和饱和度信息，同时保留亮度，将RGB图像转换为灰度图。通过计算R、G和B 分量的加权和：0.2989*R+0.5870*G+0.1140*B，将RGB值转换为灰度值，图像像素处理为32*32。将提取的轴承时频域故障特征输入GMVCDBN模型里构建其HI。

S2：基于MVGCDBN的健康状态表征指标构建

标准的CDBN是一种新的分层生成模型，由多个卷积受限玻尔兹曼机(convolutional restricted Boltzmann machines，CRBMs)构造而成。CRBM是标准受限玻尔兹曼机(restricted Boltzmann machines，RBM)的扩展，采用局部连接和共享权值的方式，考虑输入数据的二维(2-D)结构和周期特性，适合进行数据的空间扩展，克服了所有可视单元必须通过不同的权值参数与所有隐藏单元相关联的缺点，可以直接从原始振动数据中通过卷积层和池化层交替使用来获取有效表征，自动提取隐藏在数据中的局部特征，组合抽象成高层特征，并通过全连接层构建HI。

标准的CDBN是由一个可视层V和一个隐藏层H构成的两层结构。可视层是一个N_V×N_V二进制单元矩阵，且包含K个N_W×N_W(N_W＝N_V-N_H+1)的卷积核 W_k(即滤波器)。隐藏层是由K个N_H×N_H二进制单元矩阵的特征映射面组成，标记其中第k个特征映射面H^k的卷积核为W^k，偏置为b_k，所有的可视层共享偏置a。

如果可视单位为二进制，则能量函数可定义为式(3)：

如果可视单位为实值，则能量函数可定义为式(4)：

式中v和h分别代表所有可视层和隐藏层的状态向量。v_i,j是矩阵v第i行第j 列的元素，

是第k个隐藏单元第i行第j列的元素，

是第k个滤波器第r 行第s列的元素。

是

的翻转矩阵，

*表示卷积运算，·表示元素乘积后求和，即A·B＝trA^TB。

在特征映射曲面H^k上，隐藏层单元在第m行第n列的条件激活概率可表示为式(5)：

同理，可视层单元的条件激活概率可表示为式(6)：

以上采用高斯可视单元和二进制隐藏层构建新的一维CDBN模型，由此模型得到的隐藏层输出在一定程度上损失了原始数据的信息。由Welling等人于 2004年发表在会议期刊《Proceedings of the Neural Information Processing Systems》上的论文中指出CRBM中的隐藏层节点和可视层节点的分布可以为任意的指数族分布。本发明令可视层和隐藏层均满足高斯分布，此模型假设可视层输入为高斯分布，能够更加准确的拟合输入数据的同时可以提取到原始样本的更加丰富的信息。另外，输入为图像、视频等高维数据时，需要先对输入数据进行向量化，这样做会破坏样本的空间内部结构。所以借鉴矩阵变量受限玻尔兹曼机(Matrix Variate Restricted Boltzmann Machine，MVRBM)的思想，提出了矩阵变量高斯卷积深度置信网络(Gaussian Matrix Variate Convolutional Deep BeliefNetwork，MVGCDBN)用于健康状态指标的构建。假设可视层和隐藏层每个单元分别服从方差为σ₁ ²和γ₁ ²的高斯分布，则改进模型的能量函数如式(7)：

隐藏层和可视层单元条件激活概率分别为：

模型训练的目的是确定最优参数集，与标准的CDBN相似，采用基于Gibbs 采样的对比散度(contrastive divergence，CD)算法更新模型的参数。根据计算得到的隐藏层条件概率对隐藏层节点进行采样，对隐藏层节点按均值为

方差为σ₁ ²的高斯分布进行采样，不再是将计算得到的概率值与0到1直接的随机数进行比较，得到0或1值，此方法得到的采样数据不再是0或1，而是在均值附近采样的实数。可视层节点的采样与隐藏层采样方法相同。

输入：训练样本(可视层向量)v，隐藏层中滤波器个数K，学习率、正则化常数ζ、稀疏度系数ρ、最大迭代次数Q。

输出：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k。

随机初始化：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k为随机的较小数值。

通过实验确定了该模型的体系结构，MVGCDBN模型的主要参数如表3 所示：

表3用于滚动轴承HI构建的参数

注：第一层和第二层CRBMs的符号不作区分

MVGCDBN构造的测试轴承数据集健康状态曲线的结果如图6所示。

S3：自学习故障预警控制限

S31：分布拟合

(1)将“正常”阶段的HI数据做归一化处理。

(2)采用最大似然估计计算“正常”阶段的HI数据的正态分布及伽马分布的形状参数。

若随机变量x(HI)服从一个位置参数为μ、尺度参数为σ的概率分布，其概率密度函数为：

则称随机变量x服从参数为μ，σ²的正态分布，通常记作：x～N(μ,σ²)。

若随机变量x(HI数据)服从一个形状参数为α、逆尺度参数为β的概率分布，其概率密度函数为：

令x～Γ(α,β)，且令

即

其中α＞0，β＞0。则称随机变量 X服从参数为α，β的伽马分布，通常记作：x～G(α,β)。

S32：拟合优度检验

采用K-S检验对拟合的分布进行检验，确定HI所服从的分布。此处令显著性水平α₁＝0.05。

(1)令f(x)为频率分布，g(x)理论分布，假设H₀：数据符合理论分布。

(2)计算D＝max|f(x)-g(x)|，当实际观测值D＞D(n，α₁)则拒绝H₀，否则则接受H₀假设，即当h＝0，

接受原假设。

S33：计算HI上侧分位数对应的阈值来确定归一化的控制限，再对其进行反归一化求得故障预警控制限，实现对故障点的早期预测。

因篇幅有限，此处仅展示每种工况下的第一个测试数据集(即轴承B1-3、 B2-3、B3-3)预警控制限的计算结果。为了验证所发明方法的优越性和泛化性，将发明的方法与常用的量化轴承从健康到故障的统计参数(峭度和RMS)进行比较，其它结果在表4做统计及对比分析。

在自学习控制限构建过程中，外部影响在采集过程中产生的尖峰误差一般情况下为5％，所以双侧分位数取0.05。如图7所示，是自学习的轴承B1-3故障预警控制限结果。采用K-S检验分别对伽马分布和正态分布进行拟合优度检验。首先检验伽马分布，由于h＝0，p＝0.5486，所以在显著性水平0.05下接受原假设，即认为B1-3的HI服从参数为α＝5.6695，β＝0.034572的伽马分布。再检验正态分布，由于h＝0，p＝0.6321，所以在显著性水平0.05下接受原假设，即认为B1-3的HI服从参数为μ＝0.19601，σ²＝0.082291的正态分布。综上，认为B1-3 的HI服从参数为μ＝0.19601，σ²＝0.082291的正态分布。如图7(b)所示，正常数据与故障数据已经发生了分离，前2110组数据除了刚开始运行的不平稳数据以外，均处于自学习控制限(HI＝0.33)以下，判断2210点为早期故障点，即轴承运转至约368分钟(轴承B1-3全寿命约396分钟)。基于算例B1-3数据，由图1-4所知，通过计算峭度、均方根无法确定其故障预警点。

如图8所示，是自学习的轴承B2-3故障预警控制限结果。经K-S检验认为 B2-3的HI服从参数为μ＝0.63004，σ²＝0.012808的正态分布。如图8(b)所示，正常数据与故障数据已经发生了分离，前411组数据除了刚开始运行的不平稳数据以外，均处于自学习控制限(值为0.84)以下，判断411点为早期故障点，即轴承运转至69分钟(轴承B2-3全寿命约326分钟)。基于算例B2-3数据，由图1-4所知，RMS确定判断1947点为故障预警点，峭度无法确定其故障预警点。

如图9所示，是自学习的轴承B3-3故障预警控制限结果。经K-S检验认为 B3-3的HI服从参数为的α＝3.0943，β＝0.10952伽马分布。如图9(a)所示，正常数据与故障数据已经发生了分离，前310组数据除了刚开始运行的不平稳数据以外，均处于自学习控制限(值为0.70)以下，判断310点为早期故障点，即轴承运转至约52分钟(轴承B3-3全寿命约72分钟)。基于算例B3-3数据，由图1-4所知，RMS确定判断310点为故障预警点，峭度无法确定其故障预警点。

表4不同故障预测方法的结果对比

注：“/”表示早期故障不能被预测出来

从表4和图10可以看出，本发明所提方法在早期故障预测上的整体表现优于峭度和RMS，这具体表现在：

(1)对于轴承B1-1、B1-3峭度和RMS均无法预测出它们的早期故障；而所发明的方法判断出B1-1和B1-3的早期故障点分别是：2433点(即轴承运行到约406分钟)、2210点(即轴承运行到约368分钟)。

(2)对于轴承B1-2、B2-1、B2-3、B2-4、B2-7和B3-3，峭度无法判断出它们的早期故障点，RMS判断出来它们的早期故障点分别为：827点(轴承运行到约138分钟)、875点(轴承运行到约146分钟)、1947点(轴承运行到约 325分钟)、744点(轴承运行到约124分钟)、226点(轴承运行到约38分钟)、 310点(轴承运行到约52分钟)；而所发明的方法判断出它们的早期故障点分别为：720点(轴承运行到约120分钟)、790点(轴承运行到约132分钟)、 411点(轴承运行到约69分钟)、720点(轴承运行到约120分钟)、180点(轴承运行到约30分钟)、310点(轴承运行到约52分钟)。与RMS方法相比，所发明的方法与B3-3轴承的RMS检测结果一致。对于剩余轴承，所发明的方法比RMS分别可提前18分钟、14分钟、4分钟、8分钟预测出轴承B1-2、 B2-1、B2-4和B2-7的早期故障。

(3)对于轴承B2-5，RMS无法判断出它的早期故障点，所发明的方法与峭度均可判断出它的早期故障点是2269点(轴承运行到约378分钟)。

(4)对于轴承B2-2和B2-6，峭度、RMS和本发明的方法可同时判断出来它们的早期故障点分别是752点(轴承运行到约125分钟)和688点(轴承运行到约115分钟)。

(5)对于轴承B1-4、B1-5、B1-6、B1-7、B3-1和B3-2，峭度判断出来它们的早期故障点分别为：1083点(轴承运行到约181分钟)、2412点(轴承运行到约402分钟)、1631点(轴承运行到约272分钟)、2210点(轴承运行到约368分钟)、491点(轴承运行到约82分钟)、1433点(轴承运行到约239分钟)；RMS判断出来它们的早期故障点分别为：1090点(轴承运行到约182分钟)、2450点(轴承运行到约408分钟)、1631点(轴承运行到约272分钟)、 2210点(轴承运行到约368分钟)、491点(轴承运行到约82分钟)、1614点(轴承运行到约269分钟)；所发明的方法判断出它们的早期故障点分别为： 1030点(轴承运行到约172分钟)、2410点(轴承运行到约402分钟)、1600 点(轴承运行到约267分钟)、2136点(轴承运行到约356分钟)、480点(轴承运行到约80分钟)、1330点(轴承运行到约223分钟)，均早于前两种方法。与峭度相比，对于轴承(B1-4、B1-5、B1-6、B1-7、B3-1和B3-2)，所发明的方法可提前9分钟、0分钟、5分钟、12分钟、2分钟和16分钟预测出其早期故障。与RMS相比，所发明的方法可提前10分钟、6分钟、5分钟、12分钟、 2分钟和16分钟预测其早期故障。

上述的实施例仅例示性说明本发明创造的原理及其功效，以及部分运用的实施例，而非用于限制本发明；应当指出，对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (5)

1.基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，包括以下步骤：

S1：时频域故障特征提取

S11：获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值；

S12：对得到的时频图进行灰度处理以便后续模型输入，将三通道图像处理为单通道灰度图像，调用MATLAB中的rgb2gray(RGB)函数，消除色调和饱和度信息，同时保留亮度，将RGB图像转换为灰度图；通过计算R、G和B分量的加权和：0.2989*R+0.5870*G+0.1140*B，将RGB值转换为灰度值，图像像素处理为32*32；

S2：基于MVGCDBN的健康状态表征指标构建

采用矩阵变量高斯卷积深度置信网络MVGCDBN用于健康状态指标的构建：

模型采用由一个可视层V和一个隐藏层H构成的两层结构；可视层是一个N_V×N_V二进制单元矩阵，且包含K个N_W×N_W(N_W＝N_V-N_H+1)的卷积核W_k；隐藏层是由K个N_H×N_H二进制单元矩阵的特征映射面组成，标记其中第k个特征映射面H^k的卷积核为W^k，偏置为b_k，所有的可视层共享偏置a；

假设可视层和隐藏层每个单元分别服从方差为σ₁ ²和γ₁ ²的高斯分布，则模型的能量函数如式(1)：

式中v和h分别代表所有可视层和隐藏层的状态向量；v_i,j是矩阵v第i行第j列的元素，

是第k个隐藏单元第i行第j列的元素，

是第k个滤波器第r行第s列的元素；

是

的翻转矩阵，

*表示卷积运算，·表示元素乘积后求和，即A·B＝trA^TB；

隐藏层和可视层单元条件激活概率分别为：

模型训练的目的是确定最优参数集，采用基于Gibbs采样的对比散度CD算法更新模型的参数；根据计算得到的隐藏层条件概率对隐藏层节点进行采样，对隐藏层节点按均值为

方差为σ₁ ²的高斯分布进行采样，不再是将计算得到的概率值与0到1直接的随机数进行比较，得到0或1值，此方法得到的采样数据不再是0或1，而是在均值附近采样的实数；可视层节点的采样与隐藏层采样方法相同；

输入：训练样本即可视层向量v，隐藏层中滤波器个数K，学习率、正则化常数ζ、稀疏度系数ρ、最大迭代次数Q；

输出：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k；

随机初始化：连接权重参数W^k、可视层偏置a、隐藏层偏置b_k为随机的较小数值；

for q＝1 to Q do

for k＝1 to K do

Step 1:更新权重参数

Step 2:利用正则化项的梯度更新权重

Step 3:重复Step 1和Step 2，直到收敛或q＝Q

end for

end for

S3：自学习故障预警控制限

S31：分布拟合

(1)将“正常”阶段的HI数据做归一化处理；

(2)采用最大似然估计计算“正常”阶段的HI数据的正态分布及伽马分布的形状参数；

若随机变量x(HI)服从一个位置参数为μ、尺度参数为σ的概率分布，其概率密度函数为：

则称随机变量x服从参数为μ，σ²的正态分布，记作：x～N(μ,σ²)；

若随机变量x服从一个形状参数为α、逆尺度参数为β的概率分布，其概率密度函数为：

令x～Γ(α,β)，且令

即

其中α＞0，β＞0；则称随机变量X服从参数为α，β的伽马分布，记作：x～G(α,β)；

S32：拟合优度检验：对拟合的分布进行检验，确定HI所服从的分布；

S33：计算HI上侧分位数对应的阈值来确定归一化的控制限，再对其进行反归一化求得故障预警控制限，实现对故障点的早期预测。

2.如权利要求1所述的基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，其特征在于：所述步骤S11采用短时傅里叶变换STFT获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值；设g_τν(t)是一个以(τ,ν)为中心的局部化标准单窗函数，即低通滤波，在时域上集中于区间[τ-Δg,τ+Δg]，其中e^iνt是频移因子；

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R(a)

式中T_x(τ,ν)表示对窗口内的信号进行傅里叶变换，x(t)是需要被转换的原始信号，

看作是用一个以τ为中心，宽度为2Δg的窗口截取信号x(t)的一部分，当τ移动时，信号的不同部分将被依次取入窗口；STFT解释为：先将信号进行频移x(t)e^-iνt，然后再通过g_τν(t)进行低通滤波，从而得到在频率ν处的时域表示；其中，采用宽度为25ms的汉明窗函数，帧的位移大小为10ms。

3.如权利要求1所述的基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，其特征在于：所述步骤S11对原始信号x(t)进行小波变换获取滚动轴承的非平稳原始振动信号中的时频域信息，即信号频率随时间序列变化的情况和各个时刻的瞬时频率及其幅值，具体如下：

式中，a₀是尺度因子，a₀＞0；b₀是时移因子；

是小波基函数；小波变换的实质就是用小波基函数

通过内积运算*将信号x(t)分解为不同频带的子信号，得到其频谱。

4.如权利要求1-3中任一项所述的基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，其特征在于：所述步骤S32拟合优度检验采用K-S检验对拟合的分布进行检验，确定HI所服从的分布；此处令显著性水平α₁＝0.05，则置信度为1-α₁＝95％；

(1)令f(x)为频率分布，g(x)理论分布，假设H₀：数据符合理论分布；

(2)计算D＝max|f(x)-g(x)|，当实际观测值D>D(n，α₁)则拒绝H₀，否则则接受H₀假设，即当h＝0，

接受原假设。

5.如权利要求1-3中任一项所述的基于健康状态指标构建和故障预警限自学习的轴承早期故障检测方法，其特征在于：所述步骤S32拟合优度检验采用χ²检验对拟合的分布进行检验，确定HI所服从的分布；此处令显著性水平α₁＝0.05，则置信度为1-α₁＝95％，

式中，f₀表示第i个组的样本观测值的实际观测频数；f_e为理论期望的频数；nbins是组数。

Images