基于互信息的主成分分析结合支持向量回归的滚动轴承剩余寿命预测研究

欧白羽 杨 勐 韩 旭

(1.陆军航空兵学院 陆军航空兵研究所，北京 101121；2.陆军航空兵学院 航空机械工程系，北京 101123；3.93146部队，北京 100000)

引 言

随着现代工业的快速发展，旋转机械设备越来越大型化、复杂化和智能化[1]，如何保证旋转设备的高可靠性安全运行，提高其工作效率，已成为当前亟需解决的关键问题之一。滚动轴承作为旋转机械设备的重要部件，在复杂的工作条件下易发生疲劳损伤[2-3]，进而影响设备的运行效率与安全。因此，对滚动轴承进行有效、可靠、准确的剩余使用寿命预测具有重要的实际意义。基于预测的剩余寿命开展及时的维护工作，可提高设备的整体运转效率，避免了在传统的维护工作中存在的设备得不到及时维护、设备发生故障的风险增加以及健康状态下维护次数过多而造成的资源浪费等问题[4]。

鉴于轴承的振动信号中包含大量的运行状态信息[5]，当前多以设备的振动监测数据为基础进行剩余寿命的预测研究[6]。通常，采集的振动信号中包含数据异常点和噪声，而这些噪声会给预测结果带来误差，因此需要对原始振动信号进行去噪处理。Beale等[7]提出了自适应小波包阈值降噪算法用于消除风力涡轮机叶片的振动噪声。许志华等[8]研究了基于变分模态分解的风机齿轮箱振动信号降噪方法。王大兴[9]针对轴承振动信号提出了经验模态分解与小波包多阈值相结合的降噪方法，有效剔除了信号中的异常数据点。

在实现振动信号的降噪处理后，就需要提取表征轴承退化状态的特征，目前主要使用信号处理的方式提取常见的时域特征、频域特征和时频域特征等[10]。在提取的众多特征中，存在着部分冗余特征和对轴承退化状态不敏感的特征，因此需要对其进行约简，以获得有效的敏感特征作为预测模型的输入，从而提高模型的预测精度。常用的特征约简方法包括相关系数法、互信息法(mutual information, MI)、流形学习[11]和主成分分析(principal component analysis, PCA)[12]等。然而，相关系数法和互信息法只是对提取的多维特征进行筛选[13]，并没有考虑特征之间的冗余性，若特征之间存在信息重叠，则容易产生预测模型的“过拟合”现象，进而影响模型的预测精度[14]，而PCA可以消除特征之间的冗余信息[12]。近年来，智能算法的持续改进与发展为剩余寿命的预测研究提供了新思路，如双向长短时记忆神经网络[15]、改进隐马尔可夫模型[16]、基于Bayesian optimization and hyperband(BOHB)的深度置信网络[17]以及深度卷积神经网络[18]等都取得了较好的预测效果。但这些算法的结构较为复杂且超参数众多，而支持向量回归(support vector regression, SVR)算法的泛化能力和非线性映射能力强且结构简单[19]，能有效解决小样本问题[20]。

基于以上分析，本文提出了一种基于MI-PCA结合SVR的滚动轴承剩余寿命预测方法。通过小波包降噪算法实现原始振动信号的去噪处理，并提取时域、频域和时频域特征；然后计算特征与剩余寿命之间的互信息值进行特征筛选，进而结合PCA算法获取敏感特征用于SVR的输入；最后建立SVR预测模型，实现轴承剩余寿命的精准预测。滚动轴承的全寿命试验数据验证表明，与基于MI和基于PCA的SVR模型相比，所提模型具有更高的预测精度，有利于开展及时的设备维护工作，提高设备的运转效率。

1 剩余寿命预测建模

1.1 建模流程

为了实现滚动轴承剩余使用寿命的精准预测，本文以轴承的振动信号数据为基础，提出了基于MI-PCA结合SVR的剩余寿命预测建模方法。流程包括3个部分：(1)利用小波包降噪算法实现原始振动信号的预处理，并提取表征轴承退化状态的时域、频域和时频域特征；(2)基于特征与剩余寿命的互信息值进行特征筛选，再通过PCA算法获取敏感特征用于SVR模型的输入；(3)建立并训练SVR回归预测模型，实现轴承的剩余寿命预测。具体的建模流程如图1所示。

图1 剩余寿命预测建模流程

1.2 数据降噪与特征提取

由于旋转机械设备结构的复杂性以及工作条件的多样性，导致采集的滚动轴承振动信号中包含大量的噪声[21]，染噪后的信号数据会影响轴承剩余寿命预测的精度。因此，本文首先采用小波包降噪算法对原始振动信号进行去噪处理。小波包降噪的关键在于信号的分解与重构，相比于传统的小波变换，小波包分析不仅对每层中的低频部分进行分解，而且对高频部分也进行分解[22]，能够有效提高信号高频段的分辨率，从而使得信号的表达更为精确。为了进一步挖掘振动信号中的隐含状态信息，使用信号处理的方法从不同角度提取表征轴承退化状态的特征。一般可以从时域、频域和时频域中提取相关特征，以充分利用振动信号数据实现剩余寿命的精准预测。

在提取频域特征时，首先利用快速傅里叶变换将去噪后的时域振动信号x(t)(t=1,2,…,n)转换为频域振动信号y(k)(k=1,2,…,m)，然后提取相关的频域特征。

在提取时频域特征时，首先通过经验模态分解将时域振动信号从高频至低频分解为若干个本征模态函数(intrinsic mode function, IMF)，即IMF分量c(t)(t=1,2,…,n)。为了保证分量信号的有效性，选取前3阶IMF分量的有效值作为时频域特征。本文提取的时域、频域和时频域特征如表1所示。

表1 提取的时域、频域和时频域特征

1.3 敏感特征提取

1.3.1互信息计算

互信息是两个随机变量之间相互依赖的度量,用来衡量两个变量的相互关联程度，表示两个变量之间含有相同信息的部分[23]。假设两个随机变量X和Y，其各自的边缘概率密度函数和联合概率密度函数分别为p(x)，p(y)和p(x,y)，则它们的互信息值I如式(1)所示。

(1)

若X和Y的互信息值为0，表示两个变量为无关变量或相互独立变量，即二者没有相互包含的信息；反之，互信息值越大，表示两个变量之间相互关联的程度越高，包含的相同信息量也就越多。本文通过计算特征与剩余寿命之间的互信息值进行特征筛选，保留与剩余寿命关联度较高的特征。

1.3.2主成分分析

主成分分析法是一种常用的数据降维方法，其基本原理是将具有一定相关性的原始变量经线性变化转换为一组互不相关的变量，且第一变量集中了原始变量所包含的主要信息[24]。基于PCA降维的转换公式如式(2)所示。

(2)

式中，X和Z分别为原始变量和转换后的变量，L为荷载矩阵，即转换矩阵，s为原始变量维数，m为转换后的变量维数，且m≤s。本文首先基于互信息法筛选主要特征，然后利用PCA减少主要特征之间的冗余和重叠信息，从提取的多维特征中获取关键敏感特征，作为预测模型的输入。

1.4 支持向量回归预测模型

支持向量机(support vector machine, SVM)是一种用于解决小样本、非线性及高维模式识别下统计估计和预测问题的机器学习方法[25]。支持向量回归是SVM函数在回归拟合领域的应用。SVR的基本思想是在高维特征空间上寻找一个最优超平面，能够使训练集样本距离该超平面最近，并将该超平面作为回归模型对未知数据进行预测。假设训练样本为(x1,y1),(x2,y2),…,(xn,yn)，xi∈Rn为模型输入，yi∈R为模型期望输出，n为训练样本个数，则回归函数可表示为

(3)

(4)

式中，σ为高斯核函数的指数，影响着SVR模型的预测精度。

1.5 模型评估指标

本文采用平均相对泛化误差EARG(average relative generalization error, ARGE)、均方根误差ERMS(root mean square error, RMSE)和相对误差ER(relative error, RE)来评估模型的回归预测能力，计算公式分别如下。

(5)

(6)

(7)

式中，Pout,k为模型的预测输出，Eout,k为模型的期望输出，n为模型输入样本个数，k表示第k个样本。

2 试验验证

2.1 试验数据

以IEEE PHM12轴承加速性能退化试验的全寿命数据[26]为例验证所提方法的有效性。试验数据来自如图2所示的PRONOSTIA试验台。通过加速度传感器采集轴承在水平和垂直两个方向上的振动信号，每间隔10 s采集一次数据，设置采样频率为25.6 kHz，每次采样0.1 s，因此每次可采集2 560个振动加速度。该数据集共包含3种工况下17个滚动轴承的全寿命周期试验数据，工况划分如表2所示。由于轴承水平方向振动信号包含更多的状态信息[26]，因此本文以工况1下轴承1水平方向的振动信号为研究对象，构建基于MI-PCA结合SVR的滚动轴承剩余寿命预测模型。

表2 工况划分设置

图2 PRONOSTIA试验台

2.2 数据降噪与敏感特征提取

首先，采用小波包降噪算法处理轴承1的原始振动信号，其原始振动信号及降噪后的振动信号如图3所示。从图中可以看出降噪后的数据更加平滑，剔除了异常数据点，并保留了原始振动信号中的有效信息。

图3 原始信号和小波包降噪信号的对比

为了更好地突出显示小波包降噪的效果，选取某一次采集的振动信号进行分析，同时将时域振动信号转换为频域信号，以便更直观地分析降噪效果。图4为小波包降噪前后的频域信号对比，可以看出小波包降噪方法可以在一定程度上滤除噪声，而且能够保留原始信号的有效信息，即降噪后的频域信号含有峰值，可表征轴承的磨损状况[27]，从而可提取出对轴承退化敏感的特征。

图4 原始信号和小波包降噪信号的频域信号对比

图5为提取的18个时域、频域和时频域特征。从图5可以看出，部分特征的退化趋势不明显，无法表征轴承的退化状态，因此首先采用互信息法筛选对轴承退化状态较为敏感的特征。表3为提取的特征与真实剩余寿命的互信息计算结果，将互信息的平均值设置为阈值，保留互信息值大于阈值的特征。

图5 降噪信号的时域、频域和时频域特征

表3 特征与真实剩余寿命的互信息值

本文案例计算的阈值为7.28，剔除7个对轴承退化状态不敏感的特征，然后基于保留的11个特征进行PCA降维，选取累积贡献率大于90%的第一主元变量作为敏感特征[28]，结果如图6所示。从图6可以看出，敏感特征可以很好地表征轴承的退化状态，即正常状态、轻微退化、严重退化和失效状态。

图6 敏感特征

2.3 预测结果分析

将提取的敏感特征作为SVR模型的输入，将轴承的真实剩余寿命作为SVR模型的期望输出，对模型进行训练和测试。选取4种退化状态中的各50个样本进行测试，为了验证基于MI-PCA的SVR模型预测结果的准确性，分别建立基于MI和基于PCA的SVR回归预测模型对相同的样本进行测试。

表4为3种预测模型的评估指标结果。从表4可以看出，3种模型对前3种退化状态剩余寿命的预测误差较低，对失效状态的预测误差较大，然而MI-PCA-SVR模型的预测误差要远小于另外两种模型。总体来看，MI-SVR模型、PCA-SVR模型以及MI-PCA-SVR模型的平均预测准确率分别为90%、93%和97%，MI-PCA-SVR模型的预测准确率最高，预测误差最低，说明采用MI-PCA方法对输入样本进行降维处理可提高模型的预测精度。

表4 3种模型评估结果

图7为3种模型的相对误差曲线。从图7可以看出，MI-PCA-SVR模型在严重退化和失效状态下的预测误差均小于MI-SVR模型和PCA-SVR模型。此外，MI-SVR模型的相对误差曲线波动较大，其相对误差最高可达50%，而MI-PCA-SVR模型的相对误差曲线较平稳，总体误差均在13%以内。结合表4的RMSE计算结果和图7的误差曲线，说明MI-PCA-SVR模型的稳定性要优于另外两种模型。基于以上分析，采用基于MI-PCA的SVR模型能够实现滚动轴承剩余寿命的精准预测。

图7 3种模型的相对误差曲线

2.4 IMS轴承数据案例验证

为进一步验证所提方法的有效性，采用美国辛辛那提大学智能维护系统中心(Intelligent Maintenance System, IMS)的滚动轴承全寿命周期试验数据进行案例分析与验证。试验数据来自图8所示的试验台，通过加速度传感器采集轴承的振动信号，每间隔10 min采集一次数据，设置采样频率为20 kHz，采样点数为204 800。试验结束后共取得3个数据集。本文以数据集2中轴承1运行的全寿命周期加速疲劳振动数据作为研究对象，构建基于MI-PCA结合SVR的滚动轴承剩余寿命预测模型。

图8 IMS试验台

为体现本文所提方法的全局优势，分别使用MI-SVR模型、PCA-SVR模型和MI-PCA-SVR模型预测轴承1的全过程剩余使用寿命。3种方法的预测结果如图9所示。从图9可以看出，MI-SVR模型和PCA-SVR模型在轴承的正常运行及初始退化阶段预测值和实际值之间的误差较大，随着时间的推移，在轴承的退化后期，误差逐渐缩小，而MI-PCA-SVR模型的整体预测误差均较小，预测值在实际值上下波动。此外，为了衡量模型预测结果的稳定性，计算得到MI-SVR模型、PCA-SVR模型和MI-PCA-SVR模型的RMSE值分别为14.922 7、20.065 2和2.609 0，表明本文所提方法的剩余寿命预测结果更加稳定。通过分析3种模型的预测结果曲线以及RMSE评估指标值可以得出，本文提出的MI-PCA-SVR模型的预测精度更高，而且预测结果的稳定性最优。

图9 3种模型的剩余寿命预测结果

3 结论

为了提高旋转机械设备的运转效率，实现设备的智能维护，本文提出了一种基于MI-PCA结合SVR模型的滚动轴承剩余寿命预测方法。首先利用小波包降噪算法消除原始振动信号中的噪声点对预测精度的影响，并基于降噪数据提取时域、频域和时频域特征。然后计算提取的特征与轴承真实剩余寿命的互信息值进行特征筛选，进而通过PCA降维算法剔除筛选特征之间的冗余信息，获取表征轴承退化状态的敏感特征。最后以敏感特征作为模型的输入数据实现SVR模型的训练和测试。试验结果表明，基于MI-PCA的敏感特征提取方法可以充分挖掘振动信号中的状态信息，有效区分轴承退化过程中的不同状态；且基于MI-PCA构建的SVR模型可实现正常状态、轻微退化、严重退化和失效状态这4种退化状态下剩余寿命的精准预测，为开展及时有效的设备维护工作提供了有力依据。

符号说明

x(t)—时域振动信号，mm/s2

y(t)—频域振动信号，mm/s2

c(t)—IMF分量，mm/s2

Xss—标准差，mm/s2

Xpp—峰峰值，mm/s2

Xrms—均方根，mm/s2

Xs—方差，mm2/s4

CL—偏度因子

Kr—峭度因子

C—峰值因子

S—波形因子

I—脉冲因子

M—边缘因子

f1—平均频率，mm/s2

f2—频率均方根，mm/s2

f3—频率峰值，mm/s2

I1—IMF1有效值，mm/s2

I2—IMF2有效值，mm/s2

I3—IMF3有效值，mm/s2