CN111175046A - 一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于流形学习和s‑k‑means聚类的滚动轴承故障诊断方法，属于旋转机械故障诊断领域。该方法主要包括5个步骤：获取轴承在正常及故障状态下的振动信号作为原始信号；获取时域特征集和小波包能量特征集；利用散布矩阵，构造特征分类能力指数，排除部分无关特征，实现特征选取；将经过特征选择后的特征集作为维数约简算法的高维输入，使用NPE(邻域保持嵌入)的流形学习方法进行维数约简，得到经过约简后的低维特征集；采用s‑k‑means聚类方法对故障特征矩阵进行聚类分析，确定故障类型。该方法简单有效，有助于实现高维故障特征的有效约简和二次提取，得到了与故障联系密切的少量典型特征，使得特征分类时的计算量大大减少，从而准确识别故障类型。

Description

一种基于流形学习和s-k-means聚类的滚动轴承故障诊断 方法

技术领域

本发明属于旋转机械故障诊断领域，具体涉及一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法。

背景技术

滚动轴承作为旋转机械的关键零部件，它在旋转机械中起着关键的作用，其工作状态的好坏将直接影响到整台机械设备的整体性能。滚动轴承故障是导致旋转机械设备发生故障的主要原因之一，严重时甚至可能导致重大财产损失。因此，为了避免由轴承导致的机械故障，减少经济损失，对轴承进行状态监测及故障诊断从而保证其正常运行非常有必要。当滚动轴承出现局部缺陷故障时，会产生一系列的冲击信号，表现为振动信号的突变。这些变化含有丰富的故障信息特征，因而振动信号的检测在机械故障诊断领域得到广泛采用。

振动信号分析主要包括三个步骤：数据采集与预处理，特征提取，模式识别。其中，特征选择和提取是故障诊断的关键步骤，模式识别是核心。对于不同的故障只有选择合适的特征和有效的模式识别方法，才能保证故障诊断结构的可靠性。

不同的故障可以根据不同特征来识别，但是不同特征对故障的敏感程度也不相同。有些特征与故障密切相关，另一些对于故障诊断来说可能是冗余甚至无关的特征。在将特征输入分类器之前，排除对故障敏感度低的特征，保留对故障敏感度高的特征，不但可以提高分类的准确性而且有利于减少计算量，避免维数灾难。

邻域保持嵌入(NPE)是一种典型的非线性维数约简方法，在尽可能地保证数据间的几何关系和距离测度不变的前提下，获得原始观测空间的真实流形结构，保留非线性流形结构的内在低维特征，挖掘高维数据中的本质信息，适用于滚动轴承早期故障特征信息的挖掘和提取。

聚类问题是将一个数据点集合中的元素，按某种相似程度的度量，分别赋予不同的类别标号，从而实现数据区分。K-means聚类是一种无监督动态聚类方法，具有一定的适应性，可以有效进行聚类分析。但是K-means聚类算法严重依赖于初始聚类中心的选取，易受孤立点的影响，且需要事先设定聚类数。s-k-means聚类算法可以根据数据集自动给出最佳的初始聚类中心和聚类数，实现自动分类的效果，从而准确实现故障类型划分。

发明内容

针对故障轴承振动信号的非平稳性及高维特征选择问题，本发明公开了一种基于流形学习和s-k-means聚类分析的滚动轴承早期故障诊断方法。

本发明包括以下步骤：

步骤一、获取滚动轴承在不同运行状态下的振动加速度信号，得到时域信号样本集；

步骤二、计算样本集中全部样本的时域特征和小波包能量特征，构成初始高维特征集；

具体进行如下处理：

步骤2.1：分别计算每个样本的方根幅值、偏度、方差、绝对平均值、峭度、最大值、波形指标、脉冲指标、峭度指标、峰值指标、裕度指标、偏度指标，构成时域特征集。

步骤2.2：对每个样本信号进行三层小波包分解，计算信号每一个分解层次上每个小波包空间中信号分布的能量，构成了小波包频带能量特征集。

步骤三、利用散布矩阵，构造特征分类能力指数，筛选特征，实现特征选取；

具体进行如下处理：

步骤3.1：计算原始高维特征集的类内散布矩阵S_w、类间散布矩阵S_b和混合散布矩阵S_m。

设原始高维特征集为X＝{x₁,x₂,…,x_D}，其中D表示高维特征集的维数，x_i表示第i个特征。则类内协方差矩阵为：

S_j表示第j类数据协方差矩阵，x_j表示第j类数据特征值，μ_j表示第j类特征值均值。则类内散布矩阵S_w为：

其中

表示全部样本中属于第j类的样本数。

类间散布矩阵S_b为：

其中μ是全部样本的全局均值向量，

混合散布矩阵：S_m＝E[(x-μ₀)(x-μ₀)^T]，即S_m是全局均值向量的协方差矩阵。且S_m＝S_b+S_w

步骤3.2：构造分类能力指数：

通过分类能力指数J排除部分对进行故障诊断敏感度低的特征。

步骤四、将经过特征选择后的特征集作为维数约简算法的高维输入，使用邻域保持嵌入(NPE)的流形学习方法进行维数约简，得到经过约简的低维特征集；

具体进行如下处理：

步骤4.1：构建邻域图

利用K近邻(KNN)算法来构建邻域图。假设共有m个样本，则邻域图共有m个结点，其中X_i表示第i个结点。如果X_j是X_i的k个最近邻居之一，那么将这两点连接，反之不连接。

步骤4.2：计算权重矩阵

设权重矩阵为W，其中元素W_ij代表结点i和结点j之间的边的权重，如果两点之间没有连接，则对应的矩阵元素为0。矩阵W的元素值主要通过最小化如下目标函数得到：

其中W应满足归一化约束：

步骤4.3：计算映射

通过求解广义特征向量问题来计算降维的线性映射：XMX^Ta＝λXX^Ta。其中数据集X＝(x₁,…x_m)。矩阵M＝(I-W)^T(I-W)，矩阵I＝diag(1,…,1)。

按特征值从小到大的次序(λ₀≤…≤λ_d-1)将求解到的特征向量进行排列a₀,…,a_d-1，这样经过降维后的嵌入坐标则为：y_i＝A^Tx_i，其中A＝(a₀,a₁,…,a_d-1)。

步骤五、采用s-k-means聚类方法对故障特征矩阵进行聚类分析，确定故障类型。

设D维空间的n个数据点的集合为X＝{x₁,x₂,...,x_n}，对其应用新的聚类分析划分为j个类簇C_j。

具体进行如下处理：

步骤5.1：计算每一个数据点的密度指标

式中:P_i为点x_i的密度指标，若某个元素的邻近点越多，则其密度值越大。a为x_i的邻域半径，在半径以外的点对点x_i的密度指标贡献很小。选择其中密度值最大的数据对象x_c,i，则可得到的首个类中心点。

步骤5.2：设x_c，k为第k次选出的聚类中心，P_c，k为其对应的密度指标，则对于每个数据点的密度指标按照式以下公式修正

式中:b表示密度指标函数显著减少的邻域。选取密度指标最高的数据点x_c，k+1，作为新的聚类中心。

步骤5.3：判断

式中:δ是1个给定的参数。如果不成立则转步骤5.2，如果成立则退出。

步骤5.4：设I＝1，将求取的前k个数据点x_c，j，j＝1,2,…,k，记为m_j(I),j＝1,2,…,k。

步骤5.5：计算每一个数据点与这k个聚类中心的距离d(x_i,m_j(I)),i＝1,2,...n,j＝1,2,..,k。如果满足d(x_i,m(I))＝min{d(x_i,m_j(I)),j＝1,2,…，k}则x_i∈C_j。

步骤5.6：计算k个新的聚类中心

式中:j＝1,2,…,k,x_i∈C_j，N_j为属于第j个聚类中心C_j的数据集x_i的个数。

步骤5.7：判断，如果m_j(I+1)≠m_j(I),j＝1,2…,k。则I＝I+1，并返回步骤5.5；否则聚类过程结束，输出聚类结果。

本发明的优点与积极效果在于：

(1)计算量小，可有效避免维数灾难

首先通过特征分类能力指数对原始特征集进行筛选，可排除部分对故障不敏感的特征，而通过邻域保持嵌入(NPE)的流形学习降维方式对经过筛选的特征集进行学习提取数据内在的信息结构，将复杂的高维特征空间转化为简单的低维特征空间，极大程度地保留故障信号特征集合中内含的整体几何结构信息。实现高维特征的有效约简和二次提取，得到最佳滚动轴承故障特征。

(2)模式识别方便有效

s-k-means聚类方法首先应用最高密度指标寻找聚类中心，将寻找到的聚类中心作为初始聚类中心；然后再利用K-means聚类进行聚类分析。克服了传统K-means聚类对初始聚类中心敏感且容易受孤立点的影响的缺陷，且可以根据数据集自动给出最佳的聚类数，准确实现故障类型划分。

附图说明

图1是本发明的滚动轴承故障诊断方法整体步骤流程图；

图2(a)是本发明实例中不同故障聚类分类结果图；

图2(b)是本发明实例中不同故障聚类轮廓图；

图3(a)是本发明实例中内圈不同尺寸故障聚类分类结果图；

图3(b)是本发明实例中内圈不同尺寸故障聚类轮廓图。

具体实施方式

本发明是一种滚动轴承故障诊断算法，下面将结合附图，对本发明作进一步的详细说明。

以凯斯西储大学轴承振动数据库的部分实验数据为例，采用的轴承为SKF 6205单列深沟球轴承。分别对正常轴承，0.007英寸的单点内圈故障，0.007英寸的单点外圈故障，0.007英寸的单点滚动体故障和4种不同程度(0.007英寸、0.014英寸、0.021英寸、0.028英寸)的单点内圈故障(1英寸＝2.54厘米)，进行分析。

本发明公开的滚动轴承故障诊断方法的整体步骤流程图如图1所示，具体步骤如下：

1.获取滚动轴承在不同运行状态下的振动加速度信号，得到时域信号样本集；

2.计算样本集中每个样本的时域特征和小波包能量特征，构成初始高维特征集；

具体进行如下处理：

(1)：本发明以12个典型的时域统计特征参数为例，构成初始时域特征集。这12个统计参数分别为方根幅值、偏度、方差、绝对平均值、峭度、最大值、波形指标、脉冲指标、峭度指标、峰值指标、裕度指标、偏度指标。

(2)：设原始信号xt)，对信号x(t)进行三层小波包分解，计算信号x(t)每一个分解层次上每个小波包空间中信号分布的能量，构成了x(t)的14个小波包频带能量特征集。总共可以得到26个原始高维特征向量X＝{x₁,x₂,...,x₂₆}，并对原始特征向量进行归一化处理。

3.利用散布矩阵选择特征

具体进行如下处理：

(1)：计算原始高维特征集的类内散布矩阵S_w、类间散布矩阵S_b和混合散布矩阵S_m。

设原始高维特征集为X＝{x₁,x₂,…,x_D}，其中D表示原始高维特征集的维数26，x_i表示第i个特征。则类内协方差矩阵为：

S_j表示第j类数据协方差矩阵，x_j表示第j类数据特征值，μ_j表示第j类特征值均值。则类内散布矩阵S_w为：

其中

表示全部样本中属于第j类的样本数。

类间散布矩阵S_b为：

其中μ是1样本的全局均值向量，

混合散布矩阵：S_m＝E[(x-μ₀)(x-μ₀)^T]，即S_m是全局均值向量的协方差矩阵，且S_m＝S_b+S_w。

(2)：构造分类能力指数：

通过分类能力指数J排除16个对进行故障诊断敏感度低的特征，实现特征选取。

4.通过极大似然法确定不同情况下流形学习约简后的维数，将经过特征选择后的特征集作为维数约简算法的高维输入，使用NPE(邻域保持嵌入)的流形学习方法进行维数约简，得到约简后的低维特征集；

具体进行如下处理：

(1)：构建邻域图

利用K近邻(KNN)算法来构建邻域图。假设共有m个样本，则邻域图共有m个结点，其中X_i表示第i个结点。如果X_j是X_i的k个最近邻居之一，那么将这两点连接，反之不连接。

(2)：计算权重矩阵

设权重矩阵为W，其中元素W_ij代表结点i和结点j之间的边的权重，如果两点之间没有连接，则对应的矩阵元素为0。矩阵W的元素值主要通过最小化如下目标函数得到：

其中W应满足归一化约束：

(3)：计算映射

通过求解广义特征向量问题计算降维的线性映射：XMX^Ta＝λXX^Ta。其中X＝(x₁,…x_m)。矩阵M＝(I-W)^T(I-W)，矩阵I＝diag(1,…,1)。

按特征值从小到大的次序(λ₀≤…≤λ_d-1)将求解到的特征向量进行排列a₀,…,a_d-1，这样降维后的嵌入坐标则为：y_i＝A^Tx_i，其中A＝(a₀,a₁,…,a_d-1)。

5.采用s-k-means聚类方法对经过流形学习得到的低维故障特征矩阵进行聚类分析，确定故障类型。

具体进行如下处理：

设D维空间的n个数据点的集合为X＝{x₁,x₂,...,x_n}，对其应用新的聚类分析划分为j个类簇C_j具体算法流程为：

(1)：计算每一个数据点的密度指标

式中:P_i为点x_i的密度指标，若某个元素的邻近点越多，则其密度值越大。a为x_i的邻域半径，在半径以外的点对点x_i的密度指标贡献很小。选择其中密度值最大的数据对象x_c,i，则可得到的首个类中心点。

(2)：设x_c，k为第k次选出的聚类中心，P_c，k为其对应的密度指标，则对于每个数据点的密度指标按照式以下公式修正

式中:b表示密度指标函数显著减少的邻域。选取密度指标最高的数据点x_c，k+1，作为新的聚类中心。

(3)：判断

式中:δ是1个给定的参数。如果不成立则转步骤(2)，如果成立则退出。

(4)：设I＝1，将求取的前k个数据点x_c，j，j＝1,2,…,k，记为m_j(I),j＝1,2,…,k。

(5)：计算每一个数据点与这k个聚类中心的距d(x_i,m_j(I)),i＝1,2,...n,j＝1,2,..,k。如果满足d(x_i,m(I))＝min{d(x_i,m_j(I)),j＝1,2,…，k}则x_i∈C_j。

(6)：计算k个新的聚类中心

式中:j＝1,2,…,k,x_i∈C_j，N_j为属于第j个聚类中心C_j的数据集x_i的个数。

(7)：判断，如果m_j(I+1)≠m_j(I),j＝1,2…,k。则I＝I+1，并返回步骤(5)；否则聚类过程结束，输出聚类结果。

对得到的聚类结果进行分析，发现对四种不同运行状态(正常轴承、内圈故障、外圈故障、滚动体故障)和四种不同尺寸(0.007英寸、0.014英寸、0.021英寸、0.028英寸)的内圈故障轴承的故障辨识率分别为100％和99％，且故障辨识结果稳定。

通过以上滚动轴承故障诊断方法的详细描述，可见本发明的基于流形学习和

s-k-means聚类的滚动轴承故障诊断方法具有明显的优势：

(1)计算量小，可有效避免维数灾难

经过特征选择和邻域保持嵌入的流形学习降维方式得到了与故障联系密切的少量典型特征，使得特征分类时的计算量大大减少，也可有效避免特征数过高而导致的维数灾难。

(2)模式识别方便有效

s-k-means聚类分析通过减法聚类采用密度指标最高原则优化初始聚类中心，克服了传统K-means聚类严重依赖于初始聚类中心的选取的缺陷。并且可以自动确定有效的聚类类别数，实现自动分类的效果，能给出较为合理的聚类划分结果。

Claims (5)

1.一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法，其特征在于：该方法包括以下步骤：

步骤1、获取滚动轴承在不同运行状态下的振动加速度信号，得到时域信号样本集；

步骤2、计算样本集中全部样本的时域特征和小波包能量特征，构成初始高维特征集；

步骤3、利用散布矩阵，构造特征分类能力指数，筛选特征，选取最优特征子集，实现特征选取；

步骤4、将经过特征选择后的特征集作为维数约简算法的高维输入，使用邻域保持嵌入NPE的流形学习方法进行维数约简，得到经过约简的低维特征集；

步骤5、采用s-k-means聚类方法对故障特征矩阵进行聚类分析，确定故障类型；

设D维空间的n个数据点的集合为X＝{x₁,x₂,...,x_n}，对其应用新的聚类分析划分为j个类簇C_j。

2.根据权利要求1所述的一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法，其特征在于：步骤2的具体实施过程如下，

步骤2.1：分别计算每个样本的方根幅值、偏度、方差、绝对平均值、峭度、最大值、波形指标、脉冲指标、峭度指标、裕度指标、偏度指标，构成时域特征集；

步骤2.2：对每个样本信号进行三层小波包分解，计算信号每一个分解层次上每个小波包空间中信号分布的能量，构成了小波包频带能量特征集。

3.根据权利要求1所述的一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法，其特征在于：步骤3的具体实施过程如下，

步骤3.1：计算原始高维特征集的类内散布矩阵S_w、类间散布矩阵S_b和混合散布矩阵S_m；

设原始高维特征集为X＝{x₁,x₂,…,x_D}，其中D表示高维特征集的维数，x_i表示第i个特征；则类内协方差矩阵为：

S_j表示第j类数据协方差矩阵，x_j表示第j类数据特征值，μ_j表示第j类特征值均值；则类内散布矩阵S_w为：

其中

表示全部样本中属于第j类的样本数；

类间散布矩阵S_b为：

其中μ是全部样本的全局均值向量，

混合散布矩阵：S_m＝E[(x-μ₀)(x-μ₀)^T]，即S_m是全局均值向量的协方差矩阵；且S_m＝S_b+S_w

步骤3.2：构造分类能力指数：

通过分类能力指数J排除部分对进行故障诊断敏感度低的特征。

4.根据权利要求1所述的一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法，其特征在于：步骤4的具体实施过程如下，

步骤4.1：构建邻域图

利用K近邻KNN算法来构建邻域图；假设共有m个样本，则邻域图共有m个结点，其中X_i表示第i个结点；如果X_j是X_i的k个最近邻居之一，那么将这两点连接，反之不连接；

步骤4.2：计算权重矩阵

设权重矩阵为W，其中元素W_ij代表结点i和结点j之间的边的权重，如果两点之间没有连接，则对应的矩阵元素为0；矩阵W的元素值通过最小化如下目标函数得到：

其中W应满足归一化约束：

步骤4.3：计算映射

通过求解广义特征向量问题来计算降维的线性映射：XMX^Ta＝λXX^Ta；数据集X＝(x₁,…x_m)；矩阵M＝(I-W)^T(I-W)，矩阵I＝diag(1,…,1)；

按特征值从小到大的次序(λ₀≤…≤λ_d-1)将求解到的特征向量进行排列a₀,…,a_d-1，这样经过降维后的嵌入坐标则为：y_i＝A^Tx_i，其中A＝(a₀,a₁,…,a_d-1)；

5.根据权利要求1所述的一种基于流形学习和s-k-means聚类的滚动轴承故障诊断方法，其特征在于：步骤5的具体实施过程如下，

步骤5.1：计算每一个数据点的密度指标

式中:P_i为点x_i的密度指标，若某个元素的邻近点越多，则其密度值越大；a为x_i的邻域半径，在半径以外的点对点x_i的密度指标贡献很小；选择其中密度值最大的数据对象x_c,i，则可得到的首个类中心点；

步骤5.2：设x_c，k为第k次选出的聚类中心，P_c，k为其对应的密度指标，则对于每个数据点的密度指标按照式以下公式修正

式中:b表示密度指标函数显著减少的邻域；选取密度指标最高的数据点x_c，k+1，作为新的聚类中心；

步骤5.3：判断

式中:δ是1个给定的参数；如果不成立则转步骤5.2，如果成立则退出；

步骤5.4：设I＝1，将求取的前k个数据点x_c，j，j＝1,2,…,k，记为m_j(I),j＝1,2,…,k；

步骤5.5：计算每一个数据点与这k个聚类中心的距d(x_i,m_j(I)),i＝1,2,...n,j＝1,2,..,k；如果满足d(x_i,m(I))＝min{d(x_i,m_j(I)),j＝1,2,…，k}则x_i∈C_j；

步骤5.6：计算k个新的聚类中心

式中:j＝1,2,…,k,x_i∈C_j，N_j为属于第j个聚类中心C_j的数据集x_i的个数；

步骤5.7：判断，如果m_j(I+1)≠m_j(I),j＝1,2…,k；则I＝I+1，并返回步骤5.5；否则聚类过程结束，输出聚类结果。

Images