基于WP-TRP的滚动轴承故障诊断方法

Fault Diagnosis Method for Rolling Bearings Based on WP-TRP

摘要

针对故障诊断中传统时频域法提取特征时易受主观因素影响而导致冗余，且深度学习算法受训练数据影响导致计算复杂性较高的缺点，将时域和频域结合，提出一种基于小波包-无阈值递归图（WP-TRP）的滚动轴承故障诊断方法.首先，提出递减信息熵准则，以克服小波包分解的主观性，获取更准确的时频域特征；在此基础上，引入无阈值递归图思想，充分提取数据初始时域特征，并利用奇异值分解进一步降低冗余特征，提高计算效率.然后，引入海洋捕食者算法来获得支持向量机最优参数，实现故障诊断的准确分类.最后，通过标准滚动轴承数据集仿真验证了所提方法的有效性.

Abstract

In fault diagnosis， the traditional time-frequency domain methods are easily affected by subjective factors while being used for feature extraction， so that the redundancy emerges. Deep learning algorithm is highly dependent on training data and has computation complexity. Fault diagnosis method for rolling bearings based on wavelet packet-thresholdless recurrence plot （WP-TRP）is proposed by combining time with frequency domains. Firstly， the decreasing information entropy criterion is developed to overcome the subjectivity of wavelet packet decomposition for acquisition of more accurate time-frequency feature. On this basis， the idea of thresholdless recurrence plot is introduced to extract the initial time domain feature. Moreover， by adopting the singular value decomposition to decrease the redundant feature， the computational efficiency can be increased. Secondly， the marine predator algorithm is introduced to obtain the optimal parameters of supporting vector machine， by which the more accurate classification can be realized. Finally， the effectiveness of the presented method is verified by using the simulation on the benchmark rolling bearing datasets.

Graphical abstract

关键词

故障诊断 / 小波包分解 / 信息熵 / 无阈值递归图 / 奇异值分解 / 海洋捕食者算法

Key words

fault diagnosis / wavelet packet decomposition / information entropy / thresholdless recurrence plot（TRP） / singular value decomposition / marine predator algorithm

对复杂部件进行故障诊断，符合工业过程的安全运行要求，具有重要的科学研究意义和市场应用价值^［1-2］.目前故障诊断主要应用传统时频域^［3］、深度学习^［4］等方法，从原始输入、输出数据中提取有效特征，并与后续分类器结合实现准确的故障分类.因此，特征能否被充分提取直接影响所建故障模型的准确性.传统方法如经验模态分解^［3］、高分模态分解等易受到信号的非线性影响，难以提取有效的特征.小波及小波包分解法^［5］能够自适应处理高/低频信号，但仍须提前确定分解层数，主观性较强.此外，传统小波包分解（WPD）法提取特征多为一维数据，难以凸显其非线性特性.相对传统方法，深度学习法的端到端策略能够直接提取特征，如循环神经网络（RNN）^［6］和卷积神经网络（CNN）^［4］.RNN主要通过门控单元提取一维数据特征，具有记忆性，能够有效处理时序数据，但在网络训练中易出现参数的梯度爆炸和消失现象，使网络训练精度不足，导致所生成的故障特征提取不够充分；CNN主要采用卷积和池化层提取特征，能有效处理一维和二维图像数据，训练速度更快，但随层数的加深，易出现过拟合，使卷积神经网络对未建模数据的泛化能力变差，尤其是多工况下数据分布复杂时，可能会直接影响生成故障特征的精确性，导致故障诊断精度降低.

由于直接将一维原始数据作为输入会导致特征提取不充分，因此可将一维数据转化成二维图像数据，如格拉姆角场（GAF）^［7］和递归图（RP）^［8］.GAF的优势在于在极坐标系中可将一维时序数据的幅值转成角度，间接保留时间信息；RP利用相空间重构，可直接从时序数据内部获得先验知识，保留了时间信息，且可揭示故障的动力学特性，更加直观具体.但RP难以准确描述长时间序列的特征，且阈值需要人为设定，主观性较强，导致所产生特征易受到阈值影响而难以充分提取.二维图像数据虽然可以提供比一维数据更丰富和详细的信息表达方式，但二维图像的数据维度较大，计算效率较低.

综上所述，本文提出了一种基于小波包-无阈值递归图（WP-TRP）的滚动轴承故障诊断方法.首先引入递减信息熵准则来实现小波包（WP）分解层数的选取；其次利用无阈值递归图（TRP）实现从一维到二维图像特征转换，并采用奇异值分解（SVD）获取图像的核心特征；最后利用海洋捕食者算法（MPA）优化的支持向量机（SVM）实现故障分类.

1 基于WP-TRP的特征提取

要获得较为精确的故障诊断结果，应构造有效的初始特征集，来尽可能地覆盖待诊断对象的非线性特征.为此，首先利用WPD并基于所提的递减信息熵准则获得初始特征集，并与TRP结合，利用其获得更有效的二维图像特征；此外，针对上述TRP生成的特征，利用SVD进一步降低冗余特征，获得更精简的一维特征，与二维特征相比，更具简单性，能提高计算效率；最后，将所获得的一维特征输入到SVM，并利用MPA来优化SVM的最优参数，实现最终的故障诊断分类，进一步提高故障诊断模型的准确性.

1.1 WPD的基本原理

针对信号x(t)，小波包分解j层得

\left.\begin{array}{l}x(t)={\displaystyle \sum _{k=\mathrm{0}}^{{\mathrm{2}}^j-\mathrm{1}}}{x}_{j,k}(t);\\ x_{j,k}(t)=\sum _k{w}_{j,k,o}{\psi }_{j,k,o}(t);\\ w_{j,k,o}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}x(t){\psi }_{j,k,o}(t)\mathrm{d}t.\end{array}\right\}（1）

式中：x_{j,k}(t)为小波包组分信号；w_{j,k,o}为小波系数；{\psi }_{j,k,o}(t)为WP函数；层数j=\mathrm{1,2},\dots ,l，l为WPD的最大层数；k=\mathrm{0,1}，\dots ,{\mathrm{2}}^j-\mathrm{1}，表示第j层的第k个频段；o=\mathrm{1,2},\dots ,h，为第j层第k个频段的第o个位置指标，h为WPD的最大位置指标.

{\psi }_{j,k,o}(t)是一组标准正交基，具有正交性，即

{\psi }_{a,k,o}(t){\psi }_{b,k,o}(t)=\mathrm{0}.（2）

式中，a和b是j的不同取值.

由式（1）获得第j层的第k个频段的小波系数w_{j,k,o}，因此，小波系数向量为{\mathit{w}}_{j,k}=(w_{j,k,\mathrm{1}},w_{j,k,\mathrm{2}},\dots ,w_{j,k,h})最终形成，其中h为第j层第k个频段的最后一个位置指标.

1.2 基于递减信息熵准则的WPD层数确定

WPD的层数一般采用试凑确定，易导致特征分解不充分.本文引入小波能量谱公式来确定分解层数的范围，降低计算量，并在此范围内提出递减信息熵准则来确定最终层数，基本原理如下：

小波能量谱公式^［9］定义为

l=⌈\mathrm{l}\mathrm{g}\frac{f_{\mathrm{s}}}{\mathrm{4}f}⌉.（3）

式中：⌈\cdot ⌉为向上取整函数；f_{\mathrm{s}}=\mathrm{12}\mathrm{000} Hz，表示信号的采样频率；本文采用凯斯西储大学标准滚动轴承数据集^［10］，电机负载为0.745 7 kW，转速 v =1 772 r/min，f为故障特征频率，可取值为69.7，159.88，106.32 Hz，代入式（3）计算分别得l为4，3，4.取最大分解层数l=4.

在此基础上，提出递减信息熵准则，进一步确定最终层数.给定轴承数据集\mathit{X}=\{{\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\dots ,{\mathit{x}}_m\}，其中{\mathit{x}}_i=(x_{i,\mathrm{1}},x_{i,\mathrm{2}},\dots ,x_{i,j},\dots ,x_{i,n})为一维振动加速度数据，m为数据总数，i为数据序号，n为数据维数.对于第i个数据{\mathit{x}}_i，经过WPD后，第j层的递减信息熵定义如下：

A_{i,j}=E_{i,j}-E_{i,j-\mathrm{1}}.（4）

式中，E_{i,j}为第i个数据第j层所有频段对应小波系数的信息熵平均值，即原始信息熵^［11］准则，其定义如下：

E_{i,j}=\left(E_{i,j,\mathrm{1}}+E_{i,j,\mathrm{2}}+\dots +E_{i,j,k}+\dots +E_{i,j,{\mathrm{2}}^j}\right)/{\mathrm{2}}^j.（5）

式中，E_{i,j,k}为第i个数据第j层的第k个频段的子信息熵值.

E_{i,j,k}=-{\displaystyle \sum _{q=\mathrm{1}}^s}P(w_{i,j,k,q})lg(P(w_{i,j,k,q})).（6）

式中：w_{i,j,k,q}为第i个数据第j层第k个频段的第q个小波系数，q=\mathrm{1,2},\dots ,s，s=n/{\mathrm{2}}^j；P(w_{i,j,k,q})为w_{i,j,k,q}发生的概率.

\left.\begin{array}{c}P(w_{i,j,k,q})=\frac{(w_{i,j,k,q}{)}^{\mathrm{2}}}{\sum _i(w_{i,j,k,q}{)}^{\mathrm{2}}},\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{q=\mathrm{1}}^s}P(w_{i,j,k,q})=\mathrm{1}.\end{array}\right\}（7）

在式（4）~式（7）所述的递减信息熵函数基础上，对于第i个数据{\mathit{x}}_i，WPD的最终分解层数j_i由A_{i,j}的最小值对应的层数确定，定义为

j_i=\{j|\underset{j\in [\mathrm{1},l]}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(A_{i,j}\right)\}.（8）

对于给定数据集\mathit{X}=\{{\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\dots ,{\mathit{x}}_m\}，最终的分解层数j为

j=\underset{i\in [\mathrm{1},m]}{\mathrm{m}\mathrm{a}\mathrm{x}}\{j_i\}.（9）

本文所提的递减信息熵准则与原始信息熵准则相比，具有如下优点：

首先，应用式（5）定义的原始信息熵准则，最终分解层数会趋于一个较大值，收敛性较差.本文所提的递减信息熵准则，利用相邻两层之间的信息熵平均值之差作为当前层的最终准则值，并形成候选的集合，然后从该集合中选择衰减最大值所在的层，作为WPD的最终层数.这种计算方式能够考虑特征分解的充分性，获得适宜的WPD层数.这是因为分解不充分，即层数过少，会导致信号降噪不彻底，影响有效信息的提取；当分解层数过高时，由于信号经过较多的滤波器，导致信号的信噪比反而下降，同时会使有用信息丢失造成信号失真.

对于任意{\mathit{x}}_i\in \mathit{X}，i\in [\mathrm{1},m]，均获得对应的j_i，并利用式（9）对所有j_i进行比较，获得最终的WPD层数j.由表1所示，当取相邻的分解层为第2层、第3层时，对应的数据最多，为486个，由此确定适宜的分解层数为3.对\mathit{X}进行WPD的3层分解，得到小波系数矩阵集合\mathit{W}=\{{\mathit{W}}_{\mathrm{1}},{\mathit{W}}_{\mathrm{2}},\dots ,{\mathit{W}}_i,\dots ,{\mathit{W}}_m\}，{\mathit{W}}_i为{\mathit{x}}_i分解后对应的小波系数矩阵.

1.3 基于TRP的二维图像特征提取

经WPD处理后的数据仍为一维频域特征，由于一维频域特征所蕴含的非线性特征难以有效体现，且单一频域特征难以保证特征的充分性.本文从充分挖掘潜在特征并结合时域和频域的角度，进一步构造故障特征，为此，引入TRP^［12］的基本概念.相比传统RP^［8］而言，克服了采用阈值判断而仅仅生成二值图像的缺点，有效保留了递归点间的距离参数，从而在对同一时间序列进行特征提取时包含更多的信息量，能够生成蕴含信息值更丰富的多值图像，更能表达非线性的特征.

基于TRP提取二维图像特征，主要包括基于相空间重构的高维数据构建和基于TRP的二维图像特征生成.

1.3.1 基于相空间重构的高维数据构建

基于相空间重构^［13］的高维数据构建过程描述如下：给定一维数据集\left\{x_I\right\}，I=\mathrm{1,2},\mathrm{3},\dots ,N，经相空间重构得到高维空间的数据集{\mathit{y}}_e.

{\mathit{y}}_e=\{x_I,x_{I+\tau },x_{I+\mathrm{2}\tau },\dots ,x_{I+(d-\mathrm{1})\tau }\}.（10）

式中：{\mathit{y}}_e为重构数据集{\mathit{y}}_e的第e个数据点，e=\mathrm{1,2},\dots ,N_d，N_d=N-(d-\mathrm{1})\tau，为重构数据集{\mathit{y}}_e的数据点个数；x_{I+\tau }为原始一维数据集\left\{x_I\right\}中第\left(I+\tau \right)个数据；\tau为延迟时间；d为嵌入维数.

1.3.2 基于TRP的二维图像特征生成

传统RP矩阵中的元素R_{\mathrm{0}}(p,u,\epsilon )定义如下：

R_{\mathrm{0}}(p,u,\epsilon )=\Theta \left(\epsilon -‖{\mathit{y}}_p-{\mathit{y}}_u‖\right).（11）

式中：\epsilon为给定阈值；\Theta \left(\cdot \right)为Heaviside函数； y_p 为重构数据集 y_e 的第p个数据点， y_u 为重构数据集 y_e 的第u个数据点，p，u=1，2，…，N_d .

（12）

对数据{\mathit{y}}_p和{\mathit{y}}_u，TRP中的矩阵元素定义如下：

R(p,u)=‖{\mathit{y}}_p-{\mathit{y}}_u‖₂.（13）

式中，‖\cdot ‖₂为2范数.

由式（11）可见，传统RP公式的结果只能是0或1.与式（11）不同，式（13）中没有阈值\epsilon的限定，R(p,u)定义为相空间中向量{\mathit{y}}_p和{\mathit{y}}_u之间的欧氏距离，使输出结果为连续值，增强了对特征的描述，最终可求得无阈值递归矩阵 R 为

\mathit{R}=\left[\begin{array}{cccccc}R(\mathrm{1},\mathrm{1})& R(\mathrm{1},\mathrm{2})& \dots & R(\mathrm{1},u)& \dots & R(\mathrm{1},N_d)\\ R(\mathrm{2},\mathrm{1})& R(\mathrm{2},\mathrm{2})& \dots & R(\mathrm{2},u)& \dots & R(\mathrm{2},N_d)\\ ⋮& ⋮& & ⋮& & ⋮\\ R(p,\mathrm{1})& R(p,\mathrm{2})& \dots & R(p,u)& \dots & R(p,N_d)\\ ⋮& ⋮& & ⋮& & ⋮\\ R(N_d,\mathrm{1})& R(N_d,\mathrm{2})& \dots & R(N_d,u)& \dots & R(N_d,N_d)\end{array}\right].（14）

图1为给定一维数据集下，用本文所提方法得到的小波系数图（WCP）、传统RP和TRP转化后的图像特征.由图可知，传统RP仅有0和1的元素值，而TRP的元素值为一个区间内的值，因此包含比传统RP更多的信息量，更能准确描述数据非线性特征.

从矩阵 R 可见，主对角线上的元素由于两两元素相等，即p=u，则对应元素为零，即R\left(p,u\right)=\mathrm{0}.与主对角线相邻的两个位置对称元素的值也相等，则递归图矩阵可简化如下：

\mathit{R}=\left[\begin{array}{cccccc}\mathrm{0}& R(\mathrm{1},\mathrm{2})& \dots & R(\mathrm{1},u)& \dots & R(\mathrm{1},N_d)\\ R(\mathrm{1},\mathrm{2})& \mathrm{0}& \dots & R(\mathrm{2},u)& \dots & R(\mathrm{2},N_d)\\ ⋮& ⋮& & ⋮& & ⋮\\ R(\mathrm{1},p)& R(\mathrm{2},p)& \dots & \mathrm{0}& \dots & R(N_d,p)\\ ⋮& ⋮& & ⋮& & ⋮\\ R(\mathrm{1},N_d)& R(\mathrm{2},N_d)& \dots & R(N_d,u)& \dots & \mathrm{0}\end{array}\right].（15）

由式（15）可见，TRP矩阵 R 为对称矩阵，以主对角线为对称轴的上下位置的元素相同，因此用于计算时会存在较大的冗余，导致计算效率低下.

1.4 基于SVD的特征约简

针对TRP所提取特征的冗余，会导致计算量增大，并可能导致后续诊断精度降低的问题，本文利用奇异值分解来进一步提取TRP所产生二维图像的核心信息，并将其转换为一维数据，以有效保留有用信息，并提高计算效率.SVD的基本原理^［14］描述如下：

对矩阵 A，若正交矩阵\mathit{U}\in {\mathit{R}}_{m\times m}和\mathit{V}\in {\mathit{R}}_{n\times n}，使\mathit{A}=\mathit{U}\mathit{D}{\mathit{V}}^{\mathrm{T}}，其中\mathit{D}=\left[\begin{array}{cc}\mathit{S}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right]，\mathit{S}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\sigma }_{\mathrm{1}},{\sigma }_{\mathrm{2}},\dots ,{\sigma }_r)，且σ₁≥σ₂≥\dots≥σ_r >0，则{\sigma }_{\mathrm{1}},{\sigma }_{\mathrm{2}},\dots ,{\sigma }_r即为矩阵\mathit{A}的奇异值，r为奇异值个数.

SVD的本质在于有效确定奇异值的个数r，但传统方法通常通过试凑给定r，具有较多主观不确定性，为此，本文引入奇异值差分谱^［15］，利用奇异值差分谱的最大峰值来确定r.

对轴承数据集\mathit{X}转化后的TRP矩阵 R 进行奇异值分解，并得到奇异值差分谱，得到峰值的位置都为1，因此本文中r=\mathrm{1}.通过SVD进行特征提取后，原本\mathrm{8}\times \mathrm{120}\times \mathrm{120}的三维特征被约简到\mathrm{8}\times \mathrm{1}维，大大降低了识别复杂度和计算量.

1.5 基于WP-TRP的特征提取算法描述

基于WP-TRP的特征提取算法步骤如下：

步骤1 给定轴承一维振动加速度数据集\mathit{X}=\{{\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\dots ,{\mathit{x}}_m\}，根据本文提出的递减信息熵准则，确定适宜的分解层数为j=\mathrm{3}.

步骤2 对每段数据{\mathit{x}}_i进行3层WPD，得到\mathit{X}所有数据的各频段小波系数矩阵集合\mathit{W}=\{{\mathit{W}}_{\mathrm{1}},{\mathit{W}}_{\mathrm{2}},\dots ,{\mathit{W}}_m\}.

步骤3 令延迟时间\tau =\mathrm{2}，嵌入维数d=\mathrm{5}，利用\tau和d，对\mathit{W}进行相空间重构，并转换为TRP特征.

步骤4 通过SVD提取TRP特征，获得\mathit{V}=[{\mathit{v}}_{\mathrm{1}},{\mathit{v}}_{\mathrm{2}},\dots ,{\mathit{v}}_m]，作为后续分类的最终特征.

2 基于MPA的SVM参数优化

SVM^［16］分类时需要惩罚因子c和核参数g，具有较强主观性，为此本文引入MPA^［17］来寻优参数，相比粒子群算法（PSO），MPA将优化过程分为3个主要阶段，利用涡流效应，降低算法局部停滞概率，使其不易陷入局部最优.其基本原理描述如下：

1） 在迭代之前，获得猎物的初始位置矩阵{\mathit{P}}_{\mathrm{0}}.

{\mathit{P}}_{\mathrm{0}}={\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}}+{\mathit{R}}_{\mathrm{0}}\cdot \left({\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}}\right).（16）

式中：{\mathit{P}}_{\mathrm{0}}为猎物种群的初始位置矩阵；{\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{a}\mathrm{x}}和{\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}}分别为猎物种群位置矩阵元素取值的上限和下限；{\mathit{R}}_{\mathrm{0}}为［0，1］的随机数矩阵.

取数据集分类准确率的倒数f_t来评价第t次迭代时的适应度值，计算如下：

{\mathit{f}}_t=\frac{n_{\mathrm{L}}}{n_{\mathrm{R}}}.（17）

式中：{\mathit{f}}_t为第t次迭代时的适应度值；n_{\mathrm{L}}为数据集的长度；n_{\mathrm{R}}为分类正确的数据点数.

2） 在迭代初期（，t_max为最大迭代次数），猎物执行布朗运动，捕食者处于初期探索阶段，公式如下：

\left.\begin{array}{l}{\mathit{P}}_{z,t+\mathrm{1}}={\mathit{P}}_{z,t}+w\cdot {\mathit{R}}_{z,t}\otimes {\mathit{s}}_{z,t};\\ {\mathit{s}}_{z,t}={\mathit{R}}_{\mathrm{B},z,t}\otimes ({\mathit{E}}_{z,t}-{\mathit{R}}_{\mathrm{B},z,t}\otimes {\mathit{P}}_{z,t}).\end{array}\right\}（18）

式中：{\mathit{P}}_{z,t+\mathrm{1}}为第t+1次迭代时第z个猎物的位置；z=1，2，…，N，N为猎物总数；{\mathit{P}}_{z,t}和{\mathit{E}}_{z,t}分别为第t次迭代时第z个猎物和捕食者的位置；w为权重变量，控制猎物位置更新的幅度，为固定值，w=0.5；{\mathit{R}}_{z,t}为［0，1］内的随机数向量；\otimes为张量积；{\mathit{s}}_{z,t}为第t次迭代时{\mathit{P}}_{z,t}和{\mathit{E}}_{z,t}将要变化的步长；{\mathit{R}}_{\mathrm{B},z,t}为第t次布朗运动迭代时第i个猎物的随机数向量.

3） 在优化的迭代中期，猎物与捕食者速度相同，猎物种群被平均分为2部分，其中前一半种群执行莱维运动，后一半种群执行布朗运动.

莱维运动中猎物种群个体位置的更新公式如下：

\left.\begin{array}{l}{\mathit{P}}_{z,t+\mathrm{1}}={\mathit{P}}_{z,t}+w\cdot {\mathit{R}}_{z,t}\otimes {\mathit{s}}_{z,t};\\ {\mathit{s}}_{z,t}={\mathit{R}}_{L,z,t}\otimes ({\mathit{E}}_{z,t}-{\mathit{R}}_{\mathrm{L},z,t}\otimes {\mathit{P}}_{z,t}).\end{array}\right\}（19）

式中：{\mathit{R}}_{L,z,t}为t次莱维运动迭代时第z个猎物的随机数向量；z=\mathrm{1,2},\dots ,N/\mathrm{2}.

布朗运动中猎物种群个体位置的更新公式如下：

\left.\begin{array}{l}{\mathit{P}}_{z,t+\mathrm{1}}={\mathit{E}}_{z,t}+{\mathit{P}}_{z,t}\cdot C_t\otimes {\mathit{s}}_{z,t};\\ {\mathit{s}}_{z,t}={\mathit{R}}_{\mathrm{B},z,t}\otimes ({\mathit{R}}_{\mathrm{B},z,t}\otimes {\mathit{E}}_{z,t}-{\mathit{P}}_{z,t});\\ C_t=(\mathrm{1}-t/t_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{2}t/t_{\mathrm{m}\mathrm{a}\mathrm{x}}}.\end{array}\right\}（20）

式中：C_t为第t次迭代时的自适应参数，其中z=N/\mathrm{2},N/\mathrm{2}+\mathrm{2},\cdot \cdot \cdot ,N.

4） 在优化的迭代后期，捕食者和猎物共同执行莱维运动，获得第t+1次迭代后的个体位置{\mathit{P}}_{z,t+\mathrm{1}}，公式如下：

\left.\begin{array}{l}{\mathit{P}}_{z,t+\mathrm{1}}={\mathit{E}}_{z,t}+{\mathit{P}}_{z,t}\cdot C_t\otimes {\mathit{s}}_{z,t};\\ {\mathit{s}}_{z,t}={\mathit{R}}_{\mathrm{L},z,t}\otimes ({\mathit{R}}_{\mathrm{L},z,t}\otimes {\mathit{E}}_{z,t}-{\mathit{P}}_{z,t}).\end{array}\right\}（21）

5） 在第t次迭代中，令猎物矩阵{\mathit{P}}_t=[{\mathit{P}}_{\mathrm{1},t},\dots ,{\mathit{P}}_{z,t},\dots ,{\mathit{P}}_{N,t}{]}^{\mathrm{T}}，根据式（17）获得{\mathit{P}}_t的适应度值列向量{\mathit{f}}_t=[f_{\mathrm{1},t},f_{\mathrm{2},t},\dots ,f_{N,t}{]}^{\mathrm{T}}，并选取{\mathit{f}}_t中的最小值f_{\mathrm{m}\mathrm{i}\mathrm{n},t}对应的{\mathit{P}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}，作为捕食者的位置{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}={\mathit{P}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}，其中，.

6） 当达到最大迭代次数t_{\mathrm{m}\mathrm{a}\mathrm{x}}时，选取{\mathit{f}}_{t_{\mathrm{m}\mathrm{a}\mathrm{x}}}中的最小值f_{\mathrm{m}\mathrm{i}\mathrm{n},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，选取{\mathit{P}}_t中对应位置的元素{\mathit{P}}_{t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，作为捕食者的位置，即{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}={\mathit{P}}_{\mathrm{m}\mathrm{i}\mathrm{n},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，并将t_{\mathrm{m}\mathrm{a}\mathrm{x}}对应的捕食者位置向量{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}=[E_{\mathrm{1},t_{\mathrm{m}\mathrm{a}\mathrm{x}}},E_{\mathrm{2},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}]分别作为SVM的最终结构参数c_{\mathrm{b}}和g_{\mathrm{b}}，即c_{\mathrm{b}}=E_{\mathrm{1},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，g_{\mathrm{b}}=E_{\mathrm{2},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}.

在每次迭代更新猎物位置后，MPA利用基于涡流效应对猎物位置进一步处理.

{\mathit{P}}_{z,t+\mathrm{1}}=\left\{\begin{array}{l}{\mathit{P}}_{t+\mathrm{1}}+G_t\cdot {\mathit{P}}_{\mathrm{0}}\otimes {\mathit{U}}_t,r_t\le p_F;\\ {\mathit{P}}_{t+\mathrm{1}}+[p_F(\mathrm{1}-r_t)+r_t]({\mathit{P}}_{H,t+\mathrm{1}}-{\mathit{P}}_{K,t+\mathrm{1}}),r_t>p_F.\end{array}\right.（22）

式中：{\mathit{P}}_{t+\mathrm{1}}=[{\mathit{P}}_{\mathrm{1},t+\mathrm{1}},{\mathit{P}}_{\mathrm{2},t+\mathrm{1}},\dots ,{\mathit{P}}_{N,t+\mathrm{1}}]为第t+1次迭代采用涡流效应后获得的猎物种群位置矩阵；G_t =（1-t/t_max）{}^{\mathrm{2}t/t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，是第t次迭代时捕食者步长的自适应参数；{\mathit{U}}_t为随机生成二进制向量的二进制数组；r_t为［0，1］内的随机数；p_F为扰动概率因子，本文取0.2；{\mathit{P}}_{H,t+\mathrm{1}}和{\mathit{P}}_{K,t+\mathrm{1}}分别表示第t+1次迭代时，从猎物种群中随机抽取的两只猎物H和K个体的位置向量.

综上所述，基于MPA的SVM参数优化过程的步骤如下：

步骤1 设置海洋捕食者初始参数为d=\mathrm{2}，N=\mathrm{6}，t_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{50}，{\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{16}，{\mathit{P}}_{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}}=-\mathrm{16}，p_F=\mathrm{0.2}；

步骤2 令猎物矩阵{\mathit{P}}_t=[{\mathit{P}}_{\mathrm{1},t},{\mathit{P}}_{\mathrm{2},t},\cdot \cdot \cdot ,{\mathit{P}}_{z,t},\cdot \cdot \cdot ,{\mathit{P}}_{N,t}{]}^{\mathrm{T}}，根据式（17）获得{\mathit{P}}_t的适应度值列向量{\mathit{f}}_t=[f_{\mathrm{1},t},f_{\mathrm{2},t},…,f_{N,t}{]}^{\mathrm{T}}，并选取{\mathit{f}}_t中的最小值f_{\mathrm{m}\mathrm{i}\mathrm{n},t}，选取{\mathit{P}}_t中对应位置元素作为捕食者位置，即{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}={\mathit{P}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}；

步骤3 捕食者在步骤2的每个迭代阶段，利用式（18）~式（22）更新{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t}.若t=t_{\mathrm{m}\mathrm{a}\mathrm{x}}，则将t_{\mathrm{m}\mathrm{a}\mathrm{x}}对应的捕食者位置向量{\mathit{E}}_{\mathrm{m}\mathrm{i}\mathrm{n},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}=[E_{\mathrm{1},t_{\mathrm{m}\mathrm{a}\mathrm{x}}},E_{\mathrm{2},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}]分别作为SVM的最终结构参数c_{\mathrm{b}}和g_{\mathrm{b}}，即c_{\mathrm{b}}=E_{\mathrm{1},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}，g_{\mathrm{b}}=E_{\mathrm{2},t_{\mathrm{m}\mathrm{a}\mathrm{x}}}.

3 基于WP-TRP的故障诊断算法描述

基于WP-TRP的故障诊断过程如图2所示.首先通过递减信息熵准则来确定待诊断数据集 X 的分解层数j，并进行WPD得到各频段的小波系数 W，对各频段的小波系数分别进行TRP转化，获得二维特征图矩阵 R，并通过SVD算法进一步提取其核心特征，获得最终特征 V，并将 V 输入MPA中，得到SVM的最终结构参数c_{\mathrm{b}}和g_{\mathrm{b}}，以此训练SVM，得到最终的故障分类结果.

4 仿真研究

本文实验数据采用凯斯西储大学的标准滚动轴承数据集.所用的滚动轴承数据参数如表2所示，在0.745 7 kW工况、采样频率为12 000 Hz时，分别有正常轴承（N）、外圈故障（OR）、内圈故障（IR）和滚动体故障（B）4种类型，故障的直径分别为0.18，0.36和0.54 mm^［10］.故障数据为1×1 024的行向量，每种故障类型取100条数据，则数据集的长度为\mathrm{100}\times \mathrm{10}=\mathrm{1}\mathrm{000}.所有仿真均在MATLAB 2022b软件上完成，其所用计算机的CPU型号为Intel Core i7-13700H，GPU型号为NVIDIA GeForce RTX 4060，内存为16 GB.

将表2中数据集通过1.5节的步骤进行特征提取，将此特征输入到MPA-SVM中进行模型训练及测试，采用随机抽样的方法，获得训练集和测试集的数据数量之比为7∶3.

对于参数寻优适应度值曲线如图3所示.图3的参数寻优过程以式（17）作为适应度函数，在第25次迭代时，训练准确率就已经达到了100%，此时的训练时间为0.06 s，最优参数c_{\mathrm{b}}为14.452 2，g_{\mathrm{b}}为1.230 5，说明采用本文方法的计算效率较高.

以图3最优参数配置SVM，对测试数据集进行故障诊断并比较性能，如表3所示.

由表3可见，RP-SVD-SVM，TRP-SVD-SVM和RQA-SVM均利用RP或TRP来提取特征，容易导致对故障特征的提取较为单一，不够充分，即便后来与SVD相结合来进一步提取最终特征时，也会因为初始特征提取得不够充分，导致测试正确率较低.而由于TRP蕴含的信息量高于有阈值递归图，因此TRP-SVD-SVM的训练和测试正确率明显高于RP-SVD-SVM.相比而言，TRP-CNN直接将TRP生成的特征与CNN结合，与RQA-SVM和TRP-SVD-SVM相比，由于CNN学习获得的特征较为显著，从而使诊断模型得到的训练和测试正确率较高，但同时，由于采用TRP生成特征时，矩阵的计算量较大，所以计算时间更长.本文提出的WP-TRP，能够利用WPD充分提取初始特征，并利用TRP提高特征的显著性，同时利用SVD实现冗余特征的约简，从而在考虑诊断正确率的前提下，能兼顾计算效率，达到训练和测试正确率均达100.00%，且运行时间最短，为10.66 s，体现了所提方法的有效性.

测试集中，本文方法与传统SVM和CNN方法对每种故障类型的分类正确率比较如表4所示，其中1~10依次对应表2的10种类型数据.

由表4可见，相比其他方法，RP-SVD-SVM对第4类的分类正确率最低，为32.3%；该方法对第9，2，8，7，6类的分类正确率也较低，分别为33.3%，63.6%，65.5%，69.0%，72.0%；此外，TRP-SVD-SVM对第4，9，7类的分类正确率也较低，分别为63.0%，68.0%，78.8%.RQA-SVM对第4类的分类正确率为74.1%.而TRP-CNN和本文提出的WP-TRP与上述3种方法相比精度较高.

这是由于RP-SVD-SVM，TRP-SVD-SVM和RQA-SVM均在时域内利用递归图思想形成初始特征，易于受到数据自身特点和所用时域指标的主观性影响，导致所建的故障诊断模型精度较低；TRP-CNN是将递归图思想与卷积神经网络相结合，其特有的卷积和池化层能够对数据结构更具鲁棒性，更充分地提取数据的非线性特征；WP-TRP方法则利用WPD，充分考虑数据的时频域特征，其频段内的特征提取更加充分，避免了卷积和池化的复杂计算过程，计算效率更高，因此有效性更高.

5 结 语

本文提出了一种基于WP-TRP的故障诊断方法.首先，将本文提出的递减信息熵准则与WPD相结合，获得初始特征.其次将TRP和SVD相结合，获得最终特征.提高了特征的显著性，并提高了计算效率.最后，利用海洋捕食者算法优化支持向量机来进一步提高故障诊断的精度.并通过标准滚动轴承数据集验证所提方法的有效性.

参考文献

[1] 侯东晓， 穆金涛， 方成， 等. 基于GADF与引入迁移学习的ResNet34对变速轴承的故障诊断［J］. 东北大学学报（自然科学版）， 2022， 43（3）： 383-389.

[2] Hou Dong-xiao， Mu Jin-tao， Fang Cheng， et al. Fault diagnosis of variable speed bearings based on GADF and ResNet34 introduced transfer learning［J］. Journal of Northeastern University （Natural Science）， 2022， 43（3）： 383-389.

[3] Yan R Q， Shang Z G， Xu H， et al. Wavelet transform for rotary machine fault diagnosis：10 years revisited［J］. Mechanical Systems and Signal Processing， 2023， 200（5）： 110545.

[4] Sun Y J， Li S H， Wang X H. Bearing fault diagnosis based on EMD and improved Chebyshev distance in SDP image［J］. Measurement， 2021， 176（17）： 109100.

[5] Ye M Y， Yan X A， Chen N， et al. Intelligent fault diagnosis of rolling bearing using variational mode extraction and improved one-dimensional convolutional neural network［J］. Applied Acoustics， 2023， 202： 109143.

[6] Sun Y K， Cao Y， Li P. Fault diagnosis for train plug door using weighted fractional wavelet packet decomposition energy entropy［J］. Accident Analysis & Prevention， 2022， 166： 106549.

[7] 周奇才， 沈鹤鸿， 赵炯， 等. 基于改进堆叠式循环神经网络的轴承故障诊断［J］. 同济大学学报（自然科学版）， 2019， 47（10）： 1500-1507.

[8] Zhou Qi-cai， Shen He-hong， Zhao Jiong， et al. Bearing fault diagnosis based on improved stacked recurrent neural network［J］. Journal of Tongji University （Natural Science）， 2019， 47（10）： 1500-1507.

[9] Liang H P， Cao J， Zhao X Q. Average descent rate singular value decomposition and two-dimensional residual neural network for fault diagnosis of rotating machinery［J］. IEEE Transactions on Instrumentation and Measurement， 2022， 71： 3512616.

[10] 张龙， 胡燕青， 赵丽娟， 等. 采用递归图编码技术与残差网络的滚动轴承故障诊断［J］. 西安交通大学学报， 2023， 57（2）： 110-120.

[11] Zhang Long， Hu Yan-qing， Zhao Li-juan， et al. Fault diagnosis of rolling bearings using recurrence plot coding technique and residual network［J］. Journal of Xi’an Jiaotong University， 2023， 57（2）： 110-120.

[12] 武明泽. 基于小波包分解和FPA-SVM的动车组轴箱轴承故障诊断研究［D］. 北京： 北京交通大学， 2021.

[13] Wu Ming-ze. Research on fault diagnosis of EMU axle box bearing based on wavelet packet decomposition and FPA-SVM［D］. Beijing： Beijing Jiaotong University， 2021.

[14] Han T， Liu C， Yang W G， et al. Deep transfer network with joint distribution adaptation： a new intelligent fault diagnosis framework for industry application［J］. ISA Transactions， 2020， 97： 269-281.

[15] 刘涛， 梁成玉. 信息熵融合的PSO-SVC涡旋压缩机故障诊断［J］. 振动、测试与诊断， 2022， 42（1）： 141-147，200.

[16] Liu Tao， Liang Cheng-yu. PSO-SVC fault diagnosis of scroll compressor based on information entropy fusion［J］. Journal of Vibration， Measurement & Diagnosis， 2022，42（1）： 141-147，200.

[17] Kok T L， Aldrich C， Zabiri H， et al. Application of unthresholded recurrence plots and texture analysis for industrial loops with faulty valves［J］. Soft Computing， 2022， 26（19）： 10477-10492.

[18] 金江涛， 许子非， 李春， 等. 基于深度学习与混沌特征融合的滚动轴承故障诊断［J］. 控制理论与应用， 2022， 39（1）： 109-116.

[19] Jin Jiang-tao， Xu Zi-fei， Li Chun， et al. Rolling bearing fault diagnosis based on deep learning and chaotic feature fusion［J］. Control Theory & Applications， 2022， 39（1）： 109-116.

[20] Cui L L， Liu Y H， Zhao D Z. Adaptive singular value decomposition for bearing fault diagnosis under strong noise interference［J］. Measurement Science and Technology， 2022， 33（9）： 095002.

[21] 郭明军， 李伟光， 杨期江， 等. 深度卷积神经网络在滑动轴承转子轴心轨迹识别中的应用［J］. 振动与冲击， 2021， 40（3）： 233-239， 283.

[22] Guo Ming-jun， Li Wei-guang， Yang Qi-jiang，et al. Application of deep convolution neural network in identification of journal bearing rotor center orbit［J］. Journal of Vibration and Shock， 2021， 40（3）： 233-239， 283.

[23] Wang Z Y， Yao L G， Chen G， et al. Modified multiscale weighted permutation entropy and optimized support vector machine method for rolling bearing fault diagnosis with complex signals［J］. ISA Transactions， 2021， 114： 470-484.

[24] Faramarzi A， Heidarinejad M， Mirjalili S， et al. Marine predators algorithm： a nature-inspired metaheuristic［J］. Expert Systems with Applications， 2020， 152（4）： 113377.

[25] 司莉， 毕贵红， 魏永刚， 等. 基于RQA与SVM的声发射信号检测识别方法［J］. 振动与冲击， 2016， 35（2）： 97-103， 123.

[26] Si Li， Bi Gui-hong， Wei Yong-gang， et al. Detection and identification of acoustic emission signals based on recurrence quantification analysis and support vector machines［J］. Journal of Vibration and Shock， 2016， 35（2）： 97-103， 123.

基金资助

天津市重点研发计划项目(19YFHBQY00040)