基于优化变分模态分解与混沌分形融合的滚动轴承故障识别

2022-10-19金江涛许子非

动力工程学报 2022年10期

孙康，金江涛，李春,2，许子非

(1.上海理工大学能源与动力工程学院，上海 200093；2.上海市动力工程多相流动与传热重点实验室，上海 200093)

风能作为应用最为广泛的可再生能源之一，以其稳定、丰富及易于开发等优势，对碳中和和碳达峰起到关键作用[1]。轴承作为应用最广泛的机械传动零部件之一，在高性能舰船推进电机、先进航空发动机及大型风力发电机组等精密旋转机械中具有至关重要的作用，其健康状况直接影响设备性能[2]。滚动轴承是机械传动系统的重要部件之一，能否正常运转关系到整台设备的性能和使用寿命[3]。但恶劣的工作环境导致轴承随服役时间增加出现各种故障。为保证设备正常运作，需要实时监测轴承工作时产生的振动信号，确保维护人员对设备的状态作出判断并对设备进行及时维修[3]。

轴承振动信号呈现出非线性与非平稳性特征[3]，传统频谱分析法存在故障判别方式与准确性不足等问题，无法满足实际工程需求。

向玲等[4]基于小波分解提取故障轴承特征信息。隋文涛等[5]采用经验模态分解(EMD)方法获取轴承故障特征频率。李洋等[6]利用辅助经验模态分解(EEMD)方法对信号降噪，完成故障识别。但上述方法都存在缺乏自适应性、端点效应及模态混叠等问题，无法处理复杂非平稳信号，为此Dragomiretskiy等[7]提出一种变尺度信号处理方法——变分模态分解(VMD)方法，该方法避免了上述问题，具有良好的鲁棒性，被国内外学者广泛应用于故障诊断领域。

但VMD方法需人为设置惩罚因子c和模态分解数K。c取值过小，易造成分解所得本征模态函数(IMF)带宽过大；反之过小。若K取值过大，原信号出现过分解现象；反之原信号无法完全分离，存在模态混叠问题[8]。唐贵基等[9]根据包络谱特征因子选择VMD参数，但未考虑转频对包络信号峰值的影响。张俊等[10]利用粒子群算法优化VMD参数，虽避免了人为经验影响，但粒子群算法易陷入局部最优，误差较大。

当结构复杂的设备出现故障时，从统计学角度分析，其振动信号具有自相似性特征，故可用混沌分形理论来描述信号的复杂性和不规则程度，揭示隐藏于随机现象中的真实规律[11]。

笔者以轴承振动信号为研究对象，将Lyapunov指数λ视为特征值。λ>0，表明信号具有混沌特性，且混沌特性与λ线性相关，以λ最大值为迭代寻优条件，选取VMD分解参数，提出优化变分模态分解(OLVMD)方法，以分解轴承原始振动信号。不同于传统故障诊断方法，OLVMD方法完全脱离频域分析，诊断时间更迅速，形式更简洁，准确率更高，具有独特优势。原信号经OLVMD方法分解得到若干IMF分量，但故障特征仅存于少数IMF分量中，大多为噪声信号。考虑到混沌具有分维性、敏感性及随机性等特点,可量化分析振动信号[11]，将各故障Lyapunov指数最大的前3个分量作为含故障信息的特征信号，重构组成新的振动信号。为准确进行故障识别，笔者引入分形理论，计算轴承在不同状态下的混沌关联维数，并采用拟合偏差平方和方法对传统的关联维数计算方法进行改进。从分形角度分析轴承振动信号，以其关联维数为特征量，定量描述轴承工作状态的变化，并分析了损伤轴承的实验数据。

1 变分模态分解

VMD将待分析信号视为若干模态分量线性叠加而成。各IMF分量被定义为调幅调频信号，表达式如下：

uk(t)=Ak(t)cosφ(t)

(1)

式中：t为时间；uk(t)为第k个IMF；Ak(t)为瞬时幅值；φ(t)为信号相位。

通过希尔伯特变换计算各IMF分量单边谱,再估计中心频率，构造变分模型：

(2)

式中：k=1,2,…,K；σ(t)为单位脉冲函数；ωk为中心频率；f(t)为输入信号；⊗表示卷积运算；∂t表示偏导运算；j为虚数。

引入惩罚因子及Lagrange乘法算子将约束变分问题变为非约束变分问题，得到扩展Lagrange表达式。

求解过程中各分量中心频率与带宽不断更新，直至满足迭代停止条件：

(3)

迭代终止时，信号频域特性已完成自适应分离,然后通过傅里叶逆变换转换到时域。

2 算法改进

2.1 优化变分模态分解方法

基于VMD方法分解振动信号时，惩罚因子c和模态分解数K直接影响分解精度[8]。Lyapunov指数λ可定量反映非线性振动信号的特征信息，λ越大，故障特征越显著[11]。故本文以λ为特征值，选取最佳参数组合[c,K]，以分解所得IMF分量中λ最大的IMF分量为特征信号。基于最大Lyapunov指数提出OLVMD方法，流程如图1所示。

图1 OLVMD方法流程图

2.2 有效性验证

为验证OLVMD方法可有效降噪并提取信号中的特征信息，以一段噪声序列为研究对象，分别采用OLVMD方法和VMD方法对其进行预处理，以分解所得IMF分量中λ最大的IMF分量为特征信号，其时域图如图2所示。其中OLVMD方法最佳分解参数组合为c=1 200，K=10；借鉴刘洋等[12]所提方法，根据EEMD方法递归特性确定模态分解数K，c与采样频率相同，即c=1 000，K=3。

图2 噪声信号处理前后对比

由图2可知，采用OLVMD方法处理后的信号降噪效果显著，为定量表示，计算2种方法处理后信号的信噪比(SNR，SNR)，结果如表1所示。SNR越高，降噪效果越好[13]。由表1可知，原信号含有大量噪声，信号特征被淹没，信噪比最高。相比VMD方法，经OLVMD方法处理后的信号信噪比较低，降噪效果更明显。OLVMD方法在处理非线性信号时，更易获取纯净故障信息。

表1 有效性验证

2.3 关联维数计算

根据文献[13]可得到关联维数D(m,r)：

(4)

式中：m为嵌入维数；Cm(r)为关联积分；r为超球半径;r1

为此，采用拟合偏差平方和方法选取线性区域，以提高计算精确度，将lnCm(r)-lnr曲线中间部分视为直线，两侧为曲线。将整段曲线横坐标记为[d1,d2]，其中直线段为[lnr1,lnr2]，满足d1≤lnr1≤lnr2≤d2。具体拟合方法如下：

(1)d1≤lnr≤lnr1时，以多项式曲线y=f1(x)进行拟合，拟合偏差平方和记为ε1：

(5)

(2) lnr1

(6)

(3) lnr2≤lnr≤d2时，以直线y=f2(x)进行拟合，拟合偏差平方和记为ε3：

(7)

拟合偏差平方总和(即误差)ε=ε1+ε2+ε3，当其值达到最小值时，[lnr1，lnr2]即为最佳线性区域。采用多项式曲线对lnCm(r)-lnr曲线进行拟合时，一般幂次取2较为合适。

2.3.1 仿真信号验证

以经典混沌动力系统Lorenz吸引子为例，采用改进G-P算法求取关联维数，其包含3个状态变量的常微分方程：

(8)

式中：Pr为普朗特数；Ra为规范化的瑞利数；b的大小与吸引子轨迹形状有关。

基于四阶龙格库塔法对上式求解，其中Pr=10，b=3，Ra=25，起始点(x,y,z)取(0,2,9)，迭代步长h为0.01，得到一个5 000×3的时间序列，选择X变量进行分析，Lorenz吸引子轨迹图及X变量时间序列如图3所示。

(a) Lorenz吸引子轨迹图

采用互信息函数法计算得到延迟时间为10，嵌入维数m从上至下依次为2,4,6，…，30所对应的双对数曲线lnCm(r)-lnr如图4所示。

图4 ln Cm(r)-ln r曲线(无噪声)

由图4可知，尺度较小时，曲线分散且不规则；随着尺度增大，曲线与lnr轴逐渐平行，并向同一方向收缩；处于中间尺度时，曲线接近于斜直线，可描述系统分形特征，分形维数在此区间具有意义。

基于拟合偏差平方和方法对图4中线性较好的部分采用最小二乘法进行拟合，计算各嵌入维数对应的曲线斜率，以最佳嵌入维数对应的曲线斜率为该系统的关联维数D(GP)，结果如图5所示。

由图5可知，嵌入维数较小时，嵌入维数与关联维数呈正相关，当嵌入维数大于16时，关联维数趋于稳定，不再增长。因此，最佳嵌入维数为16，其对应的曲线斜率3.231 0为Lorenz吸引子关联维数。

图5 无噪声时关联维数与嵌入维数的关系

2.3.2 抗噪性能分析

考虑到外部环境影响，实际轴承振动信号难免受到噪声干扰，为验证关联维数具有强鲁棒性。在Lorenz吸引子中加入不同大小的高斯白噪声，计算系统关联维数。为准确揭示关联维数变化规律，规定若其波动范围>0.01，则表明关联维数受噪声干扰，分析过程如下：

(1) 不断减小信噪比，当SNR>6 dB时，Lorenz系统嵌入维数值稳定在16，关联维数波动极小，可忽略不计，表明关联维数不受噪声干扰。

(2) 1 dB≤SNR≤6 dB时，系统嵌入维数和关联维数发生细微改变。以SNR=4 dB为例，其lnCm(r)-lnr曲线与关联维数如图6所示。

(a) ln Cm(r)-ln r曲线

(b) D(GP)-m曲线

计算不同信噪比对应的最佳嵌入维数与关联维数，结果如表2所示。

表2 不同信噪比时最佳嵌入维数与关联维数的变化

由表2可知，随着信噪比增大，关联维数不断减小，当信噪比增大到一定值时，关联维数不再变化，表明关联维数具有一定的抗噪性。

(3) 当SNR<1 dB时，lnCm(r)-lnr曲线斜率不再收敛于一个稳定值。以SNR=-4 dB为例，双对数曲线与关联维数如图7所示。

(a) ln Cm(r)-ln r曲线

由图7可知，曲线斜率不再收敛，无法判定最佳嵌入维数及其对应的关联维数。

通过上述分析可知，关联维数具有一定抗噪性，且其存在2个临界值Q1和Q2，其中Q1≤Q2。当SNRQ2时，振动信号关联维数不会产生波动，信号不受噪声影响，关联维数可准确表达信号的状态信息。

3 结果与分析

3.1 实验数据集

利用西安交通大学轴承数据中心提供的实验数据[14]验证所提分形故障诊断方法的有效性及准确性。实验平台如图8所示，该装置以交流电动机为动力源，通过转速控制器改变轴承转速，加速度传感器位于轴承的水平及竖直方向。

图8 轴承实验平台装置图

所研究的滚动轴承相关参数如表3[14]所示。

表3 轴承参数

设置采样频率为25.6 kHz，对转速为2 100 r/min和2 250 r/min 2种状态4种工况下的轴承振动数据进行采集。具体参数如表4所示。

表4 轴承加速寿命试验工况

3.2 特征选取

原信号在4种工况下的时域及频域图见图9，其中f为频率。观察图9可知，各故障振动信号时域图和频域图虽存在一定差异，但只能定性分辨几种故障的不同，无法直接判断故障类型。为提取有效故障特征，基于相空间重构法绘制4种故障三维吸引子轨迹，分析各振动信号的混沌特性，采用Cao方法与互信息函数法计算最佳延迟时间和嵌入维数，结果如表5所示，吸引子轨迹如图10所示。

表5 最佳延迟时间和嵌入维数

由图10可知，4种故障相图均呈冗杂毛球状，吸引子轨迹形态各异，呈非完全周期性也非完全随机响应，表明各故障信号具有显著混沌特性。但因早期振动信号微弱，且受噪声干扰，无法通过吸引子轨迹判断故障类型。

为对原信号降噪提取特征信息，采用OLVMD方法分解原信号，基于Lyapunov指数选取最佳分解参数[c,K]，其中K∈[3,10],c∈[1 000,2 000]，各信号Lyapunov指数如图11所示。由图11可知，以最大Lyapunov指数为寻优条件得到的最佳分解参数组合[c,K]如表6所示。

表6 最佳分解参数

分解各故障信号，各IMF分量时域图见图12。计算各IMF分量的Lyapunov指数，结果见图13。

图13 Lyapunov指数

由图13可知，4种故障的Lyapunov指数各不相同，但均为正数，表明各故障信号具有不同程度的

混沌特性，以Lyapunov指数λ为特征值，λ越大，故障特征越丰富，取λ最大的前3个分量重构故障信号。

为定量表示OLVMD方法的降噪效果和特征提取能力，求取重组信号吸引子轨迹并采用流行学习方法进行降维，可视化故障信息，结果见图14和图15。

(d) 外圈磨损

由图14可知，相较原信号混沌图，4种故障吸引子轨迹收缩，混沌特性增强，信号特征更显著。计算各故障Lyapunov指数，可知由于保持架磨损对轴承运转情况影响最小，其Lyapunov指数最小；多故障导致的混合损伤的Lyapunov指数最大。

图15 降维可视化图

由图15可知，各故障类型间距增大，不存在模态混叠现象，4种故障可明显区分，表明所提方法能很好地提取特征信息。

3.3 关联维数

关联维数描述了振动信号中各点之间的关联性，对信号波动极为敏感，可量化分析故障的状态特征，具有简洁、高效、快速的特点，各故障关联维数与嵌入维数的关系如图16所示。

观察图16可得各故障最佳嵌入维数与关联维数。为准确进行故障识别，需建立故障关联维数区间，通过计算未知故障振动信号关联维数，判断其所处区间，完成故障诊断。

考虑到振动信号本身波动性影响，会出现偏差较大的关联维数点，需剔除这些点以提高故障识别准确率。计算10次关联维数并求其平均值以减小误差，结果如表7所示。

表7 轴承各故障关联维数

由表7可知，各故障关联维数D的区间为：

为更直观地表示各故障关联维数区间，绘制状态区间图，如图17所示。由图17可知，各故障关联维数对应的状态区间存在显著差别，互不相交,可准确完成故障识别。

图17 轴承各故障关联维数区间

综上所述，对轴承进行故障识别的具体步骤如图18所示,其中Dmin和Dmax分别为关联维数的最小值和最大值。

图18 故障识别流程图

3.4 有效性验证

为验证所提方法的有效性，从4种故障中随机选取数据进行验证，为避免出现偏差，每种故障各选取3组数据，计算关联维数，结果如表8所示。

表8 未知故障信号的关联维数

对比表8中故障类型及其关联维数区间，4种故障关联维数平均值均在对应的关联维数区间内，可准确完成故障识别。

3.5 对比验证

为验证所提方法的优越性，基于EEMD和VMD方法分解原信号并重构，采用相同方法求取各故障信号的关联维数区间，结果如图19所示。

图19 不同方法各故障关联维数区间

由图19可知，采用EEMD与VMD方法所得关联维数区间均出现混叠现象，无法准确判断轴承故障类型，具有一定的局限性。

3.6 齿轮箱验证

为验证所采用的方法在实际工程机械轴承中的故障诊断效果，采用美国国家可再生能源实验室(NREL)5 MW近海风力机齿轮箱数据进行验证[15]，对齿轮断齿面积为100%～25%(严重性分为A、B、C和D)4类不同故障情况进行分析，采用OLVMD方法进行降噪并选取敏感分量重构信号，计算各故障齿轮箱振动信号的关联维数，取10次数据的平均值，结果如表9所示。

表9 齿轮箱各故障的关联维数

由表9可知，各断齿故障关联维数区间为：100%断齿D∈[5.200 2,5.318 4]；50%断齿D∈[3.889 6，3.959 8]；25%断齿D∈[4.375 7，4.459 3]；12.5%断齿D∈[2.995 4，3.187 0]。为更直观地表示各故障关联维数区间，绘制状态区间图，如图20所示。