CN103868692B - 基于核密度估计和k-l散度的旋转机械故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于核密度估计和K-L散度的旋转机械故障诊断方法，包括如下步骤：采集被监测对象的原始振动数据，并划定训练样本集和测试样本集；从步骤1中得到原始振动数据中提取指定时频域特征；选择出少数敏感特征，并且计算这些敏感特征的分类贡献率；利用核密度估计计算训练样本中不同故障类别样本集关于敏感特征的概率密度函数，并计算加入一个未知故障类别待测样本后各类样本集新的概率密度函数；计算出在选定的特征描述下，训练样本中各类故障样本集原始概率密度函数，以及加入一个待测样本后的新概率密度函数两者的K-L散度值；计算集成K-L散度，并通过集成K-L散度的大小判断待测样本的故障类别。从而提高分类器的准确率和推广能力。

Description

基于核密度估计和K-L散度的旋转机械故障诊断方法

技术领域

本发明属于机械装备故障智能诊断领域，具体涉及基于核密度估计和K-L散度等统计学工具的故障诊断方法。

背景技术

旋转机械广泛应用于工业生产实践中，例如风力发电机、数控机床、航空航天发动机等关系到国防民生的重要领域。在生产工作中，滚动轴承、齿轮等旋转机械中的关键部件由于需要承受交变机械应力和偶然冲击，加上本身固有的制造误差，经常会产生一些早期缺陷，例如轻度磨损、点蚀等。这些缺陷如果不及时诊断发现，就会不断恶化，最终导致***失效，带来很大的财产损失，甚至对国防和人身安全带来巨大威胁。对于现代化大型复杂旋转机械设备而言，尽管可以通过改善设计、制造工艺来提高零部件的质量，但仍难以确保不出故障。所以，有必要利用先进的传感和监测技术，对关键零部件及***进行有计划、有组织、有针对性的状态监测和故障诊断，尽早发现设备运行过程中的各种隐患，从而防止***损失和灾难性事故的发生。

数据驱动方法是近年来逐渐兴起的一种故障诊断技术，计算机技术的快速发展使得大数据并行高速计算变得非常容易，推动了依靠大量数据分析的故障诊断技术的发展。从应用的角度看，基于数据驱动的故障诊断方法较基于模型的方法更为切实可行，这是由于数据采集通常要比精确建立物理模型更加容易。除此之外，数据驱动的故障诊断方法还有两个明显的优点：一是该类方法更容易实现自动诊断，这与现代工业的智能化发展是切合的；二是该类方法不需要太多参数设置和专家经验知识。一般来说，一种数据驱动的故障诊断方法应包括数据获取、特征提取、特征降维、分类器设计和结果输出等五个步骤，其中分类器设计和选择是该类方法的关键。

现有的数据驱动的故障诊断方法大多都是在样本的特征空间寻找一个最优分类超曲面，从而将不同类型的故障样本分开。例如，基于支持向量机（SupportVectorMachine，SVM）的故障诊断方法，基于BP（BackPropagation）神经网络的故障诊断方法等。然而，由于噪声、测量误差等对有效振动信号的污染，使得分类问题存在一定的不确定性，进而导致了分类错误的现象，即很难找到一个恰当的超曲面将所有样本全部归类正确。

以上常用方法产生分类错误的现象的原因与它们的分类原理是分不开的。传统的智能故障诊断方法往往忽略了样本间的统计信息和关联信息，而统计信息对于随机信号处理是极为关键的，即对于正确分类很有帮助。目前，国内外从样本统计角度开展的智能故障诊断方法的研究或报道还非常少。

发明内容

本发明的目的是为了从原始样本中提取更加全面、有效的统计学信息，从而提高分类器的准确率和推广能力，提出一种基于核密度估计（KernelDensityEstimation，KDE）和K-L散度（Kullback-LeiblerDivergence）两种统计学工具的智能故障诊断方法。

本发明的基于核密度估计和K-L散度的旋转机械故障诊断方法，包括如下步骤：

步骤1：采集被监测对象的原始振动数据，并划定训练样本集和测试样本集；

步骤2：从步骤1中得到的原始振动数据中提取指定时频域特征；

步骤3：从步骤2中得到的频域特征集中选择出少数敏感特征，并且计算这些敏感特征的分类贡献率；

步骤4：利用核密度估计计算训练样本中不同故障类别样本集关于步骤3中提取的敏感特征的概率密度函数，并计算加入一个未知故障类别待测样本后各类样本集新的概率密度函数；

步骤5：计算出在选定的特征描述下，训练样本中各类故障样本集原始概率密度函数，以及加入一个待测样本后的新概率密度函数两者的K-L散度值；

步骤6：计算集成K-L散度，并通过集成K-L散度的大小判断待测样本的故障类别。

进一步地，所述步骤2中得到的指定时频域特征是通过总体平均经验模态分解方法和希尔伯特变换等信号处理方法得到的。

进一步地，所述步骤3的计算过程如下：

第一步：计算第j个特征C个类的类内距离的平均值

$\left\{\begin{array}{c}{d}_{c,j}=\frac{1}{M_c\×(M_c-1)}\underset{l,m=1}{\overset{M_c}{\mathrm{\Σ}}}|q_{m,c,j}-q_{l,c,j}|\\ l,m=\mathrm{1,2},\·\·\·,M_c,l\≠m,j=\mathrm{1,2},\·\·\·,J,c=\mathrm{1,2},\·\·\·,C\end{array}\right.---\left(1\right)$

$d_j^{\left(w\right)}=\frac{1}{C}\underset{c=1}{\overset{C}{\mathrm{\Σ}}}{d}_{c,j}---\left(2\right)$

其中，Mc表示第c类的样本个数，J表示特征个数，C表示类别个数，q_m,c,j表示第c类第m个样本的第j个特征的特征值；

第二步：计算第j个特征C个类的类间距离的平均值

$u_{c,j}=\frac{1}{M_c}\underset{m=1}{\overset{M_c}{\mathrm{\Σ}}}{q}_{m,c,j}---\left(3\right)$

$\left\{\begin{array}{c}{d}_j^{\left(b\right)}=\frac{1}{C\×(C-1)}\underset{c,e=1}{\overset{C}{\mathrm{\Σ}}}|u_{e,j}-u_{c,j}|\\ c,e=\mathrm{1,2},\·\·\·,C,c\≠e\end{array}\right.---\left(4\right)$

其中，u_c,j、u_e,j分别表示第c和第e个类第j个特征的平均值；

第三步：数据集A的类间距离与类内距离的比值α_j：

${\mathrm{\α}}_j=\frac{d_j^{\left(b\right)}}{d_j^{\left(w\right)}}---\left(5\right)$

α_j越大表示第j个特征对于类别越敏感，更符合评估原则，也是更应该被选择使用的特征。

第四步：计算选择的前n个特征的分类贡献率F_j；定义：

$F_j=\frac{{\mathrm{\α}}_j}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i},(j=\mathrm{1,2},\·\·\·,n)---\left(6\right)$

式（6）为分类贡献率。

进一步地，所述核密度估计是概率论中一种用来估计随机变量概率密度函数的非参数方法，设X₁,X₂,…,X_n是取自一元连续总体X的样本，在任意点x处的总体密度函数f(x)的核密度估计定义为：

${\hat{f}}_h\left(x\right)=\frac{1}{\mathrm{nh}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}K\left(\frac{x-x_i}{h}\right)---\left(7\right)$

其中，K(·)称为核函数，h为窗宽。

进一步地，所述的窗宽h是核密度估计方法中唯一要优化的参数，根据最小化平均集成平方误差方法，得到下式：

$h={\left(\frac{4{\widehat{\mathrm{\σ}}}^5}{3n}\right)}^{\frac{1}{5}}\≈{1.06\widehat{\mathrm{\σ}}n}^{-\frac{1}{5}}---\left(8\right)$

其中，是样本标准差，n为各类样本集所包含的样本数。

进一步地，所述步骤5的运算过程如下：

K-L散度也称为相对熵或者信息增益。在概率论和信息论中，K-L散度用来计算两个分布的分对称性或者差异性。K-L散度值越小，表明两个分布越相似。

原始的K-L散度定义式为：

$D_{\mathrm{KL}}\left(P\right|\left|Q\right)=\underset{i}{\mathrm{\Σ}}\ln \left(\frac{P\left(i\right)}{Q\left(i\right)}\right)P\left(i\right)---\left(9\right)$

为了满足距离的对称性，重新定义为：

$D_{\mathrm{KL}}(P,Q)=\frac{1}{2}(D_{\mathrm{KL}}\left(P\right|\left|Q\right)+D_{\mathrm{KL}}\left(Q\right|\left|P\right))---\left(10\right)$

由步骤4得到各类样本集的原始密度函数和加入测试样本后的新密度函数之后，就可以根据式（10）计算出对应的K-L散度值。

进一步地，所述步骤6的计算过程如下：

定义了一个集成K-L散度值IKL_i（i=1,2,…,C）：

${\mathrm{IKL}}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{F}_j\×{\mathrm{KL}}_i^j---\left(11\right)$

其中，F_j为步骤3得到的前n个特征分类贡献率，式（11）中通过将第j个特征描述下的K-L散度值与第j个特征的分类贡献率F_j加权，可以得到各类样本集的集成K-L散度值。IKL_i越小，原始分布和加入待测样本后的新概率密度函数越相似；反之，当集IKL_i比较大时，原始概率密度函数与新概率密度函数差别越大。换句话说，待测样本应该归类于K-L散度值最小的类别。

本发明的有益效果：由于使用了核密度估计（KernelDensityEstimation，KDE）和K-L散度（Kullback-LeiblerDivergence）两种统计学工具，所以充分考虑到样本的不确定性，进而提高了分类器的准确率和推广能力。

附图说明

图1为本发明的基于核密度估计和K-L散度的旋转机械故障诊断方法的流程示意图；

图2为本发明实例所针对的来自凯斯西储大学的滚动轴承故障诊断平台；

图3数据集A所有特征类间/类内距离比值α_j；

图4基于核密度估计和K-L散度的分类器原理示意图；

图5将本发明用于表1中数据集A的实验结果：四类样本集原始概率密度曲线及对应的加入一个正常样本后新的概率密度函数；

图6为不同特征数量下，本发明方法在数据集A-E上的分类准确率；

图7为不同训练样本下，本发明方法在数据集A-E上的分类准确率。

具体实施方式

下面结合附图和具体的实施例对本发明作进一步的阐述。

如图1所示，其步骤包括：

步骤1：采集被监测对象的原始振动数据，并划定训练样本集和测试样本集。

本发明以凯斯西储大学的滚动轴承故障诊断平台为例，具体的实验参数如下：

如附图2所示，实验平台包括一个2马力的电机（左侧），一个功率计（右侧）和电子控制设备（没有显示）。被测轴承起到支撑电机轴的作用。利用电火花加工技术在驱动端轴承上人为制造了单点故障，故障直径分别为0.007、0.014、0.021和0.028英寸。这些故障单独分布在轴承内圈、滚动体和外圈上。电机转轴上加载了一个冲击力，使用两个加速度传感器测量振动，一个传感器布置在电机壳体上，另一个传感器布置在驱动端轴承外圈上。采样频率为12KHz，单采样样本长度为12000。

为了验证上文所提方法的有效性，我们进行了两种实验，第一种是不同故障类型判别实验，第二种是同故障损伤程度判别实验，详细信息如表1和表2所示。

数据集A包括轴承4种不同的运行状态（正常、外圈故障、内圈故障、滚动体故障）的280个数据样本，损伤尺寸为0.007英寸。数据集A分为两部分，140个样本用来训练，140个用来测试，很明显这是一个四种不同健康状态的模式识别问题。

数据集B同样包含280个数据样本，其中训练样本与测试样本仍然各占一半，包含内圈和滚动体两种故障类型。数据集B又包含B₁和B₂两个数据子集，每个包含140个样本。数据子集B₁包含70个训练样本，每个样本的损伤程度为0.007英寸，同时包含损伤程度为0.021英寸的70个测试样本。数据子集B₂与B₁类似，只是将B₁的训练样本和测试样本进行交换，即B₁的测试样本作为B₂的训练样本。使用数据集B的作用是验证该方法对于同种故障下、不同于训练样本损伤程度的测试样本的检测能力，即验证该方法的稳健性。

表1故障类型判别实验数据集

数据集C、D和E各包含210个样本，训练和测试样本各占一半。数据集C、D和E的故障类型分别为内圈故障、滚动体故障和外圈故障。每种故障损伤程度分为三个层次，0.007、0.021和0.028英寸。使用数据C、D和E的目的是验证该方法对于同种故障不同损伤程度的分辨能力。

表2故障损伤程度判别实验数据集

为了精简篇幅，这里仅以数据集A为例叙述本发明方法的后续步骤。

步骤2：利用EEMD和希尔伯特变换等信号处理方法从原始振动数据中提取指定时频域特征。

利用EEMD对数据集A所对应的振动信号样本进行分解，白噪声初始值为0.3，集成数为100，于是可以得到每个样本的前四阶IMF（固有模态分量）。接下来，分别从数据集A原始振动信号和四阶IMF分量中提取制定的9个时域特征和10个频域特征，从而得到每类样本集所对应的95个原始特征。所提取的时域特征和频域特征种类如表3和表4所示。

表3时域特征表

表4频域特征表

步骤3：利用基于距离的特征评估方法从步骤2中得到的特征集中选择出前10个敏感特征，并且计算这些敏感特征的分类贡献率，以供后续步骤使用。应用基于距离的特征评估方法包括以下四步：

（1）计算数据集A中四类样本集在第j个特征描述下的类内距离的平均值

数据集A中共包含四类故障，各类样本数为35，总的特征数为95，所以Mc=35，C=4，J=95；将这些参数带入（1）和（2）式，有：

$\left\{\begin{array}{c}{d}_{c,j}=\frac{1}{M_c\×(M_c-1)}\underset{l,m=1}{\overset{M_c}{\mathrm{\Σ}}}|q_{m,c,j}-q_{l,c,j}|\\ l,m=\mathrm{1,2},\·\·\·,M_c,l\≠m,j=\mathrm{1,2},\·\·\·,J,c=\mathrm{1,2},\·\·\·,C\end{array}\right.$

$d_j^{\left(w\right)}=\frac{1}{4}\underset{c=1}{\overset{4}{\mathrm{\Σ}}}{d}_{c,j}$

式中，q_m,c,j表示第c类第m个样本的第j个特征的特征值。

（2）计算数据集A中四类样本集在第j个特征描述下的类间距离的平均值

同上步骤，将相关参数带入（3）和（4）式，得：

$u_{c,j}=\frac{1}{35}\underset{m=1}{\overset{35}{\mathrm{\Σ}}}{q}_{m,c,j}$

$d_j^{\left(b\right)}=\frac{1}{4\×(4-1)}\underset{c=1}{\overset{4}{\mathrm{\Σ}}}\underset{\underset{e\≠c}{e=1}}{\overset{4}{\mathrm{\Σ}}}|u_{e,j}-u_{c,j}|$

式中，u_c,j、u_e,j分别表示第c和第e个类第j个特征的平均值。

（3）数据集A的类间距离与类内距离的比值α_j：

${\mathrm{\α}}_j=\frac{d_j^{\left(b\right)}}{d_j^{\left(w\right)}},j=\mathrm{1,2},\·\·\·,95$

95个特征的α_j值如附图3所示。

（4）在本实验中，我们选取了α_j值排名前10的特征作为敏感特征，并计算得到了这10个特征的分类贡献率。这10个特征的编号和α_j值如表5所示。进一步计算被选择的这10个敏感特征的分类贡献率F_j。

对于第j个特征，有：

$F_j=\frac{{\mathrm{\α}}_j}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i},(j=\mathrm{1,2},\·\·\·,10)$

表5敏感特征的类间、类内距离比值α_j

附图4中以两种故障为例，介绍了本发明基于核密度估计和K-L散度的分类器原理，应用于此实例，具体如以下步骤4到6所述。

步骤4：利用核密度估计计算训练样本中不同故障类别样本集关于步骤3中提取的敏感特征的概率密度函数和加入一个未知故障类别待测样本后各类样本集新的概率密度函数。

针对数据集A，各类故障样本集训练样本数是35，带入（8）式有：

$h={\left(\frac{4{\widehat{\mathrm{\σ}}}^5}{3\×35}\right)}^{\frac{1}{5}}\≈{1.06\widehat{\mathrm{\σ}}\×35}^{-\frac{1}{5}}$

则可以计算出核密度估计的最佳带宽h。

基于上述方法，可得到四类样本集的关于第j特征的原始核密度函数以及加入一个相同待测样本后四类样本集新的核密度函数示例如附图5所示。

步骤5：求出在选定的特征描述下，训练样本中各类故障样本集原始概率密度函数，以及加入一个待测样本后的新概率密度函数两者的K-L散度值。

由步骤4得到数据集A中四类样本集的原始密度函数和加入测试样本后的新密度函数之后，根据对称化K-L散度定义，有：

$D_{\mathrm{KL}}(P,Q)=\frac{1}{2}(D_{\mathrm{KL}}\left(P\right|\left|Q\right)+D_{\mathrm{KL}}\left(Q\right|\left|P\right))---\left(10\right)$

求得在步骤3中选择的10个敏感特征描述下，数据集A中各类样本集原始密度函数与新密度函数的K-L散度值（j=1,2,…,10;i=1,2,…,4）。

步骤6：计算集成K-L散度，并通过集成K-L散度的大小判断待测样本的故障类别。

由步骤5可知，我们计算得到了前10个特征描述下的K-L散度值，根据式（11）中集成K-L散度值IKL_i的定义，有：

${\mathrm{IKL}}_i=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{F}_j\×{\mathrm{KL}}_i^j,i=\mathrm{1,2},\·\·\·,4$

式中，F_j为步骤3得到的前10个特征分类贡献率。接着根据四类故障样本集IKL_i值的大小，对待测样本的类别做出判定。

根据以上步骤，可以求得所有数据集的分类结果，如表6所示，这里以基于SVM和BP神经网络的故障诊断方法作为参考。

表6三种方法的分类准确率比较

除此之外，我们进一步考虑了选择的敏感特征个数和训练样本数对本方法分类准确率的影响，结果如附图6和附图7所示。综上结果可以看出，本发明分类准确率非常理想，方法的推广性和稳健性也比较好。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.一种基于核密度估计和K-L散度的旋转机械故障诊断方法，包括如下步骤：

步骤1：采集被监测对象的原始振动数据，并划定训练样本集和测试样本集；

步骤2：从步骤1中得到的原始振动数据中提取指定时频域特征；

步骤3：从步骤2中得到的频域特征集中选择出敏感特征，并且计算这些敏感特征的分类贡献率；

步骤4：利用核密度估计计算训练样本中不同故障类别样本集关于步骤3中提取的敏感特征的概率密度函数，并计算加入一个未知故障类别待测样本后各类样本集新的概率密度函数；

步骤5：计算出在选定的特征描述下，训练样本中各类故障样本集原始概率密度函数，以及加入一个待测样本后的新概率密度函数两者的K-L散度值；

步骤6：计算集成K-L散度，并通过集成K-L散度的大小判断待测样本的故障类别；

所述步骤2中得到的指定时频域特征是通过总体平均经验模态分解方法和希尔伯特变换两种信号处理方法得到的；

所述步骤3的计算过程如下：

第一步：计算第j个特征C个类的类内距离的平均值

$\left\{\begin{array}{c}{d}_{c,j}=\frac{1}{M_c\×(M_c-1)}\underset{l,m=1}{\overset{M_c}{\Σ}}|q_{m,c,j}-q_{l,c,j}|\\ l,m=1,2,...,M_c,l\≠m,j=1,2,...,J,c=1,2,...,C\end{array}\right.---\left(1\right)$

$d_j^{\left(w\right)}=\frac{1}{C}\underset{c=1}{\overset{C}{\Σ}}{d}_{c,j}---\left(2\right)$

其中，M_c表示第c类的样本个数，J表示特征个数，C表示类别个数，q_m,c,j表示第c类第m个样本的第j个特征的特征值；

第二步：计算第j个特征C个类的类间距离的平均值

$u_{c,j}=\frac{1}{M_c}\underset{m=1}{\overset{M_c}{\Σ}}{q}_{m,c,j}---\left(3\right)$

$\left\{\begin{array}{c}{d}_j^{\left(b\right)}=\frac{1}{C\×(C-1)}\underset{c,e=1}{\overset{C}{\Σ}}|u_{e,j}-u_{c,j}|\\ c,e=1,2,...,C,c\≠e\end{array}\right.---\left(4\right)$

其中，u_c,j、u_e,j分别表示第c和第e个类第j个特征的平均值；

第三步：数据集A的类间距离与类内距离的比值α_j：

${\mathrm{\α}}_j=\frac{d_j^{\left(b\right)}}{d_j^{\left(w\right)}}---\left(5\right)$

第四步：计算选择的前n个特征的分类贡献率F_j；定义：

$F_j=\frac{{\mathrm{\α}}_j}{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\α}}_i},(j=1,2,...,n)---\left(6\right)$

式(6)为分类贡献率。

2.如权利要求1所述的基于核密度估计和K-L散度的旋转机械故障诊断方法，其特征在于：所述核密度估计是概率论中一种用来估计随机变量概率密度函数的非参数方法，设X₁,X₂,…,X_n是取自一元连续总体X的样本，在任意点x处的总体密度函数f(x)的核密度估计定义为：

${\hat{f}}_h\left(x\right)=\frac{1}{nh}\underset{i=1}{\overset{n}{\Σ}}K\left(\frac{x-x_i}{h}\right)---\left(7\right)$

其中，K(·)称为核函数，h为窗宽。

3.如权利要求2所述的基于核密度估计和K-L散度的旋转机械故障诊断方法，其特征在于：所述的窗宽h是核密度估计方法中唯一要优化的参数，根据最小化平均集成平方误差方法，得到下式：

$h={\left(\frac{4{\widehat{\mathrm{\σ}}}^5}{3n}\right)}^{\frac{1}{5}}\≈1.06\widehat{\mathrm{\σ}}{n}^{-\frac{1}{5}}---\left(8\right)$

其中，是样本标准差，n为各类样本集所包含的样本数。

4.如权利要求1所述的基于核密度估计和K-L散度的旋转机械故障诊断方法，其特征在于：所述步骤5的运算过程如下：

原始的K-L散度定义式为：

$D_{KL}\left(P\right|\left|Q\right)=\underset{i}{\Σ}h\left(\frac{P\left(i\right)}{Q\left(i\right)}\right)P\left(i\right)---\left(9\right)$

为了满足距离的对称性，重新定义为：

$D_{KL}(P,Q)=\frac{1}{2}\left(D_{KL}\right(P\left|\right|Q)+D_{KL}(Q\left|\right|P\left)\right)---\left(10\right)$

由步骤4得到各类样本集的原始密度函数和加入测试样本后的新密度函数之后，就可以根据式(10)计算出对应的K-L散度值。

5.如权利要求1所述的基于核密度估计和K-L散度的旋转机械故障诊断方法，其特征在于：所述步骤6的计算过程如下：

定义了一个集成K-L散度值IKL_i(i＝1,2,…,C)：

${\mathrm{IKL}}_i=\underset{j=1}{\overset{n}{\Σ}}{F}_j\×{\mathrm{KL}}_i^j---\left(11\right)$

其中，F_j为步骤3得到的前n个特征分类贡献率，式(11)中通过将第j个特征描述下的K-L散度值与第j个特征的分类贡献率F_j加权，可以得到各类样本集的集成K-L散度值。