Optimization and Research of a Bearing Fault Diagnosis Method Based on VMD and KFCM

VMD是一种非递归的、准正交的、自适应的信号分解算法，采用交替方向乘子法(Altermate Direction Method of Multipliers，ADMM)，其实质是变分问题的构造与求解，使得每个模态的估计带宽之和最小，其中假设每个模态是具有不同中心频率的有限带宽，为解决这一变分问题，采用了交替方向乘子法，不断更新各模态及其中心频率，逐步将各模态解调到相应的基频带，最终各个模态及相应的中心频率被一同提取出来^[15].

在VMD分解中，为使各个模态之间带宽最小，变分模态u_k(t)之和等于原始信号f(t).对于每个模态u_k(t)，用希尔伯特变换估计分析每个模态的单边频谱，然后，与预估中心频率$ \mathrm{e}^{-j \omega_{k} t}$相乘，将各个模态的单边谱频谱对应到基频带上，最后，利用梯度的平方L²范数，计算每个模态信号的带宽，得到变分约束模型：

其中，ω_k(t)为各模态对应的中心频率，δ(t)为Dirac分布，k为模态分量数.

为了将约束变分问题转换为非约束变分问题，在求解约束变分模型时，引入增广拉格朗日函数：

其中，α为二次惩罚因子，保证输入信号f(t)在高斯噪声的影响下信号的重构精度; λ(t)为拉格朗日乘法算子.

采用ADMM计算出增广拉格朗日函数的最优解，即可通过VMD将原始信号f(t)分解成K个窄带本征模态函数分量^[16-17]. VMD算法步骤：

1) 初始化$\left\{ \overset{\wedge }{\mathop{u}}\,_{k}^{1} \right\},\left\{ \omega _{k}^{1} \right\},\left\{ {{\overset{\wedge }{\mathop{\lambda }}\,}^{1}} \right\}$和n，n=0;

2) n=n+1，开始整个循环;

3) 更新每个模态的频谱;

4) 更新中心频率;

5) 更新拉格朗日乘子

其中，τ为拉格朗日乘子更新参数.

6) 重复步骤2~5，对于给定判别精度e＞0^[18]，直到满足迭代条件：

得到k个窄带本征模态函数分量，迭代结束.若不满足迭代条件，返回步骤2，最终分解结果如图 1所示.

KFCM是FCM的内核版本，该算法通过核函数将原始数据点从输入空间映射到高维特征空间，如图 2所示.考虑到原始数据中的点无法用一个线性函数进行划分，于是将其变换到一个更高维度的空间中，可以在这个高维空间中找到一个线性函数，容易对原始数据进行划分，这个高维空间就叫特征空间.从低维到高维空间的映射函数的内积就叫核函数，将核函数引入机器学习的一个重要原因是：当特征空间维数很高而核函数计算量较之特征空间内的内积运算计算量很小时，这样可以提高计算效率.

核函数的定义：设X∈R^s，定义从X到特征空间H的映射：Φ：X→H：Φ(x)=y，则

其中，x和$ {\tilde x}$为s维向量，$(x, \tilde x) $为两者的欧式内积.

其中‖Φ(x_k)-Φ(v_i)‖_H²是Φ(x_k)到Φ(v_i)之间距离的平方，K(x_k，v_i)是高斯径向函数，其形式如下：

利用拉格朗日乘子法，可求得：

KFCM算法的步骤可归纳如下：

步骤1  设定径向基函数的参数σ; 聚类个数c，模糊指数m，收敛精度ε; 令迭代次数k=0;用FCM算法初始化中心矩阵V⁽⁰⁾;

步骤2   用式(14)计算U^(k+1);

步骤3  用式(15)计算V^(k+1)，令k=k+1;

重复步骤(2)和(3)，直到满足如下的终止条件：

$\left\|U^{(k)}-U^{(k-1)}\right\| <\varepsilon $或存在i(1≤i≤c)使得$\sum\limits_{j=1}^{N} u_{i j}=0 $.

为了便于分析KFCM分类的准确程度，本文采用模糊数学的贴近度来判断，即计算测试样本T和第i个状态C_i的海明贴近度^[19]：

根据计算分析，贴近度最大者即分为一类，如果2组或2组以上测试样本得到的贴近度都较大，且相差小于0.005时，认为无法辨识或被错误分类.

对于KFCM的分类性能，本文采用分类系数F、平均模糊熵H以及标准FCM分类成功率v进行评价^[16]：

式中：w为运行10次标准FCM，成功得到聚类中心的次数.

VMD-KFCM故障诊断流程如图 3所示.故障诊断基本步骤如下：

1) 将故障振动信号经过VMD分解得到若干模态分量u_k(t).

2) 通过分析各个模态中心频率、皮尔逊相关系数(Pearson correlation coefficient，PCC)，确定模态分量.

3) 对选定分量重构提取时域特征，提取各个u_k(t)分量奇异值矩阵作为特征量，形成特征量参数矩阵.

4) 对特征量参数矩阵归一化，作为KFCM的输入，进行不同故障程度样本的模型训练，根据测试样本对模型进行测试和验证.

本文采用凯斯西储大学轴承数据中心网站的滚动轴承振动数据进行分析，如图 4所示，实验数据选用6205-2RS SKF轴承，人为设置故障点直径为0.007 mm，选用负荷0(179 7 r/min)、采样频率12 kHz的故障数据.以MATLAB为仿真平台，与EMD特征提取方法进行对比分析.

由于各个模态之间的主要区别是具有不同中心频率，本文通过计算每个分量的功率谱密度、分析中心频率，并计算PCC确定K值. PCC是衡量向量相似度的一种方法.输出范围为-1到1，0代表无相关性，负值代表负相关，正值代表正相关.假设样本可以记为(X_i，Y_i)，则样本的PCC为

其中$ \frac{X_{i}-\bar{X}}{s_{X}}, \bar{X}, s_{X}$分别为标准化变量，样本均值和样本标准差.

具有低PCC值的模态被识别为纯噪声和趋势，其余的模态被标识为有用信号.对轴承内圈故障振动信号进行VMD分解，信号采样点数设置为4 096个，VMD的各个模态分量的PCC如图 5所示，其PCC阈值设置为0.02，当K=7时，其模态分量的PCC值仅为0.01，明显低于阈值.不同K值下的功率谱密度如图 6所示.可以看出，前6个模态功率谱密度逐渐增大，当K=7时，模态7的功率谱密度相对模态6衰减明显，为了减少分解模态过多导致信息丢失，选择最佳分解模态数K=6，以保证实际信号分解的保真度.

为表明各模态分量中信息的有效性，图 7分别给出了VMD和EMD方法各模态分量的信号包络谱，其中包含了内圈故障频率161.1 Hz和电机2倍转频58.59 Hz，但EMD模态4、模态5、模态6的包络谱中没有明显的故障特征，因此EMD采用前3个互信息较大的分量进行特征提取.采用同样的方法，对其故障信号进行分解，最终确定VMD中K=6，和前面对K值分析是一致的.

选取滚动轴承的正常、内圈、外圈和滚动体故障4种状态下各40组数据，其中，每组数据截取2 048个采样点，前20组数据作为已知训练样本，分别采用基于VMD和EMD的方法求取已知故障的标准聚类中心，后20组数据作为测试样本，通过贴近度计算进行故障识别.

结果如表 1所示，KFCM方法计算速度最快，但从分类系数、故障识别率和平均模糊商3个指标不难看出该方法分类效果很差; EMD-KFCM方法虽然整体性能参数都有所提升，但相对于VMD-KFCM-PSO方法性能较差一些; VMD-KFCM-PSO方法的分类系数为0.99，且该方法运行1次即可得到正确的分类中心且分类系数较EMD高11个百分点，相比较VMD的平均模糊熵也更小，运算时间更短，对于故障在线分析有明显优势.

通过对20组测试样本的贴近度求平均，得到4类测试样本相对标准聚类中心的平均贴近度.图 8为VMD-KFCM方法的测试样本平均贴近度分布图，可以看出，基于VMD的测试样本与标准聚类中心的最大贴近度都在0.95以上，同时也明显高于其他3个贴近度.因此，本方法能够更精确地提取出各工况特征，更容易实现故障类型划分.

本文进一步分析了4种故障类型的测试样本聚类中心样本特征曲线，该曲线是将4种类型的奇异值相连.通过与原聚类中心点的变化趋势对比，如图 9所示，4类测试样本特征曲线都在标准聚类中心曲线附近，没有发生明显偏离，从而保证实验测试是成功的，分类识别率保持99.12%.

图 10为粒子群优化算法对KFCM用于训练状态的参数进行优化处理，在迭代10次以后趋于稳定状态; 隶属度矩阵值如图 11所示，依次是正常、内圈故障、外圈故障、滚动体故障数据的隶属度矩阵值，图 12为聚类中心及分布图，很明显每组划分都在聚类中心附近.

本文以电机轴承故障为例，分析验证改进的机械系统故障诊断新方法，该方法充分结合模糊聚类算法与变分模态分解法来解决数据受到噪声和异常值影响时的缺陷，从而提高分类准确率.

1) 与EMD-KFCM方法相比较，VMD-KFCM-PSO方法能够有效克服EMD分解中的模态混叠效应和虚假分量问题.

2) KFCM算法通过核函数形成一种映射关系，将原始空间中的数据点转换到高维特征空间进行计算与分析，最后得到原始空间的最优划分.从而提高了聚类性能，使算法对噪声和孤立点具有较好的鲁棒性.

3) VMD特征提取方法中K值需要通过分析中心频率和计算皮尔逊相关系数来确定，目前尚无其他理论依据，有待进一步优化和完善.