CN105181336A - 一种用于轴承故障诊断的特征选取方法 - Google Patents

Abstract

本发明涉及一种用于轴承故障诊断的特征选取方法，本发明公布了故障诊断之前，轴承振动信号的奇异值分解(SVD)特征提取方法与基于特征区分度公式的特征选择方法，其实施步骤包括：(1)，不同故障轴承振动信号(加速度)采集；(2)，SVD得到包含全部信息的振动信号特征向量(奇异值向量)与自相关矩阵维数值；(3)，定义量化不同故障类型特征差异的区分度公式，评判方法的故障诊断能力；(4)，奇异值向量特征选取，提升后续诊断正确率，平衡诊断环境适应性与实时性要求。方法为轴承故障诊断提供了低维的，有明确区分的特征，在变工况负载等轴承故障诊断应用领域具有潜在经济和社会价值。

Description

一种用于轴承故障诊断的特征选取方法

技术领域

本发明涉及一种机械故障诊断方法，尤其是一种用于轴承故障诊断的特征选取方法。

背景技术

轴承是机械系统的重要部件之一，轴承故障诊断越来越受到重视，已成为保证生产安全稳定运行和提高产品质量的重要手段和关键技术。当前轴承故障特征提取方法大多依赖于某种(类)故障的先验知识，专门针对一种或几种故障，而由于故障类型、工况环境以及故障成因的多样性，其故障特征也各有差异，故面对复杂多变的故障环境，需要多种特征提取方法配合才有所成效，造成计算量与资源浪费，其根本原因是缺乏通用的，不随故障转移的特征提取方法。

同时轴承故障特征提取方法的偶然性导致了特征选择的不确定性，甚至无法实施特征选择过程，目前特征选择方法大多基于后续诊断的正确率效果反馈，以剔除弱故障区分能力的特征(或选择强故障区分能力的特征)，并重新建模。故特征对诊断效果先期估计与评判能力不足，且反复建模使得诊断时间代价较大，模型难以适应动态变化的参数和环境变化，缺乏诊断实时性。

发明内容

为了克服现有故障诊断特征提取和选择方法的不足，本发明提供了一种用于轴承故障诊断的特征选取方法。所述方法利用奇异值分解(SVD)进行故障信号特征提取，并定义特征区分度公式进行特征选择，具体采用以下技术方案：

一种用于轴承故障诊断的特征选取方法，包括如下步骤：

步骤(1)：采集轴承不同故障类型的多个振动信号；所述不同故障类型包括不同故障直径和不同故障部位；根据所采集的振动信号构建各个振动信号自相关函数的矩阵，并对所构建的矩阵进行奇异值分解，得到最大奇异值元素；

步骤(2)：逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定；得到此时各个振动信号的奇异值向量和矩阵维数值；

步骤(3)：根据步骤(2)得到的各个振动信号的奇异值向量和矩阵维数值得到不同维数下各个振动信号的奇异值向量的平均曲率；根据各个振动信号的奇异值向量的平均曲率定义反映不同类型的轴承振动信号差异的区分度公式：

$d_{ab}=q\left|ln\left(\frac{k_a}{k_b}\right)\right|$

其中，d_ab为第a类故障类型和第b类故障类型的特征区分度(a∈{1,2,…,r}；b∈{1,2,…,r})，r为故障类型总数；k_a和k_b为这两类故障特征向量曲线的平均曲率；q为分类性能指标参数统计值；

步骤(4)：逐步减小矩阵维数，计算矩阵奇异值向量，并根据步骤(3)定义的区分度公式计算不同类型的轴承振动信号之间的区分度；当所有轴承振动信号之间的区分度最小值为1时，此时的矩阵维数即为特征选取维数，同时得到的该维数下的奇异值向量即为选取的特征向量。

所述步骤(1)采集的振动信号：

y_id＝{x_i},(i＝1,2,…,k)

其中，y_id为第id(id＝1,2,3…,r)类故障的轴承振动信号序列；r为总故障类别数；x_i＝x(t_i)为t_i时刻振动加速度数据；t_i＝t₀+(i-1)Δt(i＝1,2,…,k)；t₀为初始采样时刻，Δt＝1/f为相邻采样时间间隔；f为采样频率；k为总采样数；i为采样计数。

所述步骤(2)中，逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定过程如下；

步骤(21)：周期扩展振动信号y_id，并计算扩展信号的自相关函数 $c\left(m\right)={\mathrm{\Σ}}_{m=1}^{\mathrm{\∞}}{y}_{id}\left(m\right)y_{id}(m+n),(n=0,1,...\mathrm{\∞});$ n为扩展后的振动信号第n个序列点，m为计算自相关函数时的平移点数；选取较小的维数构造初始自相关函数Hankel矩阵C₀：

$C_0={\left[\begin{array}{cccc}c\left(0\right)& c\left(1\right)& ...& c(N_0-1)\\ c\left(1\right)& c\left(2\right)& ...& c\left(N_0\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ c(M_0-1)& c\left(M_0\right)& ...& c(N_0+M_0)\end{array}\right]}_{M_0\×N_0}$

其中，C₀为初始自相关函数Hankel矩阵，M₀>0，N₀>0为二维矩阵的初始选取维数；

步骤(22)，对C₀进行奇异值分解：C₀＝U₀D₀V₀ ^T；其中，矩阵 为M₀×M₀维的正交矩阵，为矩阵元素向量； 为N₀×N₀维的正交矩阵，为矩阵元素向量；D₀为对角阵：

$D_0=\[diag\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,......,{\mathrm{\λ}}_{Q_0}\},0\],(Q_0=min\{M_0,N_0\left\}\right)$

其中，为自相关矩阵C₀的奇异值，且故 $C_0={\mathrm{\λ}}_1{u}_1{v_1}^T+{\mathrm{\λ}}_2{u}_2{v_2}^T+...+{\mathrm{\λ}}_{Q_0}{u}_{Q_0}{v_{Q_0}}^T.$

步骤(23)，令M_t＝M₀+t，N_s＝N₀+s，(t＝1,2,3…，s＝1,2,3…),Q_p＝min{M_t,N_s}，p＝min{t,s}；其中，t与s为矩阵横向和矩阵纵向的扩展维数，每次交替增加1进行扩展；重复步骤(21)和步骤(22)，直到得到最大奇异值元素λ₁达到稳定；此时奇异值向量 $f=({\mathrm{\λ}}_1,{\mathrm{\λ}}_2,......,{\mathrm{\λ}}_{Q_t})$ 为提取的特征向量。

所述步骤(23)中得到最大奇异值元素λ₁的条件为若维数Q_p时最大奇异值元素等于维数Q_p+1时的最大奇异值元素，或两者之差小于0.01λ₁值，此时Q_p的最大奇异值元素为稳定的最大奇异值元素；其中Q_p是指扩展后的维数；Q_p+1是指Q_p再次扩展后的维数。

本发明的有益效果是：本发明提出了一种用于轴承故障诊断的特征提取和选择方法，所述方法利用奇异值分解(SVD)进行故障信号特征提取，能够有效克服依赖故障特征先验知识的局限，并提升轴承故障诊断通用性；所述方法定义特征区分度公式，能够有效评判方法的故障诊断能力；特征选取能够有效提升后续诊断正确率，平衡诊断环境适应性与实时性要求。本发明为轴承故障诊断提供了低维的，有明确区分的特征，在变工况负载等轴承故障诊断应用领域具有潜在经济和社会价值。

附图说明

图1为本发明的用于轴承故障诊断的特征提取和选择算法流程图；

图2为本发明的目标振动数据与辅助振动数采集示意图(a：故障部位；b：故障直径)；

图3为本发明的轴承振动信号示例图(健康轴承，12kHz，驱动端对象)；

图4为本发明的最大特征值与自相关矩阵维数关系曲线(故障类型：轻微，严重)；

图5为本发明的SVD提取奇异值向量的对数曲线图(a：故障直径；b：故障部位)；

图6为本发明的奇异值特征向量的平均曲率曲线(a：电机转速；b：故障直径；c：故障部位；d：试验对象)；

图7为本发明的区分度与分类性能指标参数比值曲线(故障直径)；

图8为小波变换特征提取与选择方法的特征频率的信号能量曲线(a：145.83Hz；b：157.94Hz；c：315.89Hz；d：473.83Hz)。

具体实施方式

下面结合附图及具体实施例对本发明进行详细说明。

参见图1，一种用于轴承故障诊断的特征提取和选择方法，包括如下步骤：

步骤(1)：对含不同故障类型的轴承进行振动信号(加速度)采集，用于特征提取、选择以及故障区分性能测试；

步骤(1)中，不同故障类型，包括不同故障直径与不同故障部位的信号采集，其中轴承振动信号采集示意图如图2所示：

a)例如，故障直径包括：0.000mm、0.178mm、0.356mm、0.533mm直径，分别定义为健康、轻微故障、中度故障、严重故障；

b)例如，故障部位包括：内圈、滚珠、外圈6:00、外圈3:00、外圈12:00(负载6:00方向)；

c)实测中，采集频率f＝12kHz或48kHz条件的振动信号，在本发明中，采集频率为轴承转速的200倍以上。

y_id＝{x_i},(i＝1,2,…,k)

其中，y_id为第id(id＝1,2,3…,r)类故障的含观测噪音轴承振动数据序列；r为总故障类别数；x_i＝x(t_i)为t_i时刻振动加速度数据；t_i＝t₀+(i-1)Δt(i＝1,2,…,k)；t₀为初始采样时刻，Δt＝1/f为相邻采样时间间隔；k为总采样数；i为采样计数。

轴承振动信号示例图如图3所示，可以看出，在变化工作环境下，振动信号呈现周期非平稳特性，并且混有高斯噪声，使得轴承的故障特征易被正常振动引起干扰信号以及白噪声掩盖，增大提取轴承故障特征提取的复杂性和困难度。

步骤(2)：利用Hankel矩阵构建振动信号自相关函数的自相关矩阵C，增大矩阵维数，观察其奇异值分解(SVD)得到的最大特征值元素变化，直至保持恒定，并得到稳定时的振动信号特征向量(奇异值向量)与矩阵维数值；

步骤(2)中，由于轴承故障发生时，损伤点与其它表面元件碰撞产生信号冲击，可能覆盖轴承系统和传感器的固有频率，故自相关矩阵奇异值分解的实施目的是，克服大多轴承故障诊断方法依赖故障特征先验知识(故障频率、噪音干扰等)的局限，提供不依赖先验知识，并具备物理意义的清晰特征，提升轴承故障诊断通用性。

步骤(2)中，自相关矩阵C的奇异值分解步骤如下：

步骤(21)，周期扩展振动信号y_id，并计算扩展信号的自相关函数 $c\left(m\right)={\mathrm{\Σ}}_{m=1}^{\mathrm{\∞}}{y}_{id}\left(m\right)y_{id}(m+n),(n=0,1,...\mathrm{\∞});$ n为扩展后的振动信号第n个序列点，m为计算自相关函数时的平移点数；选取较小的维数构造初始自相关函数Hankel矩阵C₀：

$C_0={\left[\begin{array}{cccc}c\left(0\right)& c\left(1\right)& ...& c(N_0-1)\\ c\left(1\right)& c\left(2\right)& ...& c\left(N_0\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ c(M_0-1)& c\left(M_0\right)& ...& c(N_0+M_0)\end{array}\right]}_{M_0\×N_0}$

其中，C₀为初始自相关函数Hankel矩阵，M₀>0，N₀>0为二维矩阵的初始选取维数；

步骤(22)，对C₀进行奇异值分解：C₀＝U₀D₀V₀ ^T；其中，矩阵 为M₀×M₀维的正交矩阵，为矩阵元素向量； 为N₀×N₀维的正交矩阵，为矩阵元素向量；D₀为对角阵：

$D_0=\[diag\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,......,{\mathrm{\λ}}_{Q_0}\},0\],(Q_0=min\{M_0,N_0\left\}\right)$

其中，为自相关矩阵C₀的奇异值，且故 $C_0={\mathrm{\λ}}_1{u}_1{v_1}^T+{\mathrm{\λ}}_2{u}_2{v_2}^T+...+{\mathrm{\λ}}_{Q_0}{u}_{Q_0}{v_{Q_0}}^T.$

步骤(23)，令M_t＝M₀+t，N_s＝N₀+s，(t＝1,2,3…，s＝1,2,3…),Q_p＝min{M_t,N_s}，p＝min{t,s}；利用M_t和N_s维数重复步骤B-1和B-2，得到最大奇异值元素λ₁，若维数Q_p时最大奇异值元素等于(或两者之差小于0.01λ₁值)维数Q_p+1时最大奇异值元素，此时Q_p为包含全部有效信息的维数，此时奇异值向量为提取特征。

参见图4和图5，图4为最大特征值与自相关矩阵维数关系曲线(故障类型：轻微，严重)，可以看出，不同故障程度曲线的无冗余奇异值维数有差异(轻微故障：137维，严重故障：152维)，增加维数对有效信息量提升较少(突变点后曲线变化较为平缓)，可认为是冗余信息。进一步研究表明，不同故障类型的有效维数在127～164之间随机波动，说明有效故障特征与冗余噪音的区分边界较为模糊，Hankel矩阵维数选取M_t＝N_s＝Q_p＝150。图5为本发明的SVD提取奇异值向量的对数曲线图(a：故障直径；b：故障部位)；可见不同故障类型(a：故障直径，b：故障部位)，其特征曲线前端具有明显差异性，所提特征能有效实施故障分类。

步骤(3)：根据步骤(2)得到的各个振动信号的奇异值向量和矩阵维数值得到不同维数下各个振动信号的奇异值向量的平均曲率；并根据各个振动信号的奇异值向量的平均曲率定义反映不同类型的轴承振动信号差异的区分度公式，并利用后续故障诊断的分类性能统计指标，矫正平衡区分度公式；

步骤(3)中，区分度定义目的是，反映不同轴承故障类型特征之间差异程度，其表征两类故障特征间的可分类(可诊断)能力。

步骤(3)中，所述区分度公式为：

$d_{ab}=q\left|ln\left(\frac{k_a}{k_b}\right)\right|$

其中，d_ab为第a类故障类型和第b类故障类型的特征区分度(a∈{1,2,…,r}；b∈{1,2,…,r})，r为故障类型总数；k_a和k_b为这两类故障特征向量曲线的平均曲率；q为分类性能指标参数统计值，使得恰好达到诊断目标时，d_ab＝1，可知诊断正确率目标越高，q值越小，针对不同的诊断目的(故障直径和故障部位)，q的统计值也不同。由上述定义可知|ln(k_a/k_b)|≥0，且|ln(k₁/k₂)|越大，曲线特征差异越大。不同的故障分类目的和指标，其分类性能具有差异，即使相同的|ln(k₁/k₂)|值，对故障直径和故障部位的区分也具有差异，同时90％分类正确率指标与95％的分类正确率指标对平均曲率要求也不同。

参见图6，图6为本发明的奇异值特征向量的平均曲率曲线(a：电机转速；b：故障直径；c：故障部位；d：试验对象)，其中不同轴承故障对应点间的平均曲率间距反映特征区分能力，故有：1)故障类型的量变与质变间存在差异，如图6(b)中，健康轴承(0mm)与故障轴承的平均曲率间距大于故障轴承内部间距；2)故障类别的本质存在差异，图6(c)中，三类外圈故障间平均曲率间距小于外圈与其它故障部位间距；同时奇异值特征向量的平均曲率随维数递减，即维数Q越小，轴承故障分类能力越强，不过维数越小，包含的有用信息量越少。

下表给出了故障直径诊断目标的分类性能指标参数q列表，可见诊断目标精度越高，统计分类性能指标参数越小，例如在轴承故障直径诊断条件下，若使的目标诊断统计正确率达99.9％，d_ab/q值(即|ln(k_a/k_b)|)应大于0.893(1/1.12)；同理，若正确率达90.0％，d_ab/q值只需大于0.282。

根据两类区分度公式构造轴承故障诊断区分度的距离矩阵D：

$D={\left[\begin{array}{cccc}{d}_{11}& d_{12}& ...& d_{1r}\\ d_{21}& d_{22}& ...& d_{2r}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ d_{r1}& d_{r2}& ...& d_{rr}\end{array}\right]}_{r\×r}$

其中，d_ab(a∈{1,2,…,r}；b∈{1,2,…,r})值越大，特征提取方法对第a类故障和第b类故障区分能力越强，且d_ab＝d_ba，d_aa＝0；当D为零矩阵时，特征提取方法对任何两类故障均无区分能力；对某一特征提取方法，由于故障标识是人为定义的，(例如，健康和轻微被定义为1和2，也可定义为2和1)，故D中元素具有可交换性和等价性，D中元素均值反映特征提取方法的优劣程度，均值越大，特征提取方法越优。

步骤(4)：逐步减小矩阵维数，即从奇异值向量末端逐步减少维数，实施特征选取，使得特征包含信息量、后续分类性能和诊断计算量得到最佳平衡。

步骤(4)中，所述轴承故障特征选取的实质是降低奇异值向量维数Q_p，逐步减少特征向量末端重叠较严重的维数并计算D值，当D中最小元素为1时，即所有轴承振动信号之间的区分度最小值为1时，此时奇异值向量维数Q_c(Q_c≤Q_p)即特征选取维数，同时得到的该维数下的奇异值向量即为选取的特征向量。特征选取有助于提升分类(诊断)性能、降低诊断计算量，但一定程度上会损失故障信息量。

参见图7，图7所示为本发明的区分度与分类性能指标参数比值曲线(故障直径)；图中可以看出，对于99.9％目标诊断正确率指标，根据上表，要求区分度与分类性能指标参数比值大于0.893，图7中2组曲线满足要求，3组曲线不能满足该要求(即无论维数如何选取，诊断效果无法达到99.9％)，此外，轻微故障和中度故障(0.178mm和0.356mm故障直径)的区分性能在维数低于96时达到目标效果，故特征选取维数为96，此时兼顾诊断正确率、包含信息量以及后续计算量的要求。

参见图8，实测中，采用小波变换的特征提取与选择方法与本发明的SVD特征提取与选择方法比较，图8所示为小波变换特征提取与选择方法的特征频率的信号能量曲线(a：145.83Hz；b：157.94Hz；c：315.89Hz；d：473.83Hz)；可以看出，在特征频率±1.5Hz范围内，存在反映该频率的能量峰值，相同故障类型不同故障程度的能量具有差异性(①和②点)，反映故障程度区分特征，不同故障类型能量更有区分性(即是否存在能量峰值)。但在未知先验特征频率情况下，即使多层小波分解放大了信号的微弱细节，其仍无法判定故障特征，一般而言，放大的仅仅为具有相同统计特征的随机信号，故该方法不具备区分性，并且不同故障类型的先验特征频率往往不同(需要事先采用其它方法获知)，故缺乏通用性。但是相比小波变换的4维特征向量而言，SVD方法的96维特征向量的后续诊断计算量大幅增加，其通用性的代价是牺牲诊断实时性。

Claims (4)

1.一种用于轴承故障诊断的特征选取方法，其特征在于，包括如下步骤：

步骤(1)：采集轴承不同故障类型的多个振动信号；所述不同故障类型包括不同故障直径和不同故障部位；根据所采集的振动信号构建各个振动信号自相关函数的矩阵，并对所构建的矩阵进行奇异值分解，得到最大奇异值元素；

步骤(2)：逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定；得到此时各个振动信号的奇异值向量和矩阵维数值；

步骤(3)：根据步骤(2)得到的各个振动信号的奇异值向量和矩阵维数值得到不同维数下各个振动信号的奇异值向量的平均曲率；根据各个振动信号的奇异值向量的平均曲率定义反映不同类型的轴承振动信号差异的区分度公式：

$d_{ab}=q\left|ln\left(\frac{k_a}{k_b}\right)\right|$

其中，d_ab为第a类故障类型和第b类故障类型的特征区分度(a∈*1,2,…,r+；b∈*1,2,…,r+)，r为故障类型总数；k_a和k_b为这两类故障特征向量曲线的平均曲率；q为分类性能指标参数统计值；

步骤(4)：逐步减小矩阵维数，计算矩阵奇异值向量，并根据步骤(3)定义的区分度公式计算不同类型的轴承振动信号之间的区分度；当所有轴承振动信号之间的区分度最小值为1时，此时的矩阵维数即为特征选取维数，同时得到的该维数下的奇异值向量即为选取的特征向量。

2.根据权利要求1所述的特征选取方法，其特征在于，所述步骤(1)采集的振动信号：

y_id＝*x_i+,(i＝1,2,…,k)

其中，y_id为第id(id＝1,2,3…,r)类故障的轴承振动信号序列；r为总故障类别数；x_i＝x(t_i)为t_i时刻振动加速度数据；t_i＝t₀+(i-1)Δt(i＝1,2,…,k)；t₀为初始采样时刻，Δt＝1/f为相邻采样时间间隔；f为采样频率；k为总采样数；i为采样计数。

3.根据权利要求1所述的特征选取方法，其特征在于，所述步骤(2)中，逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定过程如下；

步骤(21)：周期扩展振动信号y_id，并计算扩展信号的自相关函数 $c\left(m\right)={\mathrm{\Σ}}_{m=1}^{\mathrm{\∞}}{y}_{id}\left(m\right)y_{id}(m+n),(n=0,1,...\mathrm{\∞});$ n为扩展后的振动信号第n个序列点，m为计算自相关函数时的平移点数；选取较小的维数构造初始自相关函数Hankel矩阵C₀：

$C_0={\left[\begin{array}{cccc}c\left(0\right)& c\left(1\right)& ...& c(N_0-1)\\ c\left(1\right)& c\left(2\right)& ...& c\left(N_0\right)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ c(M_0-1)& c\left(M_0\right)& ...& c(N_0+M_0)\end{array}\right]}_{M_0\×N_0}$

其中，C₀为初始自相关函数Hankel矩阵，M₀>0，N₀>0为二维矩阵的初始选取维数；

步骤(22)，对C₀进行奇异值分解：C₀＝U₀D₀V₀ ^T；其中，矩阵 为M₀×M₀维的正交矩阵，为矩阵元素向量； 为N₀×N₀维的正交矩阵，为矩阵元素向量；D₀为对角阵：

$D_0=\[diag\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,......,{\mathrm{\λ}}_{Q_0}\},0\],(Q_0=min\{M_0,N_0\left\}\right)$

其中，为自相关矩阵C₀的奇异值，且故 $C_0={\mathrm{\λ}}_1{u}_1{v_1}^T+{\mathrm{\λ}}_2{u}_2{v_2}^T+...+{\mathrm{\λ}}_{Q_0}{u}_{Q_0}{v_{Q_0}}^T;$

步骤(23)，令M_t＝M₀+t，N_s＝N₀+s，(t＝1,2,3…，s＝1,2,3…),Q_p＝min*M_t,N_s+，p＝min*t,s+；其中，t与s为矩阵横向和矩阵纵向的扩展维数，每次交替增加1进行扩展；重复步骤(21)和步骤(22)，直到得到最大奇异值元素λ₁达到稳定；此时奇异值向量为提取的特征向量。

4.根据权利要求3所述的特征选取方法，其特征在于，所述步骤(23)中得到最大奇异值元素λ₁的条件为若维数Q_p时最大奇异值元素等于维数Q_p+1时的最大奇异值元素，或两者之差小于0.01λ₁值，此时Q_p的最大奇异值元素为稳定的最大奇异值元素；其中Q_p是指扩展后的维数；Q_p+1是指Q_p再次扩展后的维数。