基于参数优化VMD-MCKD的强噪声背景下滚动轴承故障特征提取

Feature extraction of rolling bearing faults under strong noise background based on parameter optimized VMD-MCKD

摘要

针对强噪声背景滚动轴承故障特征难以被提取的问题，提出了参数优化变分模态分解（variational mode decomposition，VMD）和最大相关峭度反卷积（maximum correlation kurtosis deconvolution，MCKD）提取滚动轴承故障特征的方法。首先，采用改进麻雀算法对VMD参数进行离线寻优，得到最优参数组合并对原始信号进行分解。其次，根据包络谱峰值因子和样本熵构建出一种新筛选指标，对分解各固有模态函数（intrinsic mode function，IMF）分量进行筛选与重构。然后，对重构信号经改进麻雀算法在线法优化的MCKD进行增强。最后，对增强的信号进行包络解调分析，从而提取滚动轴承故障频率信息。仿真和实验结果表明，该方法能够增强淹没在强噪声中的冲击成分，有效提取滚动轴承故障特征。

Abstract

In order to solve the problem that rolling bearing fault features were difficult to be extracted under strong noise background， parameter optimized variational mode decomposition （VMD） and maximum correlation kurtosis deconvolution （MCKD） were proposed to extract rolling bearing fault features. Firstly， the original signal was decomposed by the optimal combination of parameters obtained by offline optimization of the VMD parameters using the improved sparrow algorithm. Secondly， in order to screen and reconstruct each IMF after decomposition， a new screening metric was constructed based on the envelope spectrum peak factor and sample entropy. Then， the reconstructed signal was augmented with MCKD optimized by the online method of the improved sparrow algorithm. Finally， the bearing failure frequency information was extracted from the enhanced signal by envelope demodulation analysis. Simulation and experimental results show that the proposed method is able to enhance the shock components submerged in the strong noise and effectively extract rolling bearing fault features.

关键词

特征提取 / 滚动轴承 / 变分模态分解 / 最大相关峭度反卷积 / 信号重构

Key words

feature extraction / rolling bearing / variational mode decomposition / maximum correlation kurtosis deconvolution / signal reconstruction

作为旋转机械传动系统中重要支撑零件，滚动轴承一旦出现故障，不仅会造成设备停机，严重时会造成人员受伤。由于滚动轴承的工作环境复杂，各种信号掺杂在一起，故障信号被强噪声所掩盖，导致所采集的信号难以分析^［1］。因此，如何对所采集信号进行有效去噪，并保留故障信息，一直是研究的热点方向。

目前，滚动轴承信号处理主要采用经验模态分解（empirical mode decomposition，EMD）、集合经验模态分解（ensemble empirical mode decomposition，EEMD）和局部均值分解（local mean decomposition，LMD）等方法。Saidi等^［2］将EMD与双光谱相结合，进行轴承故障诊断。Zhao等^［3］将EEMD分解和最小二乘支持向量机相结合进行轴承复合故障诊断。Song等^［4］提出复合多尺度加权排列熵、LMD与支持向量机相结合的滚动轴承故障诊断方法。虽然这些方法在滚动轴承故障诊断领域中取得了一定效果，但上述均属于递归模态分解，存在模态混叠等问题，且缺乏严格的数学理论的支持。针对递归模态分解的缺点，Ma等^［5］提出与极限学习机相结合进行滚动轴承复合故障诊断。VMD属于完全非递归模型，可以有效避免递归模态分解产生的模态混叠等问题。

最小熵解卷积（minimum entropy deconvolution，MED）通过保留信号的冲击成分提取信号的故障特征。张龙等^［6］提出MED与谱峭度算法结合，对故障信号进行二次滤波处理来提取故障特征。为了避免MED仅能突出少数大的尖端脉冲导致的其他冲击丢失问题，Mcdonald等^［7］提出了最大相关峭度反卷积（maximum correlation kurtosis deconvolution，MCKD）算法，MCKD可以突出被噪声淹没的连续冲击脉冲。吕中亮^［8］使用MCKD增强故障信号中的冲击成分，然后利用VMD进行特征提取，最后使用支持向量机进行故障诊断。

VMD算法和MCKD算法需要人为设置参数，而参数取值会影响算法效果。VMD算法需要设置模态分解个数K和惩罚因子α ^［9］；MCKD算法需要提前设置滤波器长度L、反卷积周期T和位移数M ^［10］。因此，为了避免人为选择参数对算法的影响，本文选择改进的麻雀算法对VMD算法中K和α进行离线寻优，对MCKD算法中的L和M进行在线寻优，根据重构信号自适应计算T。通过仿真信号和凯斯西楚大学滚动轴承数据信号对本文方法进行验证。

1 基于参数优化VMD的振动信号分解与重构

1.1 VMD

VMD是一种完全非递归模式的信号分解方法，将原始信号 x \left(t\right)进行Hilbert变换并引入一个指数项 e^{-j{\omega }_{k}t}，计算每个模态 u_k \left(t\right)的信号带宽，将原始信号分解成K个离散的模态函数，其约束变分模型如式（1）所示。

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

式中： \{ u_k \}和 \{ {\omega }_k \}分别为第k（k =1，2，…，K）个模态分量及其中心频率；t为时间； {\partial }_t ( \cdot )为对t求偏导的函数； \delta ( t )为脉冲函数；j为复数单位；“ \mathrm{*}”为卷积符号。

为解决约束变分模型优化问题，引入惩罚因子α和拉格朗日乘子 \lambda，得到公式如式（2）所示。

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

式（2）采用交替方向乘子法连续更新各分量及其中心频率，得到最优解，最终将输入信号分解为K个IMF分量。

1.2 VMD参数优化

VMD信号需要设置合适的模态个数K和惩罚因子α，而K和α的取值对VMD效果影响较大。因此，需对K和α进行参数寻优。寻优过程中最关键的问题是适应度函数的构建，适应度函数用来评价分解的固有模态函数（intrinsic mode function，IMF）分量包含故障特征信息的丰富程度，进而反映了K和α取值的合理性。

包络熵的值反映该分量的概率分布均匀性，互信息的值表示该分量与原始信号的相关联程度。为了反映IMF分量的噪声水平的同时反映出特征信息的完整程度，本文提出了一种新的复合指标，其表达式为

C_{\mathrm{1}} \left(F_i\right) = E_P \left(F_i\right) + \frac{\mathrm{1}}{MI\left(F_i\right)} （3）

式中： E_p ( F_i )为包络熵^［11］，其计算表达式为

E_P \left(F_i\right) = - {\displaystyle \sum _{j=\mathrm{1}}^J} \left(\frac{a_{i,j}}{{\displaystyle \sum _{j=\mathrm{1}}^J}{a}_{i,j}}\right) \mathrm{l} \mathrm{g} \left(\frac{a_{i,j}}{{\displaystyle \sum _{j=\mathrm{1}}^J}{a}_{i,j}}\right) （4）

式中： a_{i,j}为分量IMF _i 经Hilbert解调后得到的包络信号。

MI（F_i ）为第i个F分量与原始信号 \mathit{x} \in {\mathfrak{R}}^J的互信息^［12］，其计算表达式为

M I \left(F_i\right) = \frac{MI\left(F_{\mathit{i}},\mathit{x}\right)}{\sqrt[]{H\left(F_{\mathit{i}}\right)\cdot H\left(\mathit{x}\right)}} （5）

式中： M I \left(F_{\mathit{i}},\mathit{x}\right)为第i个IMF分量与原始信号 x 的互信息；H（ x ）为原始信号的熵；H（F_i ）为F_i 的熵。

本文结合麻雀搜索算法（sparrow search algorithm，SSA）^［13］和改进麻雀搜索算法（sine-cosine and cauchy mutation sparrow search algorithm，SCSSA）^［14］对VMD参数进行寻优，参数优化具体流程如下：

步骤1 读取滚动轴承正常运行状态下的历史数据，确定样本长度J，并将样本记为 x。

步骤2 设置种群大小D、最大迭代次数DC、发现者比例PD、侦察者比例SD、警戒阈值R ₂、安全阈值ST等。设定VMD优化参数的寻优范围 \alpha \in［ {\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}， {\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}］， K \in［ K_{\mathrm{m}\mathrm{i}\mathrm{n}}， K_{\mathrm{m}\mathrm{a}\mathrm{x}}］。因此，种群中麻雀的位置为 X =［K， \alpha］，令时刻t=0。

步骤3 初始化种群，在2维空间内随机生成D个麻雀的初始位置 {\mathit{X}}_d^{\mathrm{0}}。

步骤4 计算种群的适应度函数。

首先，根据第d（d=1，2，…，D）个麻雀位置 {\mathit{X}}_d^t对振动信号进行VMD分解，得到 K_d个IMF分量 \mathrm{I} \mathrm{M} {\mathrm{F}}_{d,k}。

然后，根据式（3）计算每个IMF分量的复合指标 C_{\mathrm{1}} ( F_{d,k} )。

最后，计算第d个麻雀的适应度值 f_d^t，其计算公式为

f_d^t = \mathrm{m} \mathrm{i} \mathrm{n} \left\{C_{\mathrm{1}}(F_{d,k})\right\} \hspace{1em} ( k = \mathrm{1,2} , \cdots , K_d ) （6）

步骤5 确定最优、最差个体位置。

将适应度值 f_d^t按升序排序，选择适应度值前D ₁个麻雀为发现者（D ₁由发现者比例PD确定），其余为跟随者。同时，将最小适应度值记为 f_g^t，其位置为最优个体位置 {\mathit{X}}_{{}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}^t；将最大适应度值记为 f_w^t，其位置为最差个体位置 {\mathit{X}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}^t。

步骤6 位置更新。

发现者位置更新公式为

（7）

式中：运算符“ \cdot ”为两个长度相等向量对应元素相乘； {\mathit{X}}_d^t为t时刻第d个麻雀的位置； {\mathit{X}}_{{}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}^t为当前最优位置； {\mathit{r}}_{\mathrm{2}} \in [ \mathrm{0,2} \mathrm{\pi } ]是决定麻雀的移动距离随机数；随机向量 {\mathit{r}}_{\mathrm{3}} \in [ \mathrm{0,2} \mathrm{\pi } ]控制最优个体对麻雀后一位置的影响；随机数 R_{\mathrm{2}} \in［0，1］表示警戒值； S T \in［0.5，1］为安全值。

r_{\mathrm{1}}^{\text{'}}为非线性正余弦算法的递减搜索因子，其计算公式为

{r_{\mathrm{1}}}^{\mathrm{\text{'}}} = {\left[\mathrm{1}-{\left(\frac{t}{DC}\right)}^{\eta }\right]}^{\mathrm{1}/\mathrm{\eta }} （8）

式中： \eta ( \eta \ge \mathrm{1} )为调节系数。

\omega为非线性权重因子，其计算公式为

\omega = \frac{{\mathrm{e}}^{\frac{t}{DC}}-\mathrm{1}}{\mathrm{e}-\mathrm{1}} （9）

引入柯西变异策略后，跟随者位置更新公式为

\begin{array}{l}{\mathit{X}}_d^{t+\mathrm{1}}={\mathit{X}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t+\mathrm{c}\mathrm{a}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{y}\left(\mathrm{0,1}\right){\mathit{X}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t\\ (d=D_{\mathrm{1}}+\mathrm{1},D_{\mathrm{1}}+\mathrm{2},\cdots ,D)\end{array} （10）

式中：cauchy（0，1）为标准柯西分布函数。

根据步骤4计算更新后位置 {\mathit{X}}_d^{t+\mathrm{1}}的适应度值为 f_d^{t+\mathrm{1}}。

根据侦察者比例SD，从更新后的种群中随机选取个体作为侦察者，其位置更新公式为

\begin{array}{l}{\mathit{X}}_i^{t+\mathrm{1}}=\left\{\begin{array}{c}{\mathit{X}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t+\beta \left|{\mathit{X}}_i^{t+\mathrm{1}}-{\mathit{X}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t\right|,\hspace{1em} f_i^{t+\mathrm{1}}>f_{\mathrm{1},g}^t\\ {\mathit{X}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t+\kappa \left|\frac{{\mathit{X}}_i^{t+\mathrm{1}}-{\mathit{X}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}^t}{\left(f_g^t-f_w^t\right)+\epsilon }\right|,f_i^{t+\mathrm{1}}=f_{\mathrm{1},g}^t\end{array}\right.\\ i=\mathrm{1,2},\cdots ,SD\times D\end{array} （11）

式中： \beta为步长修正系数，服从标准正态分布； \kappa \in [ \mathrm{0,1} ]为随机数； \epsilon为极小常数。

根据步骤4计算更新后侦察者的适应度，并更新全局适应度值。

步骤7 判断是否达到最大迭代次数DC。若条件成立，则输出最佳位置 {\mathit{X}}_{{}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}}^{DC}，即最优参数组合（K ^*，α^* ）。否则，令t=t+1，并返回步骤5。

1.3 信号重构

经过1.2节的参数寻优，得到最优的参数组合［K ^*，α^* ］。对原始振动信号进行VMD分解，得到K ^* 个F分量。通过筛选指标对K ^* 个F分量进行筛选，筛选指标公式如式（12）所示。

E S_k = \frac{Ec\left(F_k\right)}{\mathrm{S}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{E}\mathrm{n}\left(m,r,F_k\right)} （12）

式中：分子为第k个模态分量F_k 的包络谱最大峰值因子^［15］，其公式为

E c \left(F_k\right) = \frac{\mathrm{m}\mathrm{a}\mathrm{x}\left(c_k\left(q\right)\right)}{\sqrt[]{{\displaystyle \sum _{q=\mathrm{1}}^Q}{c}_k{\left(q\right)}^{\mathrm{2}}/Q}} （13）

式中： c_k \left(q\right) \left(q=\mathrm{1,2},\cdots ,Q\right)为分量F_k 的包络谱幅值序列。

式（12）中，分母为第k个模态分量F_k 的样本熵^［16］，其公式为

\mathrm{S} \mathrm{a} \mathrm{m} \mathrm{p} \mathrm{l} \mathrm{e} \mathrm{E} \mathrm{n} \left(m,r,F_k\right) = - \mathrm{l} \mathrm{n} \frac{A^{m+\mathrm{1}}\left(r\right)}{B^m\left(r\right)} （14）

式中： B^m ( r )、 A^{m+\mathrm{1}} ( r )分别是由 F_k构造的m、m+1维时间序列在相似容限r下的自相似概率。

该筛选指标的构建既考虑了F分量中的周期性冲击，又衡量了信号的有序性。

选择 K^{\mathrm{*}}个分量筛选指标的平均值 \overline{ES}作为阈值，选择满足 E S_k \ge \overline{ES}的分量进行信号重构，则重构信号 \tilde{\mathit{y}}为

\tilde{\mathit{y}} = \sum _{ES(\mathrm{I}\mathrm{M}{\mathrm{F}}_k)\ge \overline{ES}} F_k （15）

2 基于参数优化MCKD的信号处理

2.1 最大相关峭度解卷积

MCKD通过寻找最佳滤波器对包含噪声的振动信号进行滤波处理，使滤波后的振动信号的相关峭度最大，从而突出冲击成分，使振动信号特征更加明显。输入信号为 \tilde{\mathit{y}} = [ {\tilde{y}}_{\mathrm{1}} , {\tilde{y}}_{\mathrm{2}} , \cdots , {\tilde{y}}_N ]；滤波后的周期性冲击信号 \mathit{y} = [ y_{\mathrm{1}} , y_{\mathrm{2}} , \cdots , y_N ]。

寻找最优滤波器参数，使滤波后的信号具有最大相关峭度，其相关峭度的公式为

C K_M ( T ) = \frac{{\displaystyle \sum _{n=\mathrm{1}}^N}({\displaystyle \prod _{m=\mathrm{0}}^M}{y}_{n-mT}{)}^{\mathrm{2}}}{{({\displaystyle \sum _{n=\mathrm{1}}^N}{y}_n^{\mathrm{2}})}^{M+\mathrm{1}}} （16）

式中：T为反卷积的周期；M为位移数； y_{n-mT}为 y_n在不同m值下滤波器输出信号的延迟。

对式（16）求导，可得到滤波器系数 {\mathit{f}}^s（s=1，2，…，S），其迭代公式为

\left\{\begin{array}{l}{\mathit{y}}^{s-\mathrm{1}}={\displaystyle \sum _{l=\mathrm{1}}^L}{f}_l^{s-\mathrm{1}}{\tilde{y}}_{n-l+\mathrm{1}}\\ {\mathit{\beta }}^s=[y_{\mathrm{1}}^{s-\mathrm{1}}{y}_{{}_{\mathrm{1}-T}}^{s-\mathrm{1}}\cdots y_{{}_{\mathrm{1}-MT}}^{s-\mathrm{1}},\cdots ,y_{{}_N}^{s-\mathrm{1}}{y}_{N-T}^{s-\mathrm{1}}\cdots y_{N-MT}^{s-\mathrm{1}}{]}^{\mathrm{T}}\\ {\mathit{A}}_m^s=[(y_{\mathrm{1}}^{s-\mathrm{1}}{y}_{\mathrm{1}-T}^{s-\mathrm{1}}\cdots y_{\mathrm{1}-MT}^{s-\mathrm{1}}{)}^{\mathrm{2}}/y_{{}_{\mathrm{1}-mT}}^{s-\mathrm{1}},\cdots ,\\ (y_N^{s-\mathrm{1}}{y}_{N-T}^{s-\mathrm{1}}\cdots y_{N-MT}^{s-\mathrm{1}}{)}^{\mathrm{2}}/y_{N-mT}^{s-\mathrm{1}}{]}^{\mathrm{T}}\\ {\mathit{f}}^s=\frac{{‖{\mathit{y}}^{s-\mathrm{1}}‖}^{\mathrm{2}}}{\mathrm{2}{‖{\mathit{\beta }}^s‖}^{\mathrm{2}}}({\mathit{Y}}_{\mathrm{0}}{\mathit{Y}}_{\mathrm{0}}^{\mathrm{T}}{)}^{-\mathrm{1}}{\displaystyle \sum _{m=\mathrm{0}}^M}{\mathit{Y}}_{mT}{\mathit{A}}_m^s\end{array}\right. （17）

式中：f =（f ₁，f ₂， \cdots , f_{\mathrm{L}} )为滤波器系数；L为滤波器长度；向量 {\mathit{\beta }}^s和 {\mathit{A}}_m^s均由s-1次迭代得到的滤波后的信号 {\mathit{y}}^{s-\mathrm{1}}得到； {\mathit{Y}}_{mT}矩阵第一行为输入信号向右偏移mT个单位，原始信号的第一位偏移到该行的第mT+1位，第N-mT位为该行最后一位，前mT位用0占位；矩阵的第L行为第L-1行向右偏移1个单位的原始信号，该矩阵的表达式为

{\mathit{Y}}_{mT} = {\left[\begin{array}{cccc}{\tilde{y}}_{\mathrm{1}-mT}& {\tilde{y}}_{\mathrm{2}-mT}& \cdots & {\tilde{y}}_{N-mT}\\ \mathrm{0}& {\tilde{y}}_{\mathrm{1}-mT}& \cdots & {\tilde{y}}_{N-mT-\mathrm{1}}\\ ⋮& ⋮& & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & {\tilde{y}}_{N-mT-L+\mathrm{1}}\end{array}\right]}_{L\times N}

2.2 MCKD参数优化

MCKD算法中滤波器长度L、反卷积周期T和位移数M的选择对处理效果有较大的影响。其中，T是采样频率f_s 与重构信号的主频f之商的整数部分（四舍五入）^［17］。因此，仅需对T和M进行参数寻优。

寻优算法中最重要的是指标函数的选择，本节选择计算滤波后信号 y 的包络谱最大峰值因子的倒数为指标函数，其公式为

C_{\mathrm{2}}\left(\mathit{y}\right)=\mathrm{1}/Ec(\mathit{y}) = \frac{\sqrt[]{{\displaystyle \sum _{z=\mathrm{1}}^Z}b{(z)}^{\mathrm{2}}/Z}}{\mathrm{m}\mathrm{a}\mathrm{x}\left(b(z)\right)} （18）

式中： b \left(z\right) \left(z=\mathrm{1,2},\cdots ,Z\right)为滤波后信号 y 包络谱幅值序列。C ₂的值越小，周期性冲击越强，故障特征越明显。

基于改进麻雀算法对MCKD参数优化的简化流程如下：

步骤1 读取1.3节重构信号 \tilde{\mathit{y}}；设置寻优范围 L \in［ L_{\mathrm{m}\mathrm{i}\mathrm{n}}， L_{\mathrm{m}\mathrm{a}\mathrm{x}}］、 M \in［ M_{\mathrm{m}\mathrm{i}\mathrm{n}}， M_{\mathrm{m}\mathrm{a}\mathrm{x}}］和种群大小D、DC、PD、SD、滤波器迭代次数S等参数和种群位置 \mathit{X} = \left[L,M\right]。

步骤2 初始化滤波器系数 {\mathit{f}}_d^{\mathrm{0}}（d=1，2，…，D）并计算 {\mathit{Y}}_{\mathrm{0},d}和 ( {\mathit{Y}}_{\mathrm{0},d} {\mathit{Y}}_{\mathrm{0},d}^{\mathrm{T}} {)}^{-\mathrm{1}}，进而得到滤波后信号 {\mathit{y}}_d^{\mathrm{0}}。

步骤3 根据 {\mathit{y}}_d^{s-\mathrm{1}}（s=1，2，…，S）计算 {\mathit{A}}_m^s和 {\mathit{\beta }}^s。

步骤4 根据式（17）更新滤波器系数 {\mathit{f}}_d^s。

步骤5 重复步骤3—5，直至达到最大迭代次数S，输出相关峭度最大的滤波信号 {\mathit{y}}_d^S。

步骤6 计算 {\mathit{y}}_d^S的适应度值C _{2， d} 并进行排序，确定发现者、跟随者和侦查者。

步骤7 根据位置更新公式更新发现者、跟随者和警戒者三类麻雀的位置，并更新全局适应度及最优位置。

步骤8 判断是否达到最大迭代次数DC，若是，则输出最佳参数组合（L^*，M ^* ），否则重复步骤2—8。

3 故障特征提取与诊断流程

基于参数优化VMD与MCKD相结合的轴承故障特征提取与诊断的具体步骤如下：

步骤1 读取待分析滚动轴承振动信号。

步骤2 根据1.2节离线寻优得到最优参数（K^*，α^* ），对该信号进行VMD分解。

步骤3 根据式（12）计算K ^* 个分量的筛选指标 E S_k（k=1，2，…，K ^* ）和阈值 \overline{ES}，并选择 E S_k大于阈值的IMF分量，根据式（15）进行信号重构。

步骤4 自适应计算T并按照2.2节对MCKD参数（L，M）进行在线寻优，再根据最优参数（L^*，M^* ）对重构信号 \tilde{\mathit{y}}进行滤波。

步骤5 对滤波后的信号 y 进行包络解调分析，确定其特征频率（主频）和轴承故障类型。

4 仿真信号分析

为验证本文所提方法正确性，构造轴承内圈故障仿真信号如式（19）所示。

\left\{\begin{array}{l}x\left(t\right)=\sum _k{A}_{k}h\left(t-k{T}_{\mathrm{1}}-{\tau }_k\right)+n\left(t\right)\\ A_k=A_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{50}\mathrm{\pi }t\right)+\mathrm{1}\\ h\left(t\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-Ct\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }{f}_{n}t\right)\end{array}\right. （19）

式中：幅值 A_{\mathrm{0}}为0.2；衰减系数C为800；共振频率 f_n为4 000 Hz； {\tau }_k为第 k次冲击相对周期T的微小波动，随机波动服从正态分布；n（t）为高斯白噪声，染噪信号的信噪比设置为-12 dB。采样频率为12 800 Hz；分析点数为4 096。 T_{\mathrm{1}} = \mathrm{1} / f_f为内圈故障周期，f_f =162 Hz为故障频率。

通过对VMD参数寻优，得到最优参数K=8、α=1 698，其适应度函数迭代曲线如图1所示。本文方法迭代在第4代收敛，使用文献［18］方法迭代，则在第24代收敛。因此，本文方法收敛相对较快。

根据最优参数对仿真信号进行VMD分解后，仿真信号各IMF分量的ES值如表1所示，其筛选指标的平均值为6.24。因此，对于该样本信号选择第5和第6个IMF分量进行信号重构。

为了验证本文信号重构方法的有效性，分别对加噪信号和重构后的信号进行包络谱分析，其包络谱对比图如图2所示。图2a中频率成分杂乱，图2b中故障特征信息162.5 Hz（主频）突出，但没有倍频现象，无法确保诊断有效。因此，需对重构信号进一步分析。

重构信号按在线优化参数（L ^*=174，M ^*=1）进行MCKD滤波，滤波后的信号包络谱如图3所示。加噪信号直接进行MCKD处理的包络谱如图4所示。

由图3和图4对比可知，原始信号经本文方法处理后的包络谱中的主频与倍频均清晰可见。然而，原始信号直接进行MCKD处理的包络谱中突出频率成分仅有主频（162.5 Hz），没有倍频现象，不能确保诊断有效。因此，证明了在MCKD处理前进行VMD分解的必要性。同时，说明了本文提出方法可以有效提取强噪声背景下的故障频率特征信息。

5 实验信号分析

为了进一步验证本文所提出方法的有效性，本文采用美国凯斯西楚大学的实验数据^［19］进行滚动轴承故障信号的分析，其滚动轴承试验装置如图5所示。

本文选用滚动轴承内圈早期故障信号（数据标签105DE、故障频率162.2 Hz）、外圈早期故障信号（数据标签130DE、故障频率107.4 Hz）和滚动体早期故障信号（数据标签291FE、故障频率117.8 Hz）对本文提出的算法进行验证。其中，样本点数为4 096。为了验证本文方法可有效提取故障特征，在故障信号中添加合适的噪声将故障特征完全掩盖。

本文采用正常轴承数据对VMD参数进行寻优，得到最优参数K ^*= 6、α ^*=1 417。

5.1 内圈故障特征提取

为了验证本文在强噪声背景下特征提取的有效性，在内圈故障信号中添加了信噪比为-8dB的噪声，将染噪信号作为原始信号，其时域波形如图6所示。

根据最优参数对原始信号进行VMD分解。内圈信号各IMF分量的ES值如表2所示，其筛选指标的平均值为5.45。因此，对该样本信号选择第2和第4个IMF分量进行信号重构。

为了验证本文方法的合理性，分别对加噪原始信号、重构信号及按在线优化参数（L ^*=435、M ^* =7）进行MCKD滤波后的信号进行包络谱分析，其包络谱如图7所示。图7a中无法找出故障特征频率，图7b中仅主频（161.1 Hz）突出，而经本文方法分析得到的最终结果即图7c，主频及其倍频成分均清晰可见。与内圈故障特征频率162.2 Hz比较后，可判断为轴承内圈故障。因此，证明了该方法的合理性。

5.2 外圈故障特征提取

为了验证本文方法的有效性，在外圈故障信号中添加了信噪比为-12 dB的噪声，将染噪信号作为原始信号。图8为外圈故障信号包络谱对比。图8a中没有任何突出频率成分，图8b中的主频108.4 Hz突出，最后经本文方法获得的图8c中故障特征频率及其倍频现象均清晰可见。与外圈故障特征频率107.4 Hz比较后，可判断为轴承外圈故障。因此，证明了该方法有效提取了强噪声背景下的轴承外圈故障特征。

5.3 滚动体故障特征提取

为了验证本文方法的鲁棒性，在滚动体故障信号中添加信噪比为-4 dB的噪声。对滚动体故障的三阶段信号进行包络谱分析，其包络谱对比如图9所示。

原始信号经本文方法处理后的包络谱中主频（117.2 Hz）及其二倍频（234.4 Hz）、三倍频（351.6 Hz）和四倍频（468.8 Hz）均清晰可见，而文献［20］出现的最大倍频仅为3倍频。滚动体故障特征频率为117.8 Hz，可判断为轴承滚动体故障。因此，本文方法也适用于被强噪声掩盖的滚动体故障特征提取。

6 结论

在强噪声背景下，仅采用优化的VMD或MCKD算法均不能准确判断滚动轴承故障类型。基于本文提出的新筛选指标，进行VMD参数寻优并重构，有效抑制了噪声成分，同时增强了信号的瞬态冲击特性。经参数在线优化的MCKD算法能够增强重构信号中的冲击成分，使倍频现象突出。因此，经改进VMD联合MCKD方法能够准确提取滚动轴承被强噪声掩盖的故障频率特征，进而诊断出滚动轴承故障类型。仿真信号和实验信号均表明本文方法的合理性、有效性和鲁棒性。

参考文献

[1] 梁子寅，刘向东，卢宗慧，滚动轴承故障诊断分析方法研究综述［J］.制造业自动化，2023，45（5）：172-175.

[2] Saidi L， Ben A J， Fnaiech F.Bi-spectrum based-EMD applied to the non-stationary vibration signals for bearing faults diagnosis［J］.ISA Transactions，2014，53（5）：1650-1660.

[3] Zhao Y X， Fan Y， Li H，et al.Rolling bearing composite fault diagnosis method based on EEMD fusion feature［J］.Journal of Mechanical Science and Technology，2022，36（9）：4563-4570.

[4] Song E Z， Gao F， Yao C，et al.Research on rolling bearing fault diagnosis method based on improved LMD and CMWPE［J］.Journal of Failure Analysis and Prevention，2021，21（5）：1714-1728.

[5] Ma J， Yu S， Cheng W.Composite fault diagnosis of rolling bearing based on chaotic honey badger algorithm optimizing VMD and ELM［J］.Machines，2022，10（6）：469.

[6] 张龙，胡俊锋，熊国良.基于MED和SK的滚动轴承循环冲击特征增强［J］.振动、测试与诊断，2017，37（1）：97-101.

[7] Mcdonald G L， Zhao Q， Zuo M J.Maximum correlated kurtosis deconvolution and application on gear tooth chip fault detection［J］.Mechanical Systems and Signal Processing，2012，33：237-255.

[8] 吕中亮.基于变分模态分解与优化多核支持向量机的旋转机械早期故障诊断方法研究［D］.重庆：重庆大学，2016.

[9] Zhou Z Y， Chen W H， Yang C.Adaptive range selection for parameter optimization of VMD algorithm in rolling bearing fault diagnosis under strong background noise［J］.Journal of Mechanical Science and Technology，2023，37（11）：5759-5773.

[10] Tang G J， Wang X L， He Y L.Diagnosis of compound faults of rolling bearings through adaptive maximum correlated kurtosis deconvolution［J］.Journal of Mechanical Science and Technology，2016，30（1）：43-54.

[11] 唐贵基，王晓龙.参数优化变分模态分解方法在滚动轴承早期故障诊断中的应用［J］.西安交通大学学报，2015，49（5）：73-81.

[12] 鲁铁定，何锦亮，贺小星，参数优化变分模态分解的GNSS坐标时间序列降噪方法［J］.武汉大学学报（信息科学版），2024，49（10）：1856-1866.

[13] Xue J K， Shen B.A novel swarm intelligence optimization approach：sparrow search algorithm［J］.Systems Science & Control Engineering，2020，8（1）：22-34.

[14] 李爱莲，全凌翔，崔桂梅，融合正余弦和柯西变异的麻雀搜索算法［J］.计算机工程与应用，2022，58（3）：91-99.

[15] 张俊，张建群，钟敏，基于PSO-VMD-MCKD方法的风机轴承微弱故障诊断［J］.振动、测试与诊断，2020，40（2）：287-296.

[16] 朱俊杰，张清华，朱冠华，基于EEMD与GWO-SVM的石化机组轴承故障诊断［J］.自动化与仪表，2023，38（11）：60-65.

[17] 任良，甄龙信，赵云，基于SSA-VMD-MCKD的强背景噪声环境下滚动轴承故障诊断［J］.振动与冲击，2023，42（3）：217-226.

[18] 刘迎松，魏志刚，束海星，基于参数自适应VMD和MCKD的滚动轴承微弱故障特征提取［J］.噪声与振动控制，2023，43（3）：102-109.

[19] Smith W A， Randall R B.Rolling element bearing diagnostics using the case western reserve university data：a benchmark study［J］.Mechanical Systems and Signal Processing，2015，64：100-131.

[20] Zhang J， Zhang J Q， Zhong M，et al.Detection for incipient damages of wind turbine rolling bearing based on VMD-AMCKD method［J］.IEEE Access，2019，7：67944-67959.

基金资助

国家自然科学基金(62003223)