CN107831013B - 利用概率主分量分析增强循环双谱的轴承故障诊断方法 - Google Patents

Abstract

本发明属于机械设备故障诊断领域，本发明公开了一种概率主分量增强循环双谱的轴承故障诊断方法，首先建立能反映原始信号与主元向量潜在联系的概率主元模型；其次利用概率主元模型对原始信号进行消噪，完整地保留信号的有用成分，并大幅度地提高信噪比，解决故障信号易受噪声干扰问题；最后对消噪信号作循环双谱分析，并从单一循环频率双谱的等高线图中诊断出轴承故障。本发明方法一方面利用概率主分量分析消噪能有效地提高信噪比，抑制噪声，使轴承共振频率附近的故障频率调制现象凸显出来，从而增加循环双谱形成的六边形结构清晰度；另一方面，从单一循环频率双谱的等高线图中可直观检测出轴承故障类型。

Description

利用概率主分量分析增强循环双谱的轴承故障诊断方法

技术领域

本发明属于机械设备故障诊断技术领域，尤其是一种利用概率主分量分析增强循环双谱的轴承故障诊断方法。

背景技术

旋转机械一般由支承部件(轴承)、旋转部件(轴系)和传动部件(齿轮)组成，是一种现代设备不可或缺的重要的动力装置，其结构形式复杂、多样，被广泛应用于现代机械加工领域、核电、风能、水电能源领域、大型机组设备和航空装备上。据统计，旋转机械发生的故障中，轴承为易损部件，其故障占旋转机械故障的19％。因此，对轴承的工作状态进行实时监测，以便早期发现异常故障是十分有价值和现实意义的。

然而，旋转机械轴承的早期故障特征十分微弱，极易被各种强背景噪声所掩盖。另外，轴承的故障信号通常是一种具有冲击特征的循环平稳信号，具有非高斯性和较宽的频带。要想从强背景噪声下准确地提取出故障特征，采用传统的频谱分析法几乎不可能，必须采用能处理非线性、非高斯和非白色的加性噪声的现代信号处理方法。而现代信号处理方法在机械故障诊断领域迫切需要解决的问题之一就是可靠的轴承故障特征提取方法。

概率主分量分析(Probabilistic Principal Component Analysis,PPCA)，本质上是一种因子分析隐含变量模型，它将信号的所有变量都用概率分布表达，通过构建概率主元模型，并用最大期望(Expectation Maximization,EM)算法对模型相关参数估计，最后由概率主元模型重新生成观测数据，从而实现消噪并提高信噪比。循环双谱分析需要对循环平稳信号进行非线性变换，并产生有限强度的正弦波。然而，信号本身并不含有任何有限强度的加性正弦波分量，为了产生正弦波分量，需在非线性变换的基础上进行正弦波抽取运算，产生正弦波的频率就是循环频率。由于循环双谱分析易受非高斯噪声影响，直接影响到故障诊断的成功和准确与否。为此，本文为了克服循环双谱易受非高斯噪声影响的不足，对信号采用概率主分量分析增强，提出了一种利用概率主分量分析增强循环双谱的轴承故障诊断方法，有关这方面研究，目前尚无报道。

发明内容

为了克服现有技术的不足，本发明提供了一种利用概率主分量分析增强循环双谱的轴承故障诊断方法，该诊断方法可以准确检测轴承故障情况。

为了实现上述目的，本发明采用的技术方案是：一种利用概率主分量分析增强循环双谱的轴承故障诊断方法，其特征在于，包括以下步骤：①建立能反映原始信号与主元向量潜在联系的概率主元模型；②利用概率主元模型对原始信号进行消噪，获得高信噪比信号，完整地保留信号的有用成分，并大幅度地提高信噪比，解决故障信号易受噪声干扰问题；③对消噪信号作循环双谱分析，通过计算三阶累积量、正弦抽取运算、二维FFT变换，从而获得单一循环频率双谱的等高线图，从单一循环频率双谱的等高线图中规则的六边形顶点坐标与故障特征频率比较，最终诊断出轴承故障；

其中，步骤①包括以下内容：

概率主元模型是在因子分析法隐含变量模型的基础上，假设高斯噪声的方差各

向同性，由因子向量(主元向量)叠加噪声产生的生成模型，其数学形式为：

X＝P·u+E (1)

式(1)中，X＝{s₁,s₂,…,s_m}∈R^n×m是由原始信号s生成得n×m阶矩阵，其中：n 为嵌入维数(原始变量个数)，m为采样数；P＝{p₁,p₂,…,p_k}∈R^n×k为待求的n×k 阶负荷矩阵(主元矩阵)，k为主元个数，且k<n；u＝{u₁,u₂,…,u_m}∈R^k×m为主元变量，u_i(i＝1,2,…m)为主元变量的分量；E为服从N(0,σ²I)分布的高斯噪声，σ²为噪声方差；因此，X为服从N(0,PP^T+σ²I)的高斯分布；

确定概率模型后，分别计算主元变量、原始变量在给定主元变量、主元变量在给定原始变量下的分布，则可获得概率主元模型；具体计算过程如下：

假设主元变量u服从多元标准正态分布，其概率分布为：

原始变量在给定主元变量下的条件分布为：

原始变量的分布为：

式(4)中，C＝PP^T+σ²I是n×n阶原始变量的方差，根据贝叶斯定理，得到主元变量在给定原始变量下的条件分布:

式(5)中，M＝P^TP+σ²I，是已降维的k×k阶方阵；由式(2)～(5)可知，只要确定参数P和σ²，就可求得概率模型；由隐含模型经典参数估计EM算法估算参数：

式(6)和(7)中，为原始变量的协方差矩阵，tr(·)是矩阵的迹，即对主对角线上的所有元素取和，反复迭代式(6)和(7)直至收敛，即可求得负荷向量P和高斯噪声的方差σ²；

步骤②包括以下内容：

根据建立的概率主元模型得到的负荷向量P，将原始信号X分解为P的外积之和，得到原始信号的主元模型，其数学表示为：

式(8)中，每一个主元实际上是原始信号X在和这个主元分量对应的负荷向量方向上的投影；主元向量的对应的特征值反映了原始信号X在p_i方向的覆盖程度；将每一个主元向量按其特征值大到小排序，则其对应的负荷向量就表示出原始信号X变化方向的大小，选择包含主要信息量的前(1,2,…,k)个负荷向量重构数据，最终得到降维的原始信号矩阵，实现信号消噪的目的；

概率主元模型对原始信号进行消噪取决于嵌入维数(原始变量个数)n和主元个数k；k可以由方差累积贡献率(cumulative variance contribution rate,CVCR)计算；单一主元贡献率计算公式如下：

式(9)中，λ为原始变量的协方差矩阵S的特征值；因此，全部主元贡献率计算公式为：

一般而言，k为的阈值达到85％时的最小值(文献：I.T.Jolliffe,Principal Component Analysis(2nd edition).Springer-Verlag,New York,2008.)；

当k确定了之后，由k<n确定最小的嵌入维数(原始变量个数)；

步骤③包括以下内容：

通过计算三阶累积量、正弦抽取运算、二维FFT变换，从而获得单一循环频率双谱的等高线图，具体为：

在频域中计算三阶累积量公式：

式(11)中，u(f)表示u(n)的傅里叶变换，E(·)表示求期望，u^*表示求共轭，α为引入的循环频率；由式(11)计算循环双谱可分为两个步骤：双谱估计；正弦抽取运算；

经典的双谱估计有非参数化法和参数化法，非参数化法又有直接法和间接法两种；本发明首先中采用间接法进行双谱估计，具体算法如下：

令u(1),u(2),...,u(m)是一组主元变量分量的采样信号，并设f_s是采样频率，N₀是总的频率采样点数，则Δf＝f_s/N₀是双谱区域沿水平和垂直方向上的频率采样间隔；

第1步：将信号u(1),u(2),...,u(m)分成长度为m的k段，即N＝k×m,并减去各段的均值；

第2步：设{u⁽ⁱ⁾(j),j＝1,2,...,m}是第i段数据，估计各段的三阶累积量：

式(12)中，m₁＝max(0,-k,-l),m₂＝min(m-1,m-1-k,m-1-l)；

第3步：计算所有段的三阶累积量的均值作为整个观测数据组的三阶累积量估计值

第4步：进行双谱估计：

式(14)中，L＜M-1,ω(k,l)是二维窗函数；

其次，进行正弦抽取运算；

从式(11)定义式可知，在频域定义的循环双谱相当于将双谱的频率向左平移循环频率α后得到；根据傅里叶移频定理，在时域，相当于将信号的三阶循环累积量乘以因子e^jαt后得到的，这就是正弦抽取运算，得到的循环双谱的统计量定义为：

式(15)中,T为周期，在时域中计算三阶累积量公式：

式(16)中，一、二、三阶矩分别为：

进行正弦抽取运算的关键是要预先估计循环频率α，它直接决定了双谱分析的成功与故障频率诊断的准确性；实际应用循环双谱时，采用线性搜索的方法确定循环频率：

第1步：在T＝[0,2π)内选择特定的循环频率集{α}，以定步长Δα搜索，按式 (16)计算每一循环频率α_j∈{α}的三阶循环累积量

第2步：对循环频率集{α}内计算的所有三阶累积量求和，得到

第3步：对进行二维傅里叶变换，得到属于循环频率区间(α₁,α₂)的循环双谱：

步骤④包括以下内容：根据获得单一循环频率双谱的等高线图，从单一循环频

率双谱的等高线图中规则的六边形顶点坐标与故障特征频率比较，最终诊断出轴承故障。

采用上述方案，本发明的优点如下：本方法一方面利用概率主分量分析消噪能有效地提高信噪比，抑制噪声，使轴承共振频率附近的故障频率调制现象凸显出来，从而增加循环双谱形成的六边形结构清晰度；另一方面，从单一循环频率双谱的等高线图中可靠地检测出轴承故障类型。

下面结合附图对本发明作进一步描述。

附图说明

附图1为概率主分量分析增强循环双谱的轴承故障诊断流程图；

附图2为原始仿真信号波形和频谱图，(a)波形图，(b)频谱图；

附图3为概率主分量分析消噪后的信号波形和频谱图，(a)波形图，(b)频谱图；

附图4为消噪后信号在α＝1000Hz循环双谱等高线和三维图，(a)等高线图，(b)三维图；

附图5为原始信号在α＝1000Hz循环双谱等高线和三维图，(a)等高线图，(b)三维图；

附图6为轴承内圈故障原始信号波形图，(a)时域图，(b)频谱图；

附图7为概率主分量分析消噪后信号波形和频谱图，(a)时域图，(b)频谱图；

附图8为循环双谱等高线和三维图，(a)等高线图，(b)等高线局部放大图，(c)三维图；

附图9为未消噪信号循环双谱等高线和三维图，(a)等高线图，(b)三维图。

具体实施方式

本发明的保护范围不局限于下述具体实施方式，本领域一般技术人员根据本发明公开的内容，可以采用其他多种具体实施方式实施本发明的，或者凡是采用本发明的设计结构和思路，做简单变化或更改的，都落入本发明的保护范围。

本发明的具体实施例如图1＝9所示。

具体实施例1：为验证概率主分量分析增强循环双谱的轴承故障诊断方法的正确性，根据轴承故障特征，建立数字仿真模型：

式(20)中，n(t)为零均值的平稳随机信号，h(t)为由轴承故障产生的一系列重复脉冲冲击串，A_i为随机调制幅值，T为相邻两次脉冲的冲击间隔，τ_i表示第i冲击相对于故障平均周期T的微小波动，i为冲击数。

轴承对每个脉冲的响应，定义为：

h(t)＝e^-Bt sin 2πf_nt (21)

式(21)中，B是由轴承共振频率、质量系数和应力松弛时间决定的系数，f_n是轴承的中心频率。

轴承故障的严重程度受到各个滚动体上的径向载荷分布、装配的动刚度、滚动体和套圈的圆度误差、滚动体组差等随机变量的影响，所以冲击脉冲的幅值A_i应为随机调制幅值，表示为：

A_i＝cos(2πf_At+φ_A)+C_A+randn(t) (22)

式(22)中，φ_A和C_A为任意常数，f_A为调制频率(外圈故障时f_A＝0；内圈故障时f_A等于轴的转频；滚动体故障时f_A等于保持架转频)，高斯随机噪声 randn(t)模拟幅值的随机性。

零均值的平稳随机信号n(t)用正弦型随机相位信号模拟，定义为：

n(t,φ_i)＝Acos(ω₀t+φ_i) φ_i∈(-π,π) (23)

式(23)中，A、ω₀为常数，φ_i为(-π,π)上均匀分布的随机变量。

设定轴承故障冲击周期T＝0.02，即故障频率为50Hz，调制频率 f_A＝1480/60Hz，采样频率f_s＝4000Hz，相对故障平均周期的微小波动τ_i＝0.0155T，轴承的中心频率f_n＝1000Hz。

从式(23)可知，仿真信号既含有高斯白噪声又含有平稳的非高斯噪声，高斯白噪声出现在调幅部分，而平稳的非高斯噪声作为整个信号的加性噪声。该实施例通过对概率主分量分析消噪前后的信号计算循环双谱，验证PPCA消噪的有效性和循环双谱诊断的准确性。

在总采样点中选择样本数m＝4076作为预处理的原始数据，并将其它分成20 段，即原始变量个数n＝20；设定主元变量个数k＝2，则原始数据被表示成X_20×4076的矩阵形式。概率主分量分析消噪前的原始信号波形和频谱图如图2所示。概率主分量分析消噪后的信号波形和频谱图如图3所示。对比图2(a)和图3(a)波形图可看出：概率主分量分析消噪能有效地抑制噪声。消噪前后信号的频谱图分别如图 2(b)、图3(b)所示。对比消噪前后信号的频谱图，同样可看出：P概率主分量分析能有效地消除噪声，使信号的频谱特征更加清晰和准确。

采用线性搜索的方法确定循环频率为了确定最优的循环频率α＝1000Hz，然后进行正弦抽取运算获得双谱分析的等高线图和三维图计算结果，如4所示。从图4 可知，六变形的顶点坐标分别为(0,-50)、(50,-50)、(50,0)、(0,50)、(-50,50)、(- 50,0)。这说明：对于循环平稳信号，中心频率位置的信噪比最高，不易受噪声干扰，而且中心频率与其它频谱成份关系相对固定，正好等于故障频率为50Hz，且稳定性最强。为验证概率主分量分析增强循环双谱效果，对原仿真信号直接计算循环双谱，且取循环频率仍为α＝1000Hz。其等高线图和三维图如图5所示。从图 5(a)可看出，由于平稳的非高斯噪声干扰，不能由故障频率值构成规则的循环双谱六边形结构。

具体实施例2：轴承内圈故障诊断

取某真实传动系统轴承内圈故障信号，实验用轴承的外圈带有故障，型号为 ER-12K，其中，滚动体个数为8个，滚动体直径为7.9375㎜，轴承的节圆直径为 33.4772mm，滚动轴承内圈故障特征频率理论计算结果197.1Hz。实验轴承空载运行，转速2990r/min，采样频率25.6kHz，采样点数327680，在采样点中选择样本数m＝327660作为预处理的原始数据，并将其分成20段，即原始变量个数 n＝20；设定主元变量个数k＝2，则原始数据被表示成X_20×327660的矩阵形式。

其原始信号波形及频谱图如图6所示，同样对原始数据信号先进行概率主分量分析消噪预处理，及消噪信号波形及频谱图如图7所示。从图6可知，有用信息完全被淹没在噪声中。而从图7可知，冲击成分突显存在，噪声得到一定的抑制，但还不能诊断出故障特征。

为了能进一步提取轴承内圈的故障特征信息，对概率主分量分析消噪后的信号计算循环双谱。采用线性搜索的方法确定循环频率为了确定最优的循环频率α＝2551Hz，然后进行正弦抽取运算获得双谱分析的等高线图和三维图计算结果，如8所示。从放大图8(b)中清楚看到围绕中心坐标(0,0)周围，生成一个规则的六边形结构，六个顶点分别由故障特征频率对组成，频率坐标分别为(0,-197)、(197,- 197)、(197,0)、(0,197)、(-197,197)、(-197,0)，说明轴承中心频率α＝2551Hz处信号调制频率为197Hz。调制频率与滚动轴承内圈故障特征频率理论计算结果 197.1Hz接近，从而诊断出轴承内圈故障。

为验证概率主分量分析增强循环双谱效果，对原始数据滚动轴承信号直接计算循环双谱，仍取循环频率为中心频率α＝2551Hz。其等高线图和三维图如图9所示。从图9(a)看出，由于平稳的非高斯噪声干扰，不能由故障频率值构成有规则的循环双谱六边形结构。

Claims (1)

1.一种利用概率主分量分析增强循环双谱的轴承故障诊断方法，其特征在于，包括以下步骤：①建立能反映原始信号与主元向量潜在联系的概率主元模型；②利用概率主元模型对原始信号进行消噪，获得高信噪比信号，完整地保留信号的有用成分，并大幅度地提高信噪比；③对消噪信号作循环双谱分析，通过计算三阶累积量、正弦抽取运算、二维FFT变换，从而获得单一循环频率双谱的等高线图，从单一循环频率双谱的等高线图中规则的六边形顶点坐标与故障特征频率比较，最终诊断出轴承故障；

其中，步骤①包括以下内容：

概率主元模型是在因子分析法隐含变量模型的基础上，假设高斯噪声的方差各向同性，由主元向量叠加噪声产生的生成模型，其数学形式为：

X＝P·u+E (1)

式(1)中，X＝{s₁,s₂,…,s_m}∈R^n×m是由原始信号s生成的n×m阶矩阵，称为原始变量，其中：n为嵌入维数，m为采样数；P＝{p₁,p₂,…,p_k}∈R^n×k为待求的n×k阶主元矩阵，k为主元个数，且k<n；u＝{u₁,u₂,…,u_m}∈R^k×m为主元向量，u_i(i＝1,2,…m)为主元向量的分量；E为服从N(0,σ²I)分布的高斯噪声，σ²为噪声方差；因此，X为服从N(0,PP^T+σ²I)的高斯分布；

确定概率模型后，分别计算主元向量、原始变量在给定主元向量、主元向量在给定原始变量下的分布，则可获得概率主元模型；具体计算过程如下：

假设主元向量u服从多元标准正态分布，其概率分布为：

原始变量在给定主元向量下的条件分布为：

原始变量的分布为：

式(4)中，C＝PP^T+σ²I是n×n阶原始变量的方差，根据贝叶斯定理，得到主元向量在给定原始变量下的条件分布:

式(5)中，M＝P^TP+σ²I，是已降维的k×k阶方阵；由式(2)～(5)可知，只要确定参数P和σ²，就可求得概率模型；由隐含模型经典参数估计EM算法估算参数：

式(6)和(7)中，为原始变量的协方差矩阵，tr(·)是矩阵的迹，即对主对角线上的所有元素取和，反复迭代式(6)和(7)直至收敛，即可求得主元矩阵P和高斯噪声的方差σ²；

步骤②包括以下内容：

根据建立的概率主元模型得到的主元矩阵P，将原始变量X分解为P的外积之和，得到原始信号的主元模型，其数学表示为：

式(8)中，i＝1,2,…,k，主元向量的对应的特征值反映了原始变量X在p_i方向的覆盖程度；将每一个主元向量按其特征值大到小排序，则其对应的主元矩阵就表示出原始变量X变化方向的大小，选择包含主要信息量的前i个主元矩阵重构数据，最终得到降维的原始变量X，实现信号消噪的目的；

概率主元模型对原始信号进行消噪取决于嵌入维数n和主元个数k；k可以由方差累积贡献率计算；单一主元贡献率计算公式如下：

式(9)中，λ为原始变量的协方差矩阵S的特征值；因此，全部主元贡献率计算公式为：

式(10)中，当达到85％时的最小k值，可确定主元个数；

当k确定了之后，由k<n确定最小的嵌入维数；

步骤③包括以下内容：

通过计算三阶累积量、正弦抽取运算、二维FFT变换，从而获得单一循环频率双谱的等高线图，具体为：

在频域中计算三阶累积量公式：

式(11)中，u(f₁)、u(f₂)表示主元向量u的二维傅里叶变换，E(·)表示求期望，u^*表示求共轭，α为引入的循环频率；由式(11)计算循环双谱可分为两个步骤：双谱估计；正弦抽取运算；

双谱估计采用非参数化法的间接法，所述间接法进行双谱估计的具体算法如下：

令u(1),u(2),...,u(m)是一组主元向量分量的采样信号，并设f_s是采样频率，N₀是总的频率采样点数，则Δf＝f_s/N₀是双谱区域沿水平和垂直方向上的频率采样间隔；

第1步：将主元向量u分成长度为m的k段，即u⁽ⁱ⁾(1),u⁽ⁱ⁾(2),...,u⁽ⁱ⁾(m)，其中i＝1,2,…k，并减去各段的均值；

第2步：设{u⁽ⁱ⁾(j),j＝1,2,...,m}是第i段数据，u⁽ⁱ⁾(j),i＝1,2,…,k；j＝1,2,...,m代表单个数据点，估计各段的三阶累积量：

式(12)中，m₁＝max(0,-k,-l),m₂＝min(m-1,m-1-k,m-1-l)；

第3步：计算所有段的三阶累积量的均值作为主元向量u的三阶累积量估计值

第4步：进行双谱估计：

式(14)中，L＜m-1,ω(k,l)是二维窗函数；

其次，进行正弦抽取运算；

从式(11)定义式可知，在频域定义的循环双谱相当于将双谱的频率向左平移循环频率α后得到；根据傅里叶移频定理，在时域，相当于将信号的三阶循环累积量乘以因子e^jαt后得到的，这就是正弦抽取运算，得到的循环双谱的统计量定义为：

式(15)中,T为周期，在时域中计算三阶累积量公式：

式(16)中，一、二、三阶矩分别为：

进行正弦抽取运算的关键是要预先估计循环频率α，它直接决定了双谱分析的成功与故障频率诊断的准确性；实际应用循环双谱时，采用线性搜索的方法确定循环频率：

第1步：在T＝[0,2π)内选择循环频率集{α}，以定步长Δα搜索，按式(16)计算每一循环频率α_j∈{α}的三阶循环累积量

第2步：对循环频率集{α}内计算的所有三阶累积量求和，得到

第3步：对进行二维傅里叶变换，得到属于循环频率区间(α₁,α₂)的循环双谱：

步骤④包括以下内容：根据获得单一循环频率双谱的等高线图，从单一循环频率双谱的等高线图中规则的六边形顶点坐标与故障特征频率比较，最终诊断出轴承故障。