基于PCA多模型融合的滚动轴承性能退化指标构建

Construction of performance degradation indicators for rolling bearings based on PCA multi-model fusion

摘要

单模型构建的滚动轴承性能健康指标仅能从本身的“单角度”来描述滚动轴承的性能退化状态，具有一定的局限性。为解决这个问题，提出一种基于主成分分析（principal component analysis，PCA）多模型融合的滚动轴承健康指标构建方法。该方法分别采用支持向量数据描述（support vector data description，SVDD）模型、自联想核回归（auto-associative kernel regression，AAKR）模型和高斯混合模型（Gaussian mixture module，GMM）构建相应单模型的健康指标，再将3个单模型的健康指标经主成分分析（PCA）融合，并选取第一主成分作为能够包含“多角度”性能退化信息的健康指标（SAG-HI）。试验结果表明，相比于各单模型的健康指标，SAG-HI与滚动轴承保持可靠度的灰置信水平达到98.38%，其相关性、单调性和鲁棒性也均表现为最优，且通过包络谱分析验证了其能够准确且及时监测到早期故障发生时刻。

关键词

滚动轴承 / 支持向量数据 / 自联想核回归 / 高斯混合模型 / 主成分分析 / 性能退化指标 / 多模型融合

Key words

rolling bearings / support vector data / auto-associative kernel regression / Gaussian mixture model / PCA / performance degradation indicators / multi-model fusion

准确表征轴承运行中的性能退化程度，对保障系统的稳定运行至关重要。基于单模型的指标构建方法通过融合多种特征来共同描述轴承的性能退化状态。周建民等^［1］用无故障轴承振动信号的小波包奇异谱熵训练SVDD模型，并以此构建基于距离型的轴承性能退化指标。此指标在退化阶段描述退化趋势的能力较强，但SVDD模型对超球体外的特征向量表现过于敏感，易造成早期故障虚警的发生。Rai等^［2］通过集成经验模态分解（ensemble empirical mode decomposition，EEMD）提取了健康条件下轴承振动信号的固有模态分量的奇异值特征和能矩熵特征，以构建基准GMM，并输入全寿命样本的特征向量构成模型，将输出的后验概率的杰森-瑞丽散度（Jensen-Rényi divergence，JRD）作为基于概率的性能退化指标，此退化指标在健康阶段稳定性较好，但在退化阶段描述退化趋势的能力较弱。刘路路^［3］用轴承健康样本的特征向量组成AAKR模型的健康记忆矩阵，并将全寿命样本的特征向量作为AAKR模型中的观察向量，以构建基于空间映射的性能退化指标。此指标虽基本能刻画出性能退化的趋势，但总体的鲁棒性较差。另外还可以用基于数据降维的方法来构建多特征融合指标。Widodo等^［4］将滚动轴承的峰值、峭度和熵估计值，经PCA方法进行降维，并选取第一主成分作为轴承性能退化指标。

上述方法在构建健康指标前，没有考虑到特征之间会存在冗余性且未对各特征的重要程度作出评价。故本文利用类softmax函数^［5］对经过多尺度排列熵的均值（mean of multiscale permutation entropy，MMPE）阈值^［6］筛选后的特征赋予相应的权值。针对上述单模型健康指标在表征轴承全寿命区间内的性能退化状态时，均具有各自的优缺点。故基于优势互补的思想，本文提出了一种将基于距离型的单模型健康指标、基于空间映射的单模型健康指标和基于概率的单模型健康指标经PCA方法融合，并以分散程度最大的第一主成分分量构建多模型融合的健康指标。依据滚动轴承加速寿命实验数据集的实验结果，验证了本文所提方法的有效性和优越性。

1 理论基础

1.1 SVDD模型构建健康指标

SVDD基本思想是寻求一个最优超球体，使目标样本尽可能包含在球内，非目标样本尽可能排除在球外^［7］，该球体的半径和球心分别用 R和 O表示。用健康样本的特征向量 \left\{{\mathit{x}}_h,h=\mathrm{1},\mathrm{2},\cdots ,H\right\}构建模型训练集。依据支持向量机理论，最优超球体的优化目标和需满足的约束条件为

\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}(R^{\mathrm{2}}+C{\displaystyle \sum _{h=\mathrm{1}}^H}{\xi }_h)\\ \mathrm{s}.\mathrm{t}.\Vert {\mathit{x}}_h-O{\Vert }^{\mathrm{2}}\le R^{\mathrm{2}}+{\xi }_h({\xi }_h\ge \mathrm{0},h=\mathrm{1},\mathrm{2},\cdots ,H)\end{array} （1）

式中： C为惩罚参数； {\xi }_h为松弛变量。

引入核函数 K ( {\mathit{x}}_j , {\mathit{x}}_k )（本文使用高斯核函数）可以解决在低维空间线性不可分的问题。超球体半径 R为

\begin{array}{l}{R}^{\mathrm{2}}=K({\mathit{x}}_s,{\mathit{x}}_s)-\mathrm{2}{\displaystyle \sum _{h=\mathrm{1}}^H}{\alpha }_{h}K({\mathit{x}}_h,{\mathit{x}}_s)+\\ {\displaystyle \sum _{j=\mathrm{1}}^H}{\displaystyle \sum _{k=\mathrm{1}}^H}{\alpha }_j{\alpha }_{k}K({\mathit{x}}_j,{\mathit{x}}_k)\end{array} （2）

式中： {\alpha }_{}为拉格朗日乘子； {\mathit{x}}_s为支持向量。任意一个全寿命样本中第 i个样本的特征向量 {\mathit{x}}_i到超球体中心 O的距离 D_i^{}为

\begin{array}{l}{D}_i^{\mathrm{2}}=K({\mathit{x}}_i,{\mathit{x}}_i)+{\displaystyle \sum _{j=\mathrm{1}}^H}{\displaystyle \sum _{k=\mathrm{1}}^H}{\alpha }_j{\alpha }_{k}K({\mathit{x}}_j,{\mathit{x}}_k)-\\ \mathrm{2}{\displaystyle \sum _{h=\mathrm{1}}^H}{\alpha }_{h}K({\mathit{x}}_i,{\mathit{x}}_h)\end{array} （3）

通过比较轴承的全寿命样本的特征向量 {\mathit{x}}_i到球心的距离 D_i与超球体的半径 R之间的差值来描述轴承性能退化的状况。即由SVDD模型构建的健康指标为

H I_{\mathrm{S}\mathrm{V}\mathrm{D}\mathrm{D}} ( i ) = \left\{\begin{array}{ll}\mathrm{0}& ,D_i-R\le \mathrm{0}\\ D_i-R& ,D_i-R>\mathrm{0}\end{array}\right. （4）

1.2 AAKR模型构建HI

AAKR模型构建 \mathrm{H} \mathrm{I} ^［8］的基本思想是将轴承健康状态下样本的特征向量构成健康记忆矩阵，将轴承全寿命样本的特征向量（观测向量）用健康记忆矩阵重构，可得到对应的重构估计向量，然后量化观测向量和其对应的重构估计向量之间的差异性，以构建表征轴承性能退化程度的指标。

对于含有 G个特征参数的健康样本的特征向量所构建的健康记忆矩阵 \mathit{E}为

\mathit{E} = \left[\begin{array}{cccc}{x}_{\mathrm{1,1}}& x_{\mathrm{1,2}}& \cdots & x_{\mathrm{1},G}\\ x_{\mathrm{2,1}}& x_{\mathrm{2,2}}& \cdots & x_{\mathrm{2},G}\\ ⋮& ⋮& \ddots & ⋮\\ x_{H,\mathrm{1}}& x_{H,\mathrm{2}}& \cdots & x_{H,G}\end{array}\right] （5）

式中： x_{i,j}是第 i个健康样本的特征向量中的第 j个特征参数值。

对于任意观测向量 {\mathit{x}}_i = {\left[\begin{array}{cccc}{x}_{i,\mathrm{1}}& x_{i,\mathrm{2}}& \cdots & x_{i,G}\end{array}\right]}^{\mathrm{T}}与记忆矩阵 \mathit{E}中每个健康样本的特征向量间的欧式距离为

d_z ( {\mathit{E}}_z , {\mathit{x}}_i ) = \sqrt[]{{\displaystyle \sum _{j=\mathrm{1}}^G}{\left({\mathit{E}}_{z,j}-x_{i,j}\right)}^{\mathrm{2}}} ( z = \mathrm{1,2} , \cdots , H ) （6）

其在健康空间 \mathit{E}的重构向量 {\tilde{\mathit{x}}}_i为

{\tilde{\mathit{x}}}_i = \frac{{\displaystyle \sum _{z=\mathrm{1}}^H}{w}_z{\mathit{E}}_z}{{\displaystyle \sum _{z=\mathrm{1}}^H}{w}_z} （7）

式中： w_z为权重，计算公式为

w_z = \frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{q}^{\mathrm{2}}}} e^{\frac{-{d_z}^{\mathrm{2}}}{\mathrm{2}{q}^{\mathrm{2}}}} （8）

式中： q是核函数的宽度系数。因此基于全寿命样本中的第 i个观测向量x_i 与其对应的重构估计向量 {\tilde{\mathit{x}}}_i之间的欧氏距离所构建的健康指标为

H I_{\mathrm{A}\mathrm{A}\mathrm{K}\mathrm{R}} ( i ) = {‖{\mathit{x}}_i-{\tilde{\mathit{x}}}_i‖}^{\mathrm{2}} （9）

式中： {‖•‖}^{\mathrm{2}}表示两向量间欧氏距离运算。

1.3 GMM构建HI

GMM是一种将多个单高斯概率密度函数加权组合的参数模型，可描述数据样本在高维空间中的非椭球体分布。对于健康状态下样本的特征向量 \left\{{\mathit{x}}_h,h=\mathrm{1},\mathrm{2},\cdots ,H\right\}，其GMM概率密度 p ( \mathit{x} )定义为

p ( \mathit{x} ) = {\displaystyle \sum _{j=\mathrm{1}}^M} b_j l ( \mathit{x} | {\mathit{\mu }}_j , {\mathit{\sigma }}_j ) （10）

式中： M为高斯模型的数量； b_j为各模型的混合权值，且 {\displaystyle \sum _{j=\mathrm{1}}^M} b_j = \mathrm{1}； {\mathit{\mu }}_j为第 j个高斯分布函数的均值向量； {\mathit{\sigma }}_j为协方差矩阵； l ( \mathit{x} | {\mathit{\mu }}_j , {\mathit{\sigma }}_j )为单个高斯概率密度函数。

以最大化GMM的似然函数为目标，利用EM算法^［9］实现对GMM模型中的均值向量、混合权重和协方差矩阵的参数估计如式（11）所示。

J \left(b_j,{\mathit{\mu }}_j,{\mathit{\sigma }}_j\right) = {\displaystyle \sum _{h=\mathrm{1}}^H} ln ( {\displaystyle \sum _{j=\mathrm{1}}^M} b_j(l ( {\mathit{x}}_j | {\mathit{\mu }}_j , {\mathit{\sigma }}_j ) ) （11）

对于第 i个全寿命样本的特征矢量 {\mathit{x}}_i，计算其GMM概率密度 P ( {\mathit{x}}_i )，以表示 {\mathit{x}}_i属于由健康信号构成的GMM的概率大小，引入多域对数概率^［10］作为健康指标，其计算公式为

H I_{\mathrm{G}\mathrm{M}\mathrm{M}} ( i ) = \frac{\mathrm{l}\mathrm{n}(P({\mathit{x}}_i)+{\mathrm{10}}^{-\mathrm{22}}}{\mathrm{l}\mathrm{n}({\mathrm{10}}^{-\mathrm{22}})} （12）

2 基于PCA多模型融合的性能退化指标构建

2.1 多域特征提取

本文分别从时域、频域、时频域和熵的角度对轴承的振动加速度信号共提取27个原始特征来构建候选特征集，如表1所示。其中，时域特征^［11］可分为有量纲和无量纲特征共12个，频域特征共4个。由于相比于小波包分解系数，小波包节点能量比抗干扰的能力更强，故本文选用db10正交小波基，对滚动轴承振动信号进行三层小波包分解，将得到的8个子频带的能量比作为时频域特征^［12］，熵特征^［13］共3个。

2.2 特征筛选和赋权

本文选用MMPE作为特征筛选的标准，MMPE即计算不同尺度下时间序列的排列熵的均值，其能够准确刻画出特征参数在多尺度上的复杂性与随机性，特别是对突变信息敏感性较强。同时为正确评价重要的特征成分，本文采用类softmax函数进行特征赋权，权值计算方法如式（13）所示。

{\mathit{W}}_i = \frac{e^{\frac{\mathrm{1}}{\mathrm{M}\mathrm{M}\mathrm{P}{\mathrm{E}}_i}}}{{\displaystyle \sum _{j=\mathrm{1}}^n}{e}^{\frac{\mathrm{1}}{\mathrm{M}\mathrm{M}\mathrm{P}{\mathrm{E}}_j}}} （13）

式中： {\mathit{W}}_i为第 i个特征的权重； \mathrm{M} \mathrm{M} \mathrm{P} {\mathrm{E}}_i为第 i个特征的MMPE值。

2.3 多模型融合的性能退化指标构建方法与步骤

PCA以方差来衡量信息量的多少，在降低数据维度的同时最大程度保留了3种单模型健康指标中的有效信息，去除了3种单模型健康指标中具有较强相关性的冗余信息。PCA具有强化关键信息、计算开销较小和稳健性较强等优点，更加符合实际工程应用的需求。鉴于3种单模型的健康指标在全寿命区间内表征轴承性能退化状态时各自表现出的优势与不足之处，本文基于优势互补的思想，以3种单模型的健康指标为基础，采用PCA数据融合方法构建多模型融合的性能退化指标，该方法的主要步骤如下：

步骤1：收集轴承全寿命振动数据，确定样本长度 n和样本数量 m，原始数据样本记为 {\mathit{s}}_i \in {\mathbf{R}}^n， i = \mathrm{1} , \mathrm{2} , \cdots , m。

步骤2：多域特征筛选与赋权。假设从表1中共筛选出 g种特征，并根据式（13）对这 g种特征进行赋权，对原始数据样本 {\mathit{s}}_i提取这g种特征，组成的特征向量记为 {\mathit{x}}_i \in {\mathbf{R}}^g。

步骤3：构建单模型健康指标向量为

{\mathit{h}}_i = {\left[\begin{array}{ccc}H{I}_{\mathrm{S}\mathrm{V}\mathrm{D}\mathrm{D}}(i)& H{I}_{\mathrm{A}\mathrm{A}\mathrm{K}\mathrm{R}}(i)& H{I}_{\mathrm{G}\mathrm{M}\mathrm{M}}(i)\end{array}\right]}^{\mathrm{T}} （14）

步骤4：由m组样本构建矩阵 H

\mathit{H} = {\left[\begin{array}{ccc}{\mathit{h}}_{\mathrm{1}}& {\mathit{h}}_{\mathrm{2}}& \begin{array}{cc}\cdots & {\mathit{h}}_m\end{array}\end{array}\right]}^{\mathrm{T}} （15）

步骤5：构建PCA模型。首先对 H 进行标准化处理，即计算样本均值向量 {\overline{\mathit{h}}}_{}和标准差向量 \mathit{\sigma }，每个样本 {\mathit{h}}_i标准化处理后组成的矩阵记为 \tilde{\mathit{H}}，根据 \tilde{\mathit{H}}求出协方差矩阵 A，求出 A 特征值的最大值为 {\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}和其对应的特征向量 {\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}，则变换矩阵 \mathit{B} = {\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}。计算所有样本的主元得分向量 \mathit{t}，如式（16）所示，找出 \mathit{t}中的最大值 t_{\mathrm{m}\mathrm{a}\mathrm{x}}和最小值 t_{\mathrm{m}\mathrm{i}\mathrm{n}}。

\mathit{t} = \tilde{\mathit{H}} \mathit{B} （16）

步骤6：对任意一个原始数据样本的特征向量 {\mathit{x}}_i \in {\mathbf{R}}^g，根据步骤3构建单模型健康指标向量 {\mathit{h}}_i，对其进行标准化处理，处理后的健康指标向量记为 {\tilde{\mathit{h}}}_i，则多模型融合的健康指标定义为

H I_{\mathrm{S}\mathrm{A}\mathrm{G}} ( i ) = \frac{{{\tilde{\mathit{h}}}_i}^{\mathrm{T}}\mathit{B}-t_{\mathrm{m}\mathrm{i}\mathrm{n}}}{t_{\mathrm{m}\mathrm{a}\mathrm{x}}-t_{\mathrm{m}\mathrm{i}\mathrm{n}}} （17）

3 实验验证

3.1 实验数据描述

采用西安交通大学机械装备健康监测联合实验室的滚动轴承加速寿命测试平台试验数据对本文提出的方法进行验证，轴承加速寿命实验台如图1所示。实验轴承型号为LDKUER204。实验中设置采样频率为25.6 kHz，采样间隔为1 min，每次采样时长为1.28 s。实验共设计了3类工况，每类工况下有5个轴承。本案例采用工况1的Bearing1_1数据，实际寿命123 min，最终失效位置为轴承外圈，共采集123组样本。本文每组样本选取2 000个代表数据。

3.2 特征预处理和SAG-HI构建

计算得到27个特征的MMPE值如图2所示。选取筛选阈值MMPE=2.6，共筛选出6个优选特征，它们分别是方根幅值F2、标准差F3、均方根值F6、平均频率F13、第4个子频带能量比F20和模糊熵F25。并依据式（13）计算得到各优选特征的权值如表2所示。

选取轴承全寿命样本中的前30组样本作为各模型健康训练样本集，设定各模型参数为： q = \mathrm{1}、 C = \mathrm{0.3}和 M = \mathrm{1}，构建 H I_{\mathrm{S}\mathrm{A}\mathrm{G}}，经计算得到第一主成分的贡献率^［14］为95.45%。这说明第一主成分已几乎全部包含了3个单模型HI中全部的轴承性能退化信息。同时本文采用 \mathrm{3} \sigma准则确定故障阈值，SAG-HI及其故障阈值如图3所示。

从图3可知，SAG-HI判定的初始退化时刻对应第65个样本。为验证其对初始退化点监测的准确性，分别对第64和第65个样本进行包络谱分析，结果分别如图4和图5所示。在图4中未发现与外圈故障特征频率 f = \mathrm{107.9} \mathrm{H} \mathrm{z}接近的频率成分。但在图5中可找到与 f接近的频率成分 \mathrm{108.8} \mathrm{H} \mathrm{z}，其幅值大小不太突出且包络谱中“地毯”噪声较大，可以确定轴承此时发生了早期外圈故障。这充分验证了SAG-HI能够准确地监测到早期故障的发生。

3.3 各模型健康指标对比与分析

为了定量评价各模型的健康指标在全寿命区间表征轴承性能退化状况的能力水平，将各模型的健康指标与轴承保持可靠性做灰关系评估^［15］。本文在置信水平为95%条件下，求得本征序列（健康样本序列）的最大熵估计区间为［-1.255，0.894］。以此可求得基于泊松过程的轴承保持可靠度序列 R， R以及各模型的健康指标经归一化处理结果如图6所示。

在灰关系权值系数 \omega = \mathrm{0.5}条件下^［16］，求得各模型的健康指标与保持可靠度序列 R之间的灰置信水平，如表3所示。

从图6可观察到保持可靠度序列 R与SAG-HI的趋势走向最为贴近，且由表3可知，SAG-HI与保持可靠度序列R的灰置信水平最高。这充分表明了SAG-HI更符合实际的轴承性能退化演变规律。

采用相关性、单调性与鲁棒性指标来评价各模型健康指标在全寿命区间上的综合表现^［17］。计算得到各模型HI的相关性、趋势性与鲁棒性大小如表4所示。从表4可观察到SAG-HI的相关性（0.97）、单调性（0.42）以及鲁棒性（0.86）均表现为最优。证明了本文提出的多模型融合的性能退化指标构建方法的优越性。

4 结论

本文所提出的SAG-HI可从距离、空间映射和概率3个角度共同表征轴承的性能退化状态。从上述实验可知，SAG-HI克服了单模型HI在表征性能退化指标时的本身局限性，更符合实际的轴承性能退化演变规律，同时能够及时准确地监测到早期故障的发生，这在工程实践中具有重要意义。

参考文献

[1] 周建民，徐清瑶，张龙，结合小波包奇异谱熵和SVDD的滚动轴承性能退化评估［J］.机械科学与技术，2016，35（12）：1882-1887.

[2] Rai A， Upadhyay S.An integrated approach to bearing prognostics based on EEMD-multi feature extraction，Gaussian mixture models and Jensen-Renyi divergence［J］. Applied Soft Computing Journal，2018，71（2）：36-50.

[3] 刘路路.基于特征融合和LSTM的滚动轴承剩余寿命预测［D］.成都：电子科技大学， 2021.

[4] Widodo A， Yang B.Application of relevance vector machine and survival probability to machine degradation assessment［J］.Expert Systems with Application，2011，38（3）：2592-2599.

[5] 张康.基于全寿命信息与深度时序特征的滚动轴承健康指标构建方法研究［D］.上海：上海大学，2022.

[6] 董治麟.多尺度排列熵方法及其在滚动轴承智能故障诊断中的应用研究［D］.马鞍山：安徽工业大学，2020.

[7] Alam S， Sonbhadra K S， Agarwal S，et al.Sample reduction using farthest boundary point estimation （FBPE） for support vector data description （SVDD）［J］.Pattern Recognition Letters，2020，131（3）：268-276.

[8] 田雯雯，吕丽霞，冯雪凯，基于改进AAKR的风电机组齿轮箱状态监测［J］.电子测量技术，2022，45（15）：158-165.

[9] 魏居辉，王炯琦，穆京京，基于高斯混合模型的卫星电源系统异常检测方法［J］.空间控制技术与应用，2022，48（4）：104-114.

[10] 张龙，黄文艺，熊国良，基于多域特征与高斯混合模型的滚动轴承性能退化评估［J］.中国机械工程，2014，25（22）：3066-3072.

[11] Stack J R， Habetler T G， Harley R G.Fault-signature modeling and detection of inner-race bearing faults［J］.IEEE Transactions on Industry Applications，2006，42（1）：61-68.

[12] 孙建.滚动轴承振动故障特征提取与寿命预测研究［D］.大连：大连理工大学，2016.

[13] Jin de Z， Tao L， Rui M，et al.Generalized composite multiscale permutation entropy and PCA based fault diagnosis of rolling bearings［J］.Journal of Vibration and Shock，2018，37（20）：61-66.

[14] 耿俊超，王天真，韩金刚，基于多级特征提取的并网逆变器故障诊断策略［J］.电源学报，2022，20（4）：28-36.

[15] 程立，马文锁，夏新涛，基于改进模糊熵和灰关系的滚动轴承性能退化评估［J］.机械传动，2022，46（1）：56-64，86.

[16] 叶亮，夏新涛，常振.滚动轴承振动性能的混沌特性与不确定性之间的灰关系评估［J］.机械传动，2021，45（1）：9-16，103.

[17] 常永祺. 机械旋转部件性能退化跟踪与剩余使用寿命预测方法研究［D］.哈尔滨：哈尔滨工业大学，2020.

基金资助

国家自然科学基金(62003223)