目前，智能预测理论在机械设备的运行状态评估及性能退化趋势预测中取得了大量研究成果，对机械设备早期故障的预防性维修起到了积极作用^[1–2]。传统的智能预测方法以神经网络和支持向量机（support vector machine, SVM）等预测模型最为典型^[3]，然而，神经网络需要大量均匀分布的训练样本方可提高预测精度，支持向量机虽然能够克服神经网络的局限性，但是SVM的核函数必须满足Mercer条件，并且其正则化参数难以确定。

为了克服传统智能预测模型存在的问题，Tipping提出了一种基于贝叶斯理论的概率学习模型—相关向量机（relevance vector machine, RVM）^[4–5]。较之神经网络预测模型，RVM需要的训练样本数量明显减少，预测精度更高；较之SVM预测模型，RVM克服了SVM正则化系数确定困难、核函数受Mercer条件限制等固有局限。在基于RVM的预测过程中，RVM的核函数对应不同的预测函数，进而直接影响RVM的预测精度^[6–7]。传统RVM研究中采用单个核函数，能够有效解决常规时序数据的趋势预测问题，但是对于高维非线性、小子样、强干扰等复杂数据源的机械设备运行状态数据趋势预测存在一定的局限性。

针对以上问题，提出了一种组合核相关向量机（combined kernel relevance vector machine, C-RVM）机械设备运行状态的趋势预测方法。采用归一化EEMD信息熵评估机械设备服役过程中的运行状态，依据机械设备运行状态指标参数的特性，自适应调整RVM的核函数参数，获得最优的组合核相关向量机预测模型，以此实现多特征空间的有效融合，最终实现机械设备运行状态趋势的有效预测。滚动轴承运行状态指标的预测结果表明，相对于单一核函数相关向量机以及其他智能预测模型，组合核RVM提高了滚动轴承运行状态指标的预测精度，可以满足机械设备运行状态的在线评估和性能退化趋势有效预测的实际应用需求；同时，组合核RVM还能够对滚动轴承性能退化指标的误差范围进行表征，具有较强的实用性。

1 组合核相关向量机可靠度预测 1.1 相关向量机

给定样本数据集 $D = \left\{ {\left( {{x_i},{t_i}} \right)|i = 1,2, \cdots ,m} \right\}$ ，其中， ${x_i} \in {R^n}$ 为输入变量， ${t_i} \in R$ 为相应输出， $m$ 为样本数，则相关向量机回归函数为：

$f\left( {{{x}};{{w}}} \right) = \sum\limits_{i = 1}^M {\left[ {{w_i}K\left( {{{x}},{x_i}} \right) + {w_0}} \right]} $

(1)

式中， ${{w}} = {\left[ {{w_0},{w_1}, \cdots ,{w_M}} \right]^{\rm T}}$ 为对应权值， $K\left( {{{x}},{x_i}} \right)$ 为核函数。针对普遍存在的噪声、误差等问题，相关向量机在目标值中加入了噪声影响，输入矢量与目标值之间的关系变为：

${t_i} = f\left( {{{x}};{{w}}} \right) + {\varepsilon _i}$

(2)

式中， ${\varepsilon _i}$ 为随机噪声，且 ${\varepsilon _i}{\text{～}} N\left( {0,{\sigma ^2}} \right)$ 。设 ${t_i} \in R$ 为独立的随机变量，则对应整个训练样本集的似然函数为：

$p\left( {t\left| {w,{\sigma ^2}} \right.} \right) = {\left( {2{\text{π}} {\sigma ^2}} \right)^{ - \frac{N}{2}}}\exp \left\{ { - \frac{1}{{2{\sigma ^2}}}{{\left\| {{ t} - { \varPhi} { w}} \right\|}^2}} \right\}$

(3)

式中， ${ t} \!=\! {\left[ {{t_1},t_2, \cdots ,{t_N}} \right]^{\rm T}}$ ， ${ w} \!=\! {\left[ {{w_1}, w_2,\cdots ,{w_N}} \right]^{\rm T}}$ ， ${ \varPhi} \!=\![\varPhi \left( {{x_1}} \right),$ $\varPhi \left( {{x_2}} \right), \cdots\!,\varPhi \left( {{x_N}} \right)]^{\rm T}$ ， $\varPhi \left( {{x_i}} \right) \!=\! ( 1,k\left( {{x_i},{x_1}} \right),k\left( {{x_i},{x_2}} \right), \cdots\!,k( {x_i},$ $ {x_N} ) )^{\rm T}$ 。根据概率预测公式，所求的目标值 ${t_ * }$ 的条件概率为：

$p\left( {{t_ * }\left|{ t} \right.} \right) = \int {p\left( {{t_ * }\left| {w,{\sigma ^2}} \right.} \right)} p\left( {w,{\sigma ^2}\left| { t} \right.} \right){\rm d}w{\rm d}{\sigma ^2}$

(4)

为了求解 $ w$ 和 ${\sigma ^2}$ ，根据贝叶斯理论，假设 $ w$ 服从标准正态分布，则 $ w$ 的条件概率分布为：

$p\left( {{ w}\left| \alpha \right.} \right) = \prod\limits_{i = 0}^N {N\left( {{w_i}\left| {0,\sigma _i^{ - 1}} \right.} \right)} $

(5)

由此，式（4）可更改为：

$ p\left( {{t_ * }\left| t \right.} \right) = \int {p\left( {{t_ * }\left| {{ w},{\sigma ^2}} \right.} \right)} p\left( {{ w},\alpha ,{\sigma ^2}\left| t \right.} \right){\rm d}w{\rm d}\alpha {\rm d}{\sigma ^2} $

(6)

根据贝叶斯公式，对 $p\left( {{ w},\alpha ,{\sigma ^2}\left| t \right.} \right)$ 进行整理可得：

$p\left( {{ w},\alpha ,{\sigma ^2}\left|{ t} \right.} \right) = p\left( {{ w}\left| {{ t},\alpha ,{\sigma ^2}} \right.} \right)p\left( {\alpha ,{\sigma ^2}\left|{ t} \right.} \right)$

(7)

权值矢量 $ w$ 的后验概率为：

$\begin{aligned}[b] & \quad p\left( {{ w}\left| {t,\alpha ,{\sigma ^2}} \right.} \right) = \frac{{p\left( {{ w},\alpha ,{\sigma ^2}\left|{ t} \right.} \right)}}{{p\left( {\alpha ,{\sigma ^2}\left| { t} \right.} \right)}} = \\ & {\left( {2{\text{π}} } \right)^{ - \frac{{N + 1}}{2}}}{\left| {\varSigma } \right|^{ - \frac{1}{2}}}\exp \left( { - \frac{1}{2}{{\left( {{ w} - \mu } \right)}^{\rm T}}\varSigma^{ - 1} \left( {{ w} - \mu } \right)} \right) \end{aligned} $

(8)

式中， $\varSigma \! = \!{\left( {{\sigma ^{ - 2}}{{ \phi} ^{\rm T}}{ \phi} \!+\! A} \right)^{ - 1}}$ ， $\mu = {\sigma ^{ - 2}}\varSigma {{{ \phi} ^{\rm T}}t} $ ， $A \!=\! {\rm diag}( {\alpha _0},{\alpha _1}, \cdots ,$ ${\alpha _N})$ 。引入delta近似函数 $p\left( {\alpha ,{\sigma ^2}\left| t \right.} \right) \approx \delta \left( {{\alpha _{\rm MP}},\sigma _{\rm{MP}}^2} \right)$ ，通过近似运算来求解后验模式问题式（7）。此时采用最大似然估计，用delta近似函数的峰值来逼近后验概率 $p\left( {\alpha ,{\sigma ^2}\left|{ t }\right.} \right)$ ，通过最大化 $p\left( {{ t}\left| {\alpha ,{\sigma ^2}} \right.} \right)$ ，得：

$ p\left( {{ t}\left| {\alpha ,{\sigma ^2}} \right.} \right) = {\left( {2{\text{π}} } \right)^{ - \frac{N}{2}}}{\left| {\varSigma } \right|^{ - \frac{1}{2}}}\exp \left( { - \frac{1}{2}{{ t}^{\rm T}}\varSigma ^{ - 1}{ t}} \right) $

(9)

分别对 $\alpha $ 和 ${\sigma ^2}$ 求偏导并令其等于0，得：

$\left\{ {\begin{aligned} & {\alpha _i^{\rm new} = \frac{{{\gamma _i}}}{{\mu _i^2}}}, \\ & {{{\left( {{\sigma ^2}} \right)}^{\rm new}} = \frac{{{{\left\| {t - { \varPhi} \mu } \right\|}^2}}}{{N - \displaystyle \sum\limits_{i = 0}^N {{\gamma _i}} }}} \end{aligned}} \right.$

(10)

式中， ${\gamma _i} = 1 - \varSigma _{i,i} $ ，通过迭代计算 $\alpha _i^{\rm new}$ 和 ${\left( {{\sigma ^2}} \right)^{\rm new}}$ 直到满足收敛要求为止。

1.2 相关向量机的组合核

相关向量机通过引入核函数实现样本数据从原始空间到特征空间的复杂映射过程，从而避免了高维空间中的维数灾难问题，核函数的合理与否直接影响着相关向量机的预测性能。目前，用于相关向量机的核函数主要分为局部核函数和全局核函数两类：局部核函数学习能力强，泛化能力弱；全局核函数泛化能力强，学习能力弱^[8–9]。为了获得一个学习能力和泛化推广能力均比较强的预测模型，需要考虑设计一种新的组合核函数。

常见的局部核函数包括高斯径向基核、逆多元二次核等，常见的全局核函数包括样条核、多项式核等。由于高斯核在众多工程应用中表现优异，应用广泛，在此以高斯核为基础，结合全局核中的多项式核，构造组合核相关向量机：

$\left\{ \!{\begin{aligned} & {{K_{\rm mix}}\left( {{x_i},{x_j}} \right) \!=\! \lambda {K_{\rm RBF}}\left( {{x_i},{x_j}} \right) \!+\! \left( {1\! -\! \lambda } \right){K_{\rm poly}}\left( {{x_i},{x_j}} \right)},\\ &\qquad {{\rm s.t.}\;\;0 \le \lambda \le 1,\;\;\;\lambda \in R};\\ & {{K_{\rm RBF}}\left( {{x_i},{x_j}} \right) = \exp \left( { - \frac{{{{\left\| {{x_i} - {x_j}} \right\|}^2}}}{{2{\sigma ^2}}}} \right)};\\ & {{K_{\rm poly}}\left( {{x_i},{x_j}} \right) = {{\left[ {\left( {{x_i} \cdot {x_j}} \right) + 1} \right]}^d}} \end{aligned}} \right.$

(11)

式中： ${x_i}$ 和 ${x_j}$ 为原始数据空间中的样本点； ${K_{\rm RBF}}\left( {{x_i},{x_j}} \right)$ 和 ${K_{\rm poly}}\left( {{x_i},{x_j}} \right)$ 分别为径向基核函数和多项式核函数； $\sigma $ 为核函数的参数； $\lambda $ 为核函数的权重，当 $\lambda = 0$ 或 $\lambda = 1$ 时， ${K_{\rm mix}}\left( {{x_i},{x_j}} \right)$ 退化为单一核函数，通过调节权重 $\lambda $ ，使组合核函数 ${K_{\rm mix}}\left( {{x_i},{x_j}} \right)$ 适应了不同的样本输入，兼具了核函数的学习能力和泛化推广能力。

2 运行状态评估的归一化EEMD信息熵

服役条件下的机械设备，其性能退化过程可看作是一个不确定性系统。信息熵是一种用来反映不确定性系统运行状态的量化指标，机械设备运行过程中所包含的运行状态信息可以用熵来度量^[10]。运行状态评估的归一化EEMD信息熵，通过采用集成经验模式分解处理机械设备运行时的振动信号，获得机械设备全寿命周期上不同时间节点的EEMD信息熵，并将其归一化为[0, 1]区间的数值，以此作为定义在[0, 1]区间的运行状态指标量值评估机械设备的服役性能^[10]。

2.1 EEMD信息熵原理

在机械设备运行过程中，所获取的振动信号序列为 $x\left( t \right)$ ，为了减少原始振动信号 $x\left( t \right)$ 在EMD分解过程中的模式混叠，EEMD通过在信号 $x\left( t \right)$ 中多次加入频率均匀分布的高斯白噪声，使得不同时间尺度上的特征信号 $s\left( t \right)$ 自动分布到参考尺度上，具体实现过程如下：

步骤1　在原始信号 $x\left( t \right)$ 中加入一个幅值系数为 $k$ 的高斯白噪声 $n\left( t \right)$ ，得到加噪后的信号 $X\left( t \right)$ ：

$X\left( t \right) = x\left( t \right) + n\left( t \right)$

(12)

步骤2　对加噪后的信号 $X\left( t \right)$ 进行EMD分解，得到 $n$ 个IMF分量和余项 $r\left( t \right)$ ：

$X\left( t \right) = \sum\limits_{j = 1}^n {{c_j}\left( t \right)} + r\left( t \right)$

(13)

步骤3　给原始信号 $x\left( t \right)$ 加入不同的高斯白噪声 ${n_i}\left( t \right)$ ，其中 $i = 1,2, \cdots ,M,$ 得到不同的含噪信号 ${X_i}\left( t \right)$ ，对含噪信号 ${X_i}\left( t \right)$ 分别进行EMD分解，得到 $n$ 个IMF分量 ${c_{ij}}\left( {j = 1,2, \cdots ,n} \right)$ ，其中 ${c_{ij}}$ 表示含噪信号 ${X_i}\left( t \right)$ 的第 $j$ 个IMF分量。

由于所添加的高斯白噪声 ${n_i}\left( t \right)$ 各不相同，利用不相关随机序列统计均值为0的特性，将上述对应的IMF进行集成平均，相互抵消高斯白噪声对真实IMF的影响，最终得到EEMD分解后的IMF为：

${c_j}\left( t \right) = \frac{1}{M}\sum\limits_{i = 1}^M {{c_{ij}}\left( t \right)} $

(14)

式中， ${c_j}\left( t \right)$ 即为对原始信号进行EEMD分解得到的第 $j$ 个IMF， $j = 1,2, \cdots ,n$ 。

2.2 归一化EEMD信息熵的运行状态评估指标量

通过对机械设备振动信号 $x\left( t \right)$ 的EEMD分解可以得到 $n$ 个IMF，第 $j$ 个IMF的频带能量 ${E_j}$ 和相对能量 ${\hat E_j}$ 分别为：

$\left\{ {\begin{aligned} & {{E_j} = {{\int_0^{ + \infty } {\left| {{c_j}\left( t \right)} \right|} }^2}{\rm d}t} ,\\ & {{{\hat E}_j} = {{{E_j}} \Bigg/ {\left( {\sum\limits_{j = 1}^n {{E_j}} } \right)}}} \end{aligned}} \right.$

(15)

式中， $j = 1,2, \cdots ,n$ ，相对能量 ${\hat E_j}$ 即为归一化的EEMD能量，且满足 $\displaystyle \sum\limits_{j = 1}^n {{{\hat E}_j}} = 1$ 。由此可得归一化的EEMD能量熵定义：

${w_{\rm ee}} = - \sum\limits_{j = 1}^n {{{\hat E}_j}\lg } {{\hat E} _j}$

(16)

由于本征模态分量 ${c_j}\left( t \right)\;\;\left( {j = 1,2, \cdots ,n} \right)$ 代表原始信号 $x\left( t \right)$ 从高频到低频不同的频带成分，每个本征模态分量所包含的频率成分及相对频带能量 ${\hat E_j}$ 均不相同。当 ${\hat E_j} = 0$ 时，则 ${w_{\rm ee}} = 0$ ，表明原始信号能量集中在某个或某几个频带内，此时机械设备运行状态单一确定；当各个频带能量均相等时，则 ${w_{\rm ee}}$ 最大，此时运行状态不确定性也最大。

当机械设备出现性能退化时，在相应的频带内就会出现相应的幅值波动，此时频带能量也会出现波动，使能量分布的不确定性减少，熵值减少，机械设备的服役性能降低，其对应的运行状态指标量也减小。由于归一化EEMD能量熵定义在[0, 1]区间，在得到归一化EEMD能量熵后，即可以此表征机械设备运行状态的优劣进而评估其服役性能^[10]。

2.3 归一化EEMD信息熵的轴承性能评估实例分析

采用Cincinnati大学实测的滚动轴承全寿命振动信号数据，选取某一轴承的振动信号数据进行分析，整个服役周期内，每间隔10 min进行1次振动信号采集，共采集该轴承的984个振动信号序列数据。首先对获取的振动信号序列数据分别进行EEMD分解提取各个IMF分量，以各个IMF分量和Shannon信息熵为基础计算归一化EEMD能量熵，整个服役周期内的归一化EEMD能量熵序列如图1所示。

从图1可以看出，整个服役周期内，归一化EEMD能量熵序列大致经历了5个阶段，如表1所示。表1中，第1列为轴承归一化EEMD能量熵的5个阶段变化描述，第2列为对应的时间序列区间，第3列为在该区间内某一时间序列点的归一化EEMD熵值，因为篇幅限制只列举了该区间内的代表点。由图1和表1可以看出，在熵值平稳阶段，即时间序列0～519点区间内，熵值几乎没有变化，说明在这一时间区间内，轴承性能良好。从520点到轴承最终故障失效，轴承的归一化EEMD熵值从开始出现缓慢变化到后来出现剧烈变化，这一区间对应着轴承性能已经开始退化，运行状态服役性能开始降低。由表1中第3列还可以看出，归一化EEMD能量熵值和运行状态指标量值均定义在[0, 1]区间并且变化趋势相同，如：在熵值平稳阶段代表点的归一化EEMD熵值为0.627 9，表明轴承的运行状态良好，性能指标量值较高；在熵值下降阶段代表点的归一化EEMD熵值为0.247 6，表明轴承的运行状态性能已经退化，性能指标量值降低。因此可以用轴承的归一化EEMD能量熵值来表征轴承的运行状态变化进而评估其服役性能。

3 机械设备运行状态评估及预测模型

归一化EEMD信息熵与组合核相关向量机的机械设备运行状态评估及性能退化预测流程如图2所示，主要由振动信号采集、信号预处理、计算归一化EEMD信息熵、构造输入向量集、RVM训练和预测等步骤组成。

机械设备的运行状态性能退化趋势预测模型通过平均相对误差、平均绝对误差以及预测值与真实值之间的相关系数来评价其预测性能，即：

平均相对误差 ${E_{\rm MAPE}} = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {\left| {\dfrac{{{y_i} - y'_i}}{{{y_i}}}} \right|} \times 100{\text{%}}$ ，

平均绝对误差 $MS\!E = \dfrac{1}{n}\sqrt {\displaystyle\sum\limits_{i = 1}^n {{{\left( {{y_i} - y'_i} \right)}^2}} } $ ，

相关系数 $R = \dfrac{{\displaystyle\sum {\left( {{y_i} - {{\bar y}_i}} \right)\left( {{y'}_i - {{\bar y}'}_i} \right)} }}{{\sqrt {\displaystyle\sum {{{\left( {{y_i} - {{\bar y}_i}} \right)}^2}\displaystyle\sum {{{\left( {{y'}_i - {{\bar y}'}_i} \right)}^2}} } } }}$ 。

式中， ${y_i}$ 和 ${ y_i^\prime}$ 分别为第 $i$ 个样本的预测值和真实值。

3.1 输入向量集的确定

在构建机械设备性能退化指标预测模型时，采用的振动特征参量为归一化EEMD能量熵，以机械设备前几个时段的归一化EEMD能量熵作为参数输入。考虑步长系数 $p$ ，即前 $p$ 个时刻对当前时刻机械设备服役性能的影响，设待预测设备为 $i$ ， ${w_{\rm ee}}\left( {i,t - 1} \right)$ 为待预测设备 $i$ 在 $t - 1$ 时刻的归一化EEMD能量熵值，则设备 $i$ 在 $t$ 时刻的归一化EEMD能量熵值输入/输出向量如表2所示。

依据表2，设以1点到700点之间的数据来建立模型，则初次训练模型时，总共需要采用前700个数据样本，总的训练次数为 ${\rm{700}} - p$ 。

3.2 PSO优化组合核RVM

将建立的输入向量输入到组合核RVM构建机械设备服役性能的预测模型，RVM中核函数由径向基核函数 ${K_{\rm RBF}}\left( {{x_i},{x_j}} \right)$ 和多项式核函数 ${K_{\rm poly}}\left( {{x_i},{x_j}} \right)$ 组成， $\lambda $ 为核函数的权重参数， $\sigma $ 、 $d $ 分别为径向基核和多项式核的核函数参数。

为了提高组合核RVM的预测精度，采用粒子群算法（particle swarm optimization，PSO）对组合核的参数进行自动寻优。较之其他优化算法，PSO算法结构简单且寻优速度较快，非常适合于较少参数的寻优过程。PSO算法将待优化参数初始化一群随机粒子 ${{X}} = \left[ {{x_1},{x_2}, \cdots ,{x_m}} \right]$ ，设第 $i$ 个粒子表示为 ${x_i} = \left[ {{\lambda _i},{\sigma _i},{d_i}} \right]$ ，然后通过迭代找到最优解，其算法描述如下：

$\left\{\!\!\!\!{\begin{array}{*{20}{l}} \begin{gathered} {v_i}\left( {t + 1} \right) = w\;{v_i}\left( t \right) + {c_1}\;{r_1}\left( t \right)\;\left( {{p_i}\left( t \right) - {x_i}\left( t \right)} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{c_2}\;{r_2}\left( t \right)\;\left( {{p_{\rm g}}\left( t \right) - {x_i}\left( t \right)} \right) ,\\ \end{gathered} \\ {{x_i}\left( {t + 1} \right) = {x_i}\left( t \right) + {v_i}\left( {t + 1} \right)\;} \end{array}} \right.$

(17)

式中， $t$ 表示进化代数， $w$ 表示惯性权重， ${c_1}$ 和 ${c_2}$ 分别为加速因子， ${r_1}$ 和 ${r_2}$ 均为在 $[0,1]$ 之间变化的随机数。每次迭代中，粒子 ${x_i}$ 通过跟踪两个极值 ${p_{\rm i}}$ 和 ${p_g}$ 来更新自己，其中 ${p_{\rm i}}$ 为粒子 ${x_i}$ 的历史最优解， ${p_{\rm g}}$ 为种群X中所有粒子的最优解， ${x_i}\left( t \right)$ 和 ${v_i}\left( t \right)$ 分别为粒子 $i$ 在 $t$ 时刻的位置和速度。在核参数优化过程中，种群数一般取 $[20,40]$ ，最大迭代次数为100，目标函数为混合核RVM预测结果的平均绝对误差 ${f_i}$ ：

${f_i} = \frac{1}{n}\sqrt {{{\sum\limits_{k = 1}^n {\left( {{y_{ik}} - {{\hat y}_{ik}}} \right)^2} }}} $

(18)

式中， ${y_{ik}}$ 和 ${\hat y_{ik}}$ 分别为第 $i$ 个粒子下第 $k$ 个样本的预测值和实际值， $n$ 为待预测的样本总数。

4 实验验证 4.1 数据来源

实验数据为Cincinnati大学实测的滚动轴承全寿命实验数据，实验装置的整体结构及传感器的安装位置如图3所示。在实验台的转轴上安装4个轴承，轴承型号为Rexnord公司的ZA–2115双列滚子轴承，交流电机通过带传动以2 000 r/min带动转轴旋转，实验时在待测轴承上施加6 000磅的径向载荷。通过PCB353B33加速度传感器每间隔10 min采集1次待测轴承X和Y方向的振动信号，直至轴承运行失效，采样频率为20 kHz，采样长度为20 480。以轴承1的振动信号序列为例来验证本文方法的有效性。

对实验台架中轴承1的全寿命振动信号数据，计算归一化EEMD信息熵的时间序列如图1所示，由图1可以看出，轴承的归一化EEMD信息熵序列在0到519点几乎没有变化，说明在这一时间段内轴承性能趋于稳定。从520点到轴承最终运行失效，对应的归一化EEMD信息熵序列经历了4个阶段：熵值下降阶段、熵值上升阶段、熵值再下降阶段、熵值随机阶段。

针对轴承1的全寿命振动信号数据，计算归一化EEMD信息熵并以此来表示轴承的运行状态。在全寿命的熵值变化阶段中，从熵值平稳到第1个熵值下降阶段之间非常关键，它是表征轴承性能开始退化的初始敏感信号，为了后续能够清晰地展示预测模型的预测效果，在此只显示600点到700点之间的数据，如图4所示。

4.2 优化模型参数及确定迭代步长

在此，分别构建组合核RVM和单一核RVM模型。在单一核RVM模型中，既采用高斯核函数也采用多项式核函数。组合核RVM和单一核RVM进行参数优化时，PSO算法的参数设置均相同。在每次训练阶段，考虑前3个时刻值对当前时刻滚动轴承服役性能的关联，即步长系数 $p = 3$ 。适应度函数采用平均绝对误差函数 $MSE$ ，通过粒子群算法分别优化单一核RVM（高斯核函数/多项式核函数）和组合核RVM中的参数，其结果如表3所示。图5为组合核RVM参数优化的适应度变化曲线。

在模型的训练过程中，通过迭代多步训练能够准确构建滚动轴承运行状态性能退化的预测模型，即每次只向前训练1步，每步都添加真实值来训练下一步。基于表3中优化的RVM参数值，采用组合核RVM模型，不断调整迭代步长 $p\left( {p = 2,3,4,5,6,7} \right)$ 值，进行样本训练和预测，滚动轴承在步长为3和5的迭代多步预测结果如图6所示。

优化组合核RVM在不同步长下预测结果的平均绝对百分比误差 ${E_{\rm MAPE}}$ 和样本的训练时间 $t$ 见表4。

从表4可以看出， ${E_{\rm MAPE}}$ 随着迭代步长 $p$ 的增大先变小趋于稳定后逐渐增大，预测模型的训练时间 $t$ 在不同步长系数 $p$ 下变化比较随机，综合考虑预测模型的训练效率和预测精度，在此确定训练时的步长系数 $p = 4$ 。

4.3 模型的预测性能比较

计算表征轴承1服役性能的归一化EEMD信息熵序列，采用单一核RVM和组合核RVM对信息熵序列数据进行预测，在此显示轴承1在600点到700点之间的振动信号序列数据，各种模型的实际预测结果如图7所示。

由图7可以看出，组合核RVM的预测值与实际值的贴合度较之单一核RVM方法更好。采用表2来定量评估不同预测模型的预测精度，其性能对比如表5所示。由表5可以看出，组合核RVM模型的 ${E_{\rm MAPE}}$ 和 $MSE$ 均优于单一的多项式核和径向基核RVM，预测精度最高。

RVM预测模型的平均绝对误差如图8所示，从图8可以明显看出，组合核RVM的预测模型对于缓变点预测的平均绝对误差都较小，但是对于第700点—突变点的预测有一定的滞后，因此误差有一定增大，但这并没有影响对轴承整体服役性能变化趋势的预估计。

在计算机配置环境为8G RAM和3.2 GHz Intel CPU下进行测试，用300个样本训练多核RVM，预测200个样本所需时间为2.56 s，测试结果表明组合核RVM具有快速预测能力，可以满足滚动轴承性能退化预测的快速性要求。

5 结　论

针对现役机械设备历史状态样本缺乏，运行状态难以快速评估，性能退化趋势难以有效预测的问题，提出基于组合核函数相关向量机与归一化EEMD信息熵的机械设备运行状态评估及性能退化趋势预测方法，得出如下结论：

1）从机械设备振动信号的相对能量出发，归一化EEMD信息熵可将机械设备的运行状态定义到[0, 1]之间，能够替代概率统计方法并正确评估机械设备的运行状态。

2）组合核RVM利用粒子群算法，依据机械设备运行状态性能指标参数的特性，自适应地调整RVM的核函数参数，实现了多特征空间的有效融合，提高了机械设备运行状态预测的精度和鲁棒性。

3）将归一化EEMD信息熵和组合核RVM应用于滚动轴承的运行状态评估及性能退化趋势预测，从实例上验证了所提方法的有效性和实用性。