基于参数优化VMD-KPCA和BP网络的齿轮故障诊断方法

Fault diagnosis method of gear based on parameter optimization VMD-KPCA and BP network

摘要

针对噪声环境下难以提取齿轮故障特征、诊断准确率低的问题，提出了一种将基于综合评价指标的变分模态分解（variational mode decomposition，VMD）参数寻优、核主成分分析（kernel principal component analysis，KPCA）特征融合和BP网络相结合的齿轮故障诊断方法。首先，为了有效评价VMD分解后的各固有模态函数（intrinsic mode function，IMF）分量及避免变分模态分解需要人为设定相关参数的问题，设计了一种基于包络熵与峭度的综合评价指标，用于建立VMD参数寻优中的适应度函数及筛选最优IMF分量；其次，按最优参数进行变分模态分解后对最优IMF分量提取多域特征集，再利用KPCA模型对其进行特征的融合；最后，通过BP网络模型进行故障诊断。结果表明，与其他传统方法相比，在相同实验条件下，该方法提高了齿轮故障的识别率，准确率高达98.00%，证明了该方法的有效性。

Abstract

In view of the problems of difficulty in extracting gear fault features and low diagnosis accuracy under noisy environment，a gear fault diagnosis method was proposed，which combined VMD parameter optimization based on comprehensive evaluation indicators，KPCA feature fusion and BP network.Firstly，in order to effectively evaluate the IMF components after VMD decomposition and avoid the problem of manually setting relevant parameters for VMD，a comprehensive evaluation index based on envelope entropy and kurtosis was designed to establish a fitness function for VMD parameter optimization and screen the optimal IMF components.Secondly，after performing VMD decomposition according to the optimal parameters，a multi domain feature set was extracted from the optimal IMF component，and then the KPCA model was used to fuse its features.Finally，fault diagnosis was performed using the BP network model.The experiment shows that under the same experimental conditions，compared with other traditional methods，this method improves the recognition rate of gear faults，with an accuracy of up to 98%，proving the effectiveness of this method.

Graphical abstract

关键词

齿轮故障诊断 / 综合评价指标 / 变分模态分解 / 鲸鱼优化算法 / BP网络

Key words

gear fault diagnosis / comprehensive evaluation index / variational mode decomposition / whale optimization algorithm / BP network

旋转机械设备广泛应用于航空航天、汽车、机床等行业，而齿轮是旋转设备中最具代表性的关键部件之一。由于齿轮的健康状态直接关系设备能否安全正常运行，更关系到人身安全^［1］，因此尽早识别出齿轮中存在的故障，对保证机械设备的稳定、安全运行具有重要意义。

大多数传统的故障诊断方法只通过时域、频域及时频域来描述振动信号的特征，主要方法有快速傅里叶变换^［2］等，但在复杂工况下采集的振动信号通常是复杂的非平稳信号，故障特征很难被描述。为了解决这个问题，一些信号分解的方法被提出，主要包括经验模态分解及其改进的多元经验模态分解^［3］、集成经验模态分解^［4］、小波变换^［5-6］等。复杂信号通过这些方法可以先分解为若干分量，再通过各种方法描述信号的特征。然而这些分解方法也伴随着端点效应和模态分量混叠的弊端^［7］。Dra-gomiretskiy等^［8］采用变分模态分解方法将多分量信号分解为带宽固定、中心频率不同的固有模态，不仅具有抗噪能力，且能够降低振动信号的复杂度和不稳定性，有效克服了其他分解方法固有的端点效应和模态分量混叠的问题。武英杰等^［9］针对经验模态分解和局部均值分解的模态混叠问题，提出将VMD应用到风电机组故障诊断中，发现对不平衡故障具有较好的诊断效果。VMD在旋转机械的故障诊断中已被广泛应用，并体现出良好的效果。

然而，VMD的效果与其参数K和\alpha的选择相关。在实际工程应用中，这些参数通常需要人为实验确定，耗时费力且很难找到最优参数。当参数选择不准确时，会出现过分解及欠分解的问题，时而影响最后的诊断效果。徐甜甜等^［10］使用经验值设置VMD参数，并结合谱峭度完成IMF选择，判断出齿轮箱轴承故障。Zhang等^［11］提出基于能量准则方法确定分解层数，但该方法对于最优分量的选择不具有自适应性。刘秀丽等^［12］提出一种VMD和小波分析相结合的振动信号处理方法，采用中心频率观察法确定变分模态分解的K值，但忽视了惩罚因子α对分解结果的影响。张萍等^［13］提出一种基于鲸鱼算法的VMD参数优化方法，并将包络熵作为评价指标，但该方法评价指标单一，不能从多角度综合评价各个IMF分量。

本文提出将包络熵与峭度倒数结合，作为综合评价指标，采用鲸鱼优化算法^［14］（whale optimization algorithm，WOA）优化VMD的参数，在不同状态下分别找到最优参数后取平均值作为最终的最优参数，在线信号按最终的最优参数进行分解，寻找最佳IMF分量并提取其时域、频域、时频域特征后，用训练好的KPCA模型完成特征的融合，最后导入训练好的BP网络中得到诊断结果，并通过与未经KPCA特征融合及未优化VMD的诊断结果进行对比的方式来证明本文方法的有效性。

1 基于综合评价指标的VMD参数优化

1.1 VMD算法原理

对信号采用变分模态分解要先建立变分模型，通过迭代更新找到变分模型的最优解，即原始信号f分解得到的K个IMF及其中心频率的最优解。

具体步骤如下：

（1）对每个模态分量u_k(t)经过Hilbert变换得到其单边频谱后，再添加指数项e^{-j{\omega }_{k}t}，得到相应的基频带。

（2）进行高斯平滑处理，得到各个模态分量u_k(t)的估计带宽，要求每个模态分量的估计带宽之和最小且所有的模态分量之和等于原始信号。根据上述要求建立的变分模型为

\left\{\begin{array}{l}\underset{\{u_k\},\{{\omega }_k\}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{\displaystyle \sum _{k=\mathrm{1}}^K}{‖{\partial }_t\left[\left(\delta (t)+\frac{j}{\mathrm{\pi }t}\right)\mathrm{*}{u}_k(t)\right]e^{-j{\omega }_{k}t}‖}_{\mathrm{2}}^{\mathrm{2}}\right\}\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{k=\mathrm{1}}^K}{u}_k(t)=f\end{array}\right.（1）

式中：\delta \left(t\right)为脉冲函数；j为复数单位；t为时间；“\mathrm{*}”为卷积符号；｛u_k｝=｛u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_K｝为VMD分解得到的K个模态分量；｛{\omega }_k｝=｛{\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}},\cdots ,{\omega }_K｝为每个模态分量对应的中心频率；‖‖{}_{\mathrm{2}}{}^{\mathrm{2}}为求二范数的平方；f为待分解的原始信号。

（3）求解变分模型时，为了求解式（1）的约束最优问题，先引入拉格朗日（Lagrange）算子\lambda和惩罚因子\alpha，将约束最优问题转化为无约束问题^［8］。Lagrange增广函数为

L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha {\displaystyle \sum _{k=\mathrm{1}}^K}{‖{\partial }_t\left[\left(\delta (t)+\frac{j}{\mathrm{\pi }t}\right)\mathrm{*}{u}_k(t)\right]e^{-j{\omega }_{k}t}‖}_{\mathrm{2}}^{\mathrm{2}}+{‖f(t)-{\displaystyle \sum _{k=\mathrm{1}}^K}{u}_k(t)‖}_{\mathrm{2}}^{\mathrm{2}}+〈\lambda (t),f(t)-{\displaystyle \sum _{k=\mathrm{1}}^K}{u}_k(t)〉（2）

为了求出增广Lagrange函数的鞍点，求解过程中将｛u_k^{\mathrm{1}}｝、｛{\omega }_{{}_k}^{\mathrm{1}}｝、{\lambda }^{\mathrm{1}}初始化后输入，经过交替方向乘子算法迭代更新各个模态分量及其对应的中心频率。

根据上述过程可以得出，变分模态分解的效果与模态分量个数K和惩罚因子\alpha有关，因此K和\alpha的选择是至关重要的。

1.2 综合评价指标及适应度函数

针对VMD分解参数难以确定的问题，本文采用鲸鱼优化算法优化VMD的分解层数K和惩罚因子\alpha，优化过程中的关键问题是适应度函数的构建。本文通过结合包络熵与峭度，提出一种新的综合评价指标，从而对各IMF分量进行评价，计算各IMF分量的综合评价指标值，并将其最小值作为寻优的适应度值。

包络熵是一种用于描述信号复杂程度的物理量。包络熵值越小，信号的复杂度就越小，自我相似性就越高^［15］。

峭度描述信号的冲击特性。峭度值越大，冲击特性越强，故障特征信息越丰富，也更容易提取到故障特征。

为使VMD分解后的IMF分量包含较多的特征信息，本文提出将这两种指标进行结合，构建一种新的综合评价指标C，对各IMF分量进行评价，其计算公式为

C_k=E_p(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)+\mathrm{1}/\mathrm{K}\mathrm{u}\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{s}(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)（3）

式中：E_p(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)表示第k（k=\mathrm{1,2},\cdots ,K）个IMF分量的包络熵值，其计算公式如式（4）所示。

E_p(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)=-{\displaystyle \sum _{i=\mathrm{1}}^R}\left(\frac{\epsilon (j)}{{\displaystyle \sum _{i=\mathrm{1}}^R}\epsilon (i)}\right)\mathrm{l}\mathrm{g}\left(\frac{\epsilon (j)}{{\displaystyle \sum _{i=\mathrm{1}}^R}\epsilon (i)}\right)（4）

式中：R为IMF分量序列的长度；\epsilon为第k个IMF分量的Hilbert包络信号。

\mathrm{K}\mathrm{u}\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{s}(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)为第k（k=\mathrm{1,2},\cdots ,K）个IMF分量的峭度值，其计算公式如式（5）所示。

\mathrm{K}\mathrm{u}\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{s}(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)=\frac{\mathrm{E}(\mathrm{I}\mathrm{M}{\mathrm{F}}_k-{\mathit{\mu }}_k{)}^{\mathrm{4}}}{{\sigma }_k^{\mathrm{4}}}（5）

式中，E（ ）表示数学期望；{\mathit{\mu }}_k为IMF分量的均值向量；{\sigma }_k^{}为IMF分量的标准差。

显然，C_k越小，对应的IMF分量冲击性和稳定性越好，故障特征信息越丰富。因此，将适应度函数定义为

h=\mathrm{m}\mathrm{i}\mathrm{n}(C_{\mathrm{1}},C_{\mathrm{2}},\cdots ,C_K)（6）

1.3 基于WOA的VMD参数优化流程

基于鲸鱼算法的VMD参数优化流程如图1所示，具体步骤如下：

步骤1：初始化。导入齿轮每种状态下的历史振动数据，设置VMD参数K和\alpha的最大值K_{\mathrm{m}\mathrm{a}\mathrm{x}}、{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}和最小值K_{\mathrm{m}\mathrm{i}\mathrm{n}}、{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}，以及鲸鱼种群的个体数量I、迭代次数t_{\mathrm{m}\mathrm{a}\mathrm{x}}、当前迭代次数t（t的初始值为0）等参数及初始化鲸鱼种群中个体的位置向量\mathit{X}={\left[K,\alpha \right]}^{\mathrm{T}}，其初始化公式为

{\mathit{X}}_i^{\mathrm{0}}=\mathit{r}\cdot (\mathit{u}\mathit{b}-\mathit{l}\mathit{b})+\mathit{l}\mathit{b}\hspace{1em}i=\mathrm{1,2},\cdots ,I（7）

式中：\mathit{r}为0~1的随机向量；参数边界最小值\mathit{l}\mathit{b}={\left[K_{\mathrm{m}\mathrm{i}\mathrm{n}},{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right]}^{\mathrm{T}}，最大值\mathit{u}\mathit{b}={\left[K_{\mathrm{m}\mathrm{a}\mathrm{x}},{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]}^{\mathrm{T}}；“∙”为两个向量对应元素的乘法。

步骤2：确定初值最优个体{\mathit{X}}_{\mathrm{0}}^{\mathrm{*}}。历史振动数据分别按照种群中每个个体的位置{\mathit{X}}_i^{\mathrm{0}}进行VMD分解。并根据式（4）、式（7）计算每个个体的适应度值h_i^{\mathrm{0}}，选取适应度值最小h_{\mathrm{0}}^{\mathrm{*}}=\mathrm{m}\mathrm{i}\mathrm{n}(h_{\mathrm{1}}^{\mathrm{0}},h_{\mathrm{2}}^{\mathrm{0}},\cdots ,h_I^{\mathrm{0}})的个体作为初始最优个体{\mathit{X}}_{\mathrm{0}}^{\mathrm{*}}。

步骤3：种群位置更新。根据［0，1］的随机数p及系数向量 A，每个个体按照以下3种方式对位置进行更新。当随机数，且向量 A 的模值时，选择最佳搜索代理方式进行收缩包围，其第i个个体位置{\mathit{X}}_i^t的更新公式为

{\mathit{X}}_i^{t+\mathrm{1}}={\mathit{X}}_t^{\mathrm{*}}-\mathit{A}\cdot \left|\mathit{C}\cdot {\mathit{X}}_t^{\mathrm{*}}-{\mathit{X}}_i^t\right|（8）

式中：系数向量\mathit{A}=a(\mathrm{2}{\mathit{r}}_{\mathrm{1}}-\mathrm{1})，\mathit{C}=\mathrm{2}{\mathit{r}}_{\mathrm{2}}；a=\mathrm{2}-\mathrm{2}t/t_{\mathrm{m}\mathrm{a}\mathrm{x}}，在迭代过程中从2线性衰减至0；{\mathit{r}}_{\mathrm{1}}、{\mathit{r}}_{\mathrm{2}}均为［0，1］中的随机向量；| |表示对向量中的每个元素取绝对值；{\mathit{X}}_t^{\mathrm{*}}为当前最优个体位置；\left|\mathit{C}\cdot {\mathit{X}}_t^{\mathrm{*}}-{\mathit{X}}_i^t\right|为猎物与当前鲸鱼个体间的相对距离。 1 表示单位列向量。当且\left|\mathit{A}\right|\ge \mathrm{1}时，选择随机搜索代理方式进行迭代更新。个体位置{\mathit{X}}_i^t的更新公式为

{\mathit{X}}_i^{t+\mathrm{1}}={\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}-\mathit{A}\cdot \left|\mathit{C}\cdot {\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}-{\mathit{X}}_i^t\right|（9）

式中：{\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}为随机选取鲸鱼个体的位置向量； \left|\mathit{C}\cdot {\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}-{\mathit{X}}_i^t\right|代表当前鲸鱼个体与其他任意鲸鱼个体的相对距离。当p\ge \mathrm{0.5}时，选择螺旋收缩方式进行迭代更新。个体位置{\mathit{X}}_i^t的更新公式为

{\mathit{X}}_i^{t+\mathrm{1}}={\mathit{e}}^{\mathit{b}\mathit{l}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\mathrm{\pi }l)\left|{\mathit{X}}_t^{\mathrm{*}}-{\mathit{X}}_i^t\right|+{\mathit{X}}_t^{\mathrm{*}}（10）

式中：b为定义对数螺线形状的常数；l为［-1，1］的随机数；\left|{\mathit{X}}^{\mathrm{*}}(t)-{\mathit{X}}_i(t)\right|代表当前鲸鱼个体与最优鲸鱼个体间的距离。

步骤4：计算更新后各个体{\mathit{X}}_i^{t+\mathrm{1}}的适应度值，确定当前的最优个体{\mathit{X}}_{t+\mathrm{1}}^{\mathrm{*}}及其适应度值h_{t+\mathrm{1}}^{\mathrm{*}}。当h_{t+\mathrm{1}}^{\mathrm{*}}\le h_t^{\mathrm{*}}时，{\mathit{X}}_{t+\mathrm{1}}^{\mathrm{*}}为最优个体；当h_{t+\mathrm{1}}^{\mathrm{*}}>h_t^{\mathrm{*}}时，则{\mathit{X}}_{t+\mathrm{1}}^{\mathrm{*}}={\mathit{X}}_t^{\mathrm{*}}（最优保留）。

步骤5：判断是否达到最大迭代次数，若满足则循环结束，将当前最优个体{\mathit{X}}_{t+\mathrm{1}}^{\mathrm{*}}作为最优参数结果输出。否则，返回步骤3继续循环迭代。

2 基于KPCA与BP网络的故障诊断方法

根据优化得到的最优参数进行VMD分解后，为了减少数据的维度并同时保留尽可能多的有用信息，本文采用基于KPCA与BP网络相结合的故障诊断方法。

2.1 IMF分量筛选及特征提取

如果确定好最优VMD参数K^{\mathrm{*}}、{\alpha }^{\mathrm{*}}后，就可以对齿轮振动样本数据进行VMD分解，得到K^{\mathrm{*}}个IMF分量。本文根据式（4）计算每个分量的综合评价指标，并选取最小的综合评价指标所对应的IMF分量作为最优分量\mathrm{I}\mathrm{M}{\mathrm{F}}^{\mathrm{*}}。

针对单一特征指标无法准确全面描述振动信号的故障特征，本文提出对最优IMF分量\mathrm{I}\mathrm{M}{\mathrm{F}}^{\mathrm{*}}提取时域（如方差、均值、方根幅值、均方根值等）、频域（平均频率、重心频率等）、时频域（能量熵等）特征。

2.2 基于KPCA的特征融合模型构建

KPCA的核心思想是通过核函数将原始空间中的数据通过核函数映射到高维空间中，使低维空间中线性不可分的数据在高维空间中线性可分^［16］，再通过主成分分析（PCA）将高维空间中的数据映射到另一个低维空间中，从而实现数据非线性特征融合。

应用KPCA对特征进行融合的关键问题是确定KPCA模型，即利用历史数据确定变换矩阵。构建KPCA模型的主要步骤如下：

步骤1：构建特征样本集。对齿轮历史数据（包括正常状态和故障状态）提取特征向量\mathit{s}，构建特征数据集S=\{{\mathit{s}}_{\mathrm{1}},{\mathit{s}}_{\mathrm{2}},\cdots ,{\mathit{s}}_M\}，其中M为样本数。

步骤2：标准化处理。运用特征样本均值向量\overline{\mathit{s}}=\frac{\mathrm{1}}{M}{\displaystyle \sum _{i=\mathrm{1}}^M}{\mathit{s}}_i和特征样本方差{\mathit{\sigma }}^{\mathrm{2}}=\frac{\mathrm{1}}{M}{\displaystyle \sum _{i=\mathrm{1}}^M}{\left({\mathit{s}}_i-\overline{\mathit{s}}\right)}^{\mathrm{2}}对特征集S中的样本进行标准化处理，处理后的数据集记为\widehat{S}=\{{\widehat{\mathit{s}}}_{\mathrm{1}},{\widehat{\mathit{s}}}_{\mathrm{2}},\cdots ,{\widehat{\mathit{s}}}_M\}。

步骤3：构建核矩阵 Q。使用样本特征向量集\widehat{\mathit{S}}构建核矩阵 Q，其表达式如下

\mathit{Q}={\left[\begin{array}{ccc}Q({\widehat{\mathit{s}}}_{\mathrm{1}},{\widehat{\mathit{s}}}_{\mathrm{1}})& \cdots & Q({\widehat{\mathit{s}}}_{\mathrm{1}},{\widehat{\mathit{s}}}_M)\\ ⋮& ⋮& ⋮\\ Q({\widehat{\mathit{s}}}_M,{\widehat{\mathit{s}}}_{\mathrm{1}})& \cdots & Q({\widehat{\mathit{s}}}_M,{\widehat{\mathit{s}}}_M)\end{array}\right]}_{M\times M}（11）

当选择高斯径向基核函数时，核矩阵元素Q({\widehat{\mathit{s}}}_i,{\widehat{\mathit{s}}}_j)计算公式为

Q({\widehat{\mathit{s}}}_i,{\widehat{\mathit{s}}}_j)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{‖{\widehat{\mathit{s}}}_i-{\widehat{\mathit{s}}}_j‖}_{}^{{}^{\mathrm{2}}}}{\mathrm{2}{\gamma }^{\mathrm{2}}}\right)（12）

式中：\gamma为核宽度。

步骤4：确定KPCA变换矩阵 W。对核矩阵 Q 进行特征值分解，将特征值按从大到小排列，即{\lambda }_{\mathrm{1}}\ge {\lambda }_{\mathrm{2}}\ge \cdots \ge {\lambda }_M，对应的特征向量为{\mathit{\alpha }}_{\mathrm{1}}，{\mathit{\alpha }}_{\mathrm{2}}，\cdots，{\mathit{\alpha }}_M。若前q个主元的累计贡献率{\displaystyle \sum _{i=\mathrm{1}}^q}{\lambda }_i/{\displaystyle \sum _{i=\mathrm{1}}^M}{\lambda }_i达到预设要求（通常取85%）时，KPCA变换矩阵 W 为

\mathit{W}={\left[\begin{array}{cccc}{\mathit{\alpha }}_{\mathrm{1}}& {\mathit{\alpha }}_{\mathrm{2}}& \cdots & {\mathit{\alpha }}_q\end{array}\right]}_{M\times q}（13）

显然，建立KPCA特征融合模型就是确定样本均值向量\overline{\mathit{s}}、特征样本方差{\mathit{\sigma }}^{\mathrm{2}}和变换矩阵 W。对于任一标准化后的特征样本\widehat{\mathit{s}}，其在高维空间的核向量\mathit{v}={\left[\begin{array}{cccc}Q({\widehat{\mathit{s}}}_{\mathrm{1}},\widehat{\mathit{s}})& Q({\widehat{\mathit{s}}}_{\mathrm{2}},\widehat{\mathit{s}})& \cdots & Q({\widehat{\mathit{s}}}_M,\widehat{\mathit{s}})\end{array}\right]}^{\mathrm{T}}，则\widehat{\mathit{s}}的主元得分（即融合后的特征向量）\mathit{t}={\left[\begin{array}{cccc}{t}_{\mathrm{1}}& t_{\mathrm{2}}& \cdots & t_q\end{array}\right]}^{\mathrm{T}}为

\mathit{t}={\mathit{W}}^{\mathrm{T}}\mathit{v}（14）

2.3 基于BP网络的故障诊断模型构建

为了实现对滚动轴承故障的识别，本文构建一种基于BP网络的故障诊断模型，其结构如图2所示。输入层为主元得分向量t，即融合后的特征向量；输出层为故障类别，假设齿轮共正常、断齿、裂纹3种状态，则输出层节点数为3，可将正常状态标记为100；断齿状态标记为010；裂纹状态标记为001。

训练时，带标签的输入训练集 TT 是由标准化后的特征样本{\widehat{\mathit{s}}}_i在核主元空间的得分向量{\mathit{t}}_i(i=\mathrm{1,2},\cdots ,M)构成，即

\mathit{T}{\mathit{T}}_{M\times q}={\left[\begin{array}{cccc}{\mathit{t}}_{\mathrm{1}}& {\mathit{t}}_{\mathrm{2}}& \cdots & {\mathit{t}}_M\end{array}\right]}^{\mathrm{T}}=({\mathit{W}}^{\mathrm{T}}{\mathit{Q})}^{\mathrm{T}}=\mathit{Q}\mathit{W}（15）

2.4 故障诊断流程

对一个在线样本进行故障诊断的具体步骤为：

步骤1：将该在线样本按寻优好的最优参数[K^{\mathrm{*}},{\alpha }^{\mathrm{*}}]进行VMD分解。

步骤2：筛选出最优IMF分量并对其进行特征提取，得到特征向量{\mathit{s}}_{\mathrm{n}\mathrm{e}\mathrm{w}}。

步骤3：标准化处理。通过\overline{\mathit{s}}和{\mathit{\sigma }}^{\mathrm{2}}对{\mathit{s}}_{\mathrm{n}\mathrm{e}\mathrm{w}}进行标准化处理，处理后的特征向量记为{\widehat{\mathit{s}}}_{\mathrm{n}\mathrm{e}\mathrm{w}}。

步骤4：根据式（14）计算主元得分{\mathit{t}}_{\mathrm{n}\mathrm{e}\mathrm{w}}，即：{\mathit{t}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathit{W}}^{\mathrm{T}}{\mathit{v}}_{\mathrm{n}\mathrm{e}\mathrm{w}}。其中{\mathit{v}}_{\mathrm{n}\mathrm{e}\mathrm{w}}为

{\mathit{v}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\left[Q({\widehat{\mathit{s}}}_{\mathrm{1}},{\widehat{\mathit{s}}}_{\mathrm{n}\mathrm{e}\mathrm{w}})Q({\widehat{\mathit{s}}}_{\mathrm{2}},{\widehat{\mathit{s}}}_{\mathrm{n}\mathrm{e}\mathrm{w}})\cdots Q({\widehat{\mathit{s}}}_M,{\widehat{\mathit{s}}}_{\mathrm{n}\mathrm{e}\mathrm{w}})\right]}^{\mathrm{T}}（16）

步骤5：将经过特征融合后的特征向量{\mathit{t}}_{\mathrm{n}\mathrm{e}\mathrm{w}}导入训练好的BP网络中，得到该样本的诊断结果。

3 实例验证

3.1 原始数据

本文的试验数据是利用江苏千鹏QPZZ-II型旋转机械故障模拟试验台测得，如图3所示。其中，齿轮箱由55齿的主动齿轮和75齿的从动齿轮组成，采用电子放电的方法将裂纹和断齿故障分别引入从动齿轮，齿轮箱顶部安装传感器。

本试验台的数据采样频率设置为5 kHz。采集到的数据加入-10 dB的高斯白噪声后，每种工作状态中选取100组样本数据，每组样本包括1 024个采样点。

3.2 WOA优化VMD

本文分别对3种不同状态的齿轮历史数据进行VMD优化，其中WOA种群规模为10；最大迭代次数为10；K的取值范围为［3，10］，且为整数；α的取值范围为［100，2 500］。优化过程中随着迭代次数的增加而得到的收敛曲线如图4所示，优化后的结果如表1所示。

同时，为确保实际工程意义，使最后得到的最优参数对不同状态的齿轮数据具有通用性，故对不同状态下最优的K和α取平均值作为最终的VMD参数。由表1可知，K^* =10，α^*=2 499。

3.3 KPCA特征融合

利用最优参数[K^{\mathrm{*}},{\alpha }^{\mathrm{*}}]对每个样本进行VMD分解后，根据综合评价指标筛选出最优的IMF分量。以正常状态的前3个样本为例，其对应的IMF分量的适应度值如表2所示。

在表2中，样本1的IMF10为最优分量，对应的适应度值为9.43；样本2的IMF2为最优分量，对应的适应度值为9.28；样本3的IMF3为最优分量，对应的适应度值为9.35。

对最优IMF分量进行时域、频域、时频域特征的计算，综合考虑各个特征指标在振动信号敏感度和稳定性方面的优势，本文选择方差、均值等13个时域特征，平均频率、重心频率等4个频域特征，以及1个能反映能量分布均匀性且可体现IMF分量所含频率复杂程度的时频域特征-能量熵，从而构建多域特征集。每种状态取50个样本的特征集作为训练集，每种状态取另50个样本的特征集作为测试集。

在构建多域特征集后，先进行KPCA模型的训练。其中的主元累积贡献率分布如图5所示，图中只列出了前10个主元。由图5可知，前3个主元的累积贡献率已经达到了95.00%以上，故选取前3个主元，将特征融合成三维，KPCA模型训练完毕。

测试集经过KPCA模型进行特征融合后的三维特征分布如图6所示。由图6可知，3种状态类内的聚类特性较好，类间的差异性距离较大，特征分类较明显。

3.4 BP网络诊断

先进行BP网络模型的训练。对训练集经KPCA模型后得到的融合特征集按状态划分为1、2、3共3种标签，其中1对应100；2对应010，3对应001。对带标签的融合特征集完成BP模型的训练。

测试集经KPCA模型后再输入到训练好的BP网络模型中进行故障诊断和分类识别，得到的测试集分类诊断结果如图7所示。由图7可知，诊断准确率高达98.00%。

3.5 对比试验

为验证本文所提出模型的准确性及适用性，现将本文提出的模型与中心频率法^［17］VMD-KPCA-BP及WOA-VMD-BP两种模型进行对比分析。

当采用中心频率法VMD-KPCA-BP模型时，3种状态下经中心频率法确定的分解层数K值均为3，惩罚因子选用默认值α=2 000，进而进行VMD分解，对最优IMF分量提取上文所述的多域特征后，再利用此特征集通过上文方法训练KPCA特征融合模型和BP网络诊断模型。中心频率法VMD-KPCA的特征分布如图8所示。

对比图6、图8可知，本文提出的WOA-VMD-KPCA模型对于特征分类取得了更显著的效果。测试集的分类诊断结果如图9所示。

由图7和图9可知，中心频率法VMD-KPCA-BP模型的诊断准确率为90.67%，低于本文所提出模型准确率，证明了基于综合评价指标的WOA-VMD方法的有效性。

当采用WOA-VMD-BP模型时，利用通过WOA算法得到的最优参数进行VMD分解，对最优IMF分量提取上文所述的多域特征后，直接利用此特征集通过上文方法构建BP网络诊断模型，测试集的分类诊断结果如图10所示。

由图10可知，WOA-VMD-BP模型的诊断准确率为94%，低于本文所提出模型准确率，证明了KPCA特征融合方法的有效性。

3.6 对比结果

为了更直观地体现本文模型的优势，将本文模型与中心频率法VMD-KPCA-BP模型及WOA-VMD-BP模型测试集的诊断准确率进行对比，结果如表3所示。由表3可知，本文提出的WOA-VMD-KPCA-BP模型在故障诊断过程中得到了更高的诊断准确率，取得了更显著的效果，可为后续的研究及应用提供参考。

4 结论

本文提出了一种将基于综合评价指标的VMD参数寻优、KPCA特征融合和BP网络相结合的齿轮故障诊断方法，并用实例进行了分析和验证。

（1）提出构建基于包络熵与峭度的综合评价指标并采用WOA算法对VMD参数进行寻优。

（2）提出对不同状态下最优的K和α取平均值，最终得到对任何状态下都适用的最优参数。用最优参数进行VMD分解后再利用同样的综合评价指标筛选出最优IMF分量，提取多域高维特征后经过KPCA模型完成特征的融合，最后导入BP网络中进行故障诊断。

通过对比分析，中心频率法VMD-KPCA-BP模型的准确率和WOA-VMD-BP模型的准确率分别为90.67%及94.00%，而本文所提出方法的诊断准确率为98.00%。因此，本文方法能更有效地提取特征向量，且具有更高的诊断准确率，更适用于齿轮的故障诊断。

参考文献

[1] Chen H Z， Wang J X， Tang B P，et al.An integrated approach to planetary gearbox fault diagnosis using deep belief networks［J］.Measurement Science and Technology，2017，28（2）：025010.

[2] 刘鲲鹏，白云川，吕麒鹏，等.基于改进傅里叶变换的滚动轴承故障诊断［J］.内燃机与配件，2018（19）：56-57.

[3] 文成，周传德.基于灰色关联度的形态滤波及滚动轴承故障诊断中应用［J］.振动与冲击，2015，34（14）：51-55.

[4] 杨超，杨晓霞.基于灰色关联度和Teager能量算子的轴承早期故障诊断［J］.振动与冲击，2020，39（13）：224-229.

[5] 涂志松.基于小波包和支持向量机的滚动轴承故障诊断［J］.机电工程技术，2020，49（12）：208-211.

[6] 张琛，郭俊超，甄冬，等.基于小波包时延相关解调的滚动轴承故障诊断方法［J］.机械设计，2020，37（6）：24-28.

[7] 戚晓利，叶绪丹，蔡江林，等.基于变分模态分解与流形学习的滚动轴承故障特征提取方法［J］.振动与冲击，2018，37（23）：133-140.

[8] Dragomiretskiy K， Zosso D.Variational mode decomposition［J］.IEEE Transactions on Signal Processing，2014，62（3）：531-544.

[9] 武英杰，甄成刚，刘长良.变分模态分解在风电机组故障诊断中的应用［J］.机械传动，2015，39（10）：129-132.

[10] 徐甜甜，孙文磊，王宏伟，等.变分模态分解与谱峭度法在齿轮箱故障诊断中的应用［J］.机械科学与技术，2020，39（6）：873-878.

[11] Zhang Y G， Li R X， Zhang J H.Optimization scheme of wind energy prediction based on artificial intelligence［J］.Environmental Science and Pollution Research，2021，28（29）：39966-39981.

[12] 刘秀丽，徐小力，吴国新，等.基于变分模态分解的故障弱信息提取方法［J］.华中科技大学学报（自然科学版），2020，48（7）：117-121.

[13] 张萍，张文海，赵新贺，等.WOA-VMD算法在轴承故障诊断中的应用［J］.噪声与振动控制，2021，41（4）：86-93，275.

[14] Mirjalili S， Lewis A.The whale optimization algorithm［J］.Advances in Engineering Software，2016，95（3）：51-67.

[15] 尚晓三，王式成，王振龙，等.基于样本熵理论的自适应小波消噪分析方法［J］.水科学进展，2011，22（2）：182-188.

[16] Vinay A， Shekhar V S， Murthy K，et al.Face recognition using gabor wavelet features with PCA and KPCA-A comparative study［J］.Procedia Computer Science，2015，57（12）：650-659.

[17] 陈东宁，张运东，姚成玉，等.基于变分模态分解和多尺度排列熵的故障诊断［J］.计算机集成制造系统，2017，23（12）：2604-2612.

基金资助

国家自然科学基金(62003223)