基于无迹卡尔曼滤波的离心泵剩余寿命研究*

王炳波 刘赓传

(中国石油工程建设有限公司北京设计分公司)

0 引 言

泵类设备广泛应用于石油石化行业，为液体或气体物料的输送提供动力。为保证系统的安全稳定运行，应对泵类设备性能退化和寿命预测进行研究[1]。现有的泵类设备寿命预测方法大致可以分为基于物理失效模型的方法和数据驱动的方法，基于数据驱动的方法又可分为机器学习方法、多元统计分析方法、特征量提取方法和信息融合方法等[2-3]。郭文琪等[4]综述了矿井主排水泵寿命预测方法，指出泵类寿命衰退指标应该多样化，只有多种预测模型有机结合才能保证预测模型的准确性；刘锐等[5]从离心泵流量公式出发，结合流量与间隙的物理关系模型，建立了离心泵寿命评估物理模型；韩可[6]基于变分模态分解(VMD)和支持向量数据描述(SVDD)方法，建立了液压泵性能退化模型；何庆飞等[7]建立了基于灰色支持向量机的液压泵寿命预测方法；马济乔等[8]搭建了液压泵加速退化试验平台，建立了液压泵寿命预测方法，并进行了预测方法的可靠性分析。基于物理模型的离心泵寿命预测方法，由于模型简化和失效机理的多样性，所以没有普遍适用性；基于数据驱动的泵类设备预测方法不能根据实时监测的数据状态对退化模型进行实时更新。因此，本文基于无迹卡尔曼滤波方法，提出了一种结合已有同类设备的退化数据，并有机结合特定离心泵实时退化数据的离心泵剩余寿命预测方法。

1 理论基础

1.1 离心泵性能退化模型的构建

离心泵的寿命与多种因素有关，如叶轮和轴承磨损、密封失效等，在性能退化期间其扬程、功率和效率等参数会发生改变，因此可以用扬程、功率和效率等参数来表征离心泵的健康状况。离心泵扬程H与流量Q的关系可以近似表示为如下的二次函数：

H=K1+K2Q-K3Q2

(1)

式中：K1、K2、K3为常数。

叶轮以及轴承磨损是离心泵常见的故障，磨损速率直接影响了离心泵的使用寿命[9]。叶轮磨损主要影响公式(1)中的常数项K1，轴承磨损主要影响泵转速n，上述磨损规律可用如下公式表示：

(2)

同时考虑两种磨损，即把公式(2)带入公式(1)，经过变换后可得[10]：

(3)

综上，离心泵扬程的寿命衰退关系可表示为：

H=K1+K2Q-K3Q2-At2-Bt

(4)

式中：A、B为常数。

由于离心泵的功率和效率与离心泵扬程成正比，所以其功率和效率的性能衰退关系也可表示为二次函数的形式。为了综合考虑多个性能衰退指标，引入离心泵寿命健康指标Hi。该指标根据熵值法综合考虑离心泵的扬程、功率和效率，因此离心泵运行时间t与健康指标的关系可用如下二次函数来表示：

Hi=a′t2+b′t+c′

(5)

式中：a′、b′、c′为常数，a′<0。

1.2 支持向量机基本原理

(6)

根据KKT条件[14]，上式可转化为凸二次优化问题：

(7)

当x与y为非线性关系时，可通过满足Mercer条件的核函数将该非线性问题转化为更高维度的线性问题，此时只需将式(7)中的xi·xj替换为核函数K(xi·xj)。常用的核函数为RBF核函数和Ploy核函数，本文采用RBF核函数(p为核函数的指数)，其表达式为：

(8)

1.3 卡尔曼滤波基本原理

卡尔曼滤波(Kalman filtering，KF)是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法，经过几十年的发展已拓展至非线性系统的状态评估，可应用在轴承寿命预测、交通流预测、电池剩余寿命预测和机械零件寿命预测等方面[15-18]。无迹卡尔曼滤波(Unscented Kalman Filter，UKF)是一种利用采样策略逼近非线性分布的滤波方法，其核心是无迹变换(Unscented Transform，UT)，即一种计算非线性变换中变量的统计特征的新方法，UKF既突破了卡尔曼滤波不能解决非线性问题的局限，又克服了拓展卡尔曼滤波(EKF)计算量较大的缺点。

卡尔曼滤波是一种特殊的贝叶斯滤波，是基于所有可用信息来构造系统状态的概率密度函数，即根据给定随机变量Zk=[z1……zk]T、Uk=[u0……uk]和初始状态x0的概率密度ρ0，估计出后验概率密度ρ(xk|Zk,Uk)，然后根据观测数据来计算状态量取不同值时的置信度，进而得到状态的最佳估计。其预测和更新原理为：假设k时刻状态向量的后验概率密度函数已知，则通过系统状态方程可得到k+1时刻系统状态向量的先验概率密度函数[19]：

(9)

当k+1时刻观测向量已知时，可根据贝叶斯准则用zk+1更新状态向量的先验概率密度函数：

(10)

式中：xk为系统的状态向量，uk为系统的输入向量。

2 预测模型的建立与求解

2.1 预测模型算法流程

基于UKF的离心泵寿命算法流程如图1所示。根据同类离心泵的退化历史数据，采用熵值法确定复合健康指标，并确定离心泵失效时的健康指标阈值，通过数据拟合得到离心泵卡尔曼滤波状态方程，并得到k时刻与k-1时刻状态方程的传递关系。对于需要预测寿命的离心泵，根据其历史运行数据和出厂设定的寿命，通过支持向量机回归得到剩余寿命与健康指标之间的非线性关系，将其作为卡尔曼滤波观测方程对原来的状态方程进行更新，最后把根据状态方程得到的先验估计剩余寿命和根据观测得到的剩余寿命加权(卡尔曼增益)，得到运行状态下离心泵的剩余寿命。

图1 基于UKF的离心泵寿命算法流程图Fig.1 Centrifugal pump life algorithm flowchart based on UKF

2.2 熵值法构建复合健康指标

根据寿命退化模型可知扬程H、功率P以及效率η均与离心泵的寿命有关，且这3个参数比较容易获取。采用多个指标更容易获得离心泵性能退化的真实情况，因此通过熵值法把这3个因素综合起来构建离心泵健康指标Hi。熵值法的原理是根据各个因素的离散程度来确定各因素的重要程度，熵值越大，该因素在所有因素中的权重就越大。构建的复合健康指标如下：

Hi=wHH+wPP+wηη

(11)

式中：w表示权重，其下标表示各个指标。

权重的计算方法为：

(12)

式中：ej为第j个指标的信息熵值，具体公式见文献[20]。

根据已有同类离心泵相关历史数据，可以较容易获得离心泵失效时的性能指标，该值即为健康指标的阈值。

根据该阈值和当前指标可以计算离心泵的剩余寿命。但该复合健康指标只有在符合高斯分布的情况下才能更好地应用卡尔曼滤波算法，因此如果构造的初始复合健康指标不能近似为高斯分布，则可通过Box-Cox变换将其转化为高斯分布[21]。

2.3 UKF滤波剩余寿命模型方程构建与求解

健康指标与离心泵运行时间之间的关系可表示为公式(5)所示的二次函数。当健康指标阈值确定后，对应的离心泵寿命设置为0，若已知离心泵的总寿命，则公式(5)可转化为健康指标和剩余寿命的函数关系。该二次函数的系数可通过同类离心泵的全寿命试验数据拟合得到。由于离心泵退化曲线函数光滑且单调(不考虑突发故障)，所以剩余寿命亦可表示为健康指标的二次函数，进而可以得到k时刻与k-1时刻离心泵剩余寿命与健康指标Hi之间的关系式：

Lk=Lk-1+A′·Hi+wk

(13)

公式(13)即为卡尔曼滤波的状态方程。式中的A′为常数，wk为符合正态分布的噪声，受不同厂家不同批次产品的影响。

如果要对某一特定的离心泵进行寿命预测，则除了得到上述由同类离心泵设备拟合得到的状态方程外，还要考虑该设备实时运行的参数。因此，还需建立观测方程来对状态方程进行实时更新，以获得运行中离心泵的剩余寿命。观测方程为一个非线性函数，该非线性关系通过支持向量机对离心泵的历史运行数据学习得到，具体过程如下：

对于样本序列S={(Hi(k),L(k))，Hi(k)∈Rl,L(k)∈R,k=1,2,……,N}，前n个样本数据用于支持向量机训练，后N-n个数据样本用于预测。构造如式(14)和式(15)所示的训练样本对Xtrain、Ytrain和预测样本输入Xtest。具体到本文，将收集的离心泵剩余寿命和健康指标数据分别作为输入X和输出Y，归一化和随机顺序后，挑选前n个数据代入式(7)，则可求解拉格朗日乘子，得到非线性关系h(·)，然后将拉格朗日乘子与剩余寿命预测样本输入回归函数，再计算得到预测样本，输出健康指标Ytest[22-23]。

Xtrain=[L(1),L(1),……,L(n)]

(14)

Ytrain=[Hi(1),Hi(2),……,Hi(n)]

(15)

Xtest=[L(n+1),L(n+2),……,L(N-n)]

(16)

上述训练过程得到的非线性函数即为卡尔曼滤波算法中的观测方程，可表示为：

Hi,k-1=h(Lk-1)+vk

(17)

式中：vk为SVM模型残差，假设为高斯白噪声，来源于测量仪器仪表的噪声。

2.3.1 初始向量(剩余寿命)

剩余寿命计算如下：

(18)

(19)

2.3.2 计算Sigma点

对2n+1个Sigma点进行如下计算：

(20)

2.3.3 状态更新

状态方程(公式(13))基于Sigma点下为：

χk|k-1=f(χk-1)

(21)

将Sigma点进行加权求和：

(22)

求Sigma点的加权预测方差：

(23)

式中：Qk为系统噪声向量wk的协方差矩阵。

进行Sigma点的观测方程(公式17)预测值求解。

yk|k-1=h(χk-1)

(24)

(25)

2.3.4 进行观测更新

协方差估计进行第k次更新：

(26)

观测中协方差进行第k次更新：

(27)

式中：Rk为测量噪声向量vk的协方差矩阵。

卡尔曼增益进行第k次预测更新：

(28)

最优状态估计进行第k次更新：

(29)

最优预测误差的估计进行第k次更新：

(30)

3 模型仿真与验证分析

3.1 同类离心泵退化模型数据拟合

用MATLAB实现上述模型算法，采用郭文琪等的离心泵加速退化试验数据对模型进行验证和分析。离心泵加速退化试验所用离心泵型号为IS100-80-125，扬程为20 m。选取石英砂作为磨料加入水中，石英砂颗粒尺寸为2～4目，避免颗粒较大引起管路堵塞，从而对离心泵的运行数据造成干扰，离心泵每天运行8 h，每2 h采集一次试验数据，试验时长25 d，共有200 h、100组寿命预测指标时间序列数据[7]。采用二次曲线的形式对原始数据拟合，拟合曲线如图2所示，方程如式(31)所示。

y=-0.000 203 8x2+0.001 095x+8.835

(31)

该拟合模型均方根误差为0.098 5，R-square系数为0.997 7，拟合效果较好。可见用二次函数建立离心泵剩余寿命和健康指标的方法较为可行，同时也验证了模型中卡尔曼状态方程的可行性，拟合噪声的分布可由寿命指标的真实值和拟合曲线之间的误差得到，并变换为高斯白噪声。

图2 寿命健康指标与时间的关系Fig.2 Relationship between life health index and time

3.2 基于SVM方法的观测方程可行性验证

上述拟合曲线只能反映类似机理同类型离心泵性能退化趋势，针对某一特定离心泵，除了考虑共性因素外，还需要考虑其个性因素，即对性能退化模型中的系数及时更新。对于某一处在退化期的同类但应用于不同工况的离心泵，假设已知其前期的80组退化数据，计算健康指标后，归一化处理并打乱数据顺序，选取前60组数据对SVM进行训练，采用后20组数据进行预测。模型的初始参数如下：误差初值为 0.000 1，惩罚因子初值为10 000，采用10-fold交叉验证度量SVM精度及推广能力。最终得到SVM方法预测的健康指标如图3所示。健康指标的平均绝对误差为-0.004 7，均方根误差为0.054 1。可以看出SVM方法对同类离心泵的历史数据回归和预测比较准确，可以用于模型中卡尔曼滤波观测方程的建立。

图3 寿命健康指标实际值与SVM预测值Fig.3 Actual value vs SVM predicted value of life health index

3.3 无迹卡尔曼滤波仿真结果

离心泵在衰退期内的性能退化过程可以利用上述建立的模型进行跟踪，跟踪结果如图4所示。退化初期由于可利用的测量数据较少，所以跟踪结果不理想，与真实值之间的误差较大；随着离心泵可利用数据的增加，UKF跟踪结果与真实值之间的误差慢慢变小。寿命健康指标跟踪误差如图5所示。由图5可以看出，当跟踪结果在58 h之后，所建立的方法能够较为准确地反映离心泵的性能退化情况。假设寿命健康指标<3时，离心泵失效，则进入衰退期后，离心泵可运行约178 h。上述卡尔曼滤波模型和试验方法得到的剩余寿命和衰退期内运行时间的关系如图6所示。由图6可以看出，当模型初值与真实值误差较大时，随着运行时间的进行，建立的卡尔曼滤波方法能够较快收敛到真实剩余寿命附近。但由于目前试验数据不足，所以在未来的工作中需要获取更多的试验数据来验证和完善建立模型的稳定性和预测精度。

图4 UKF方法离心泵衰退期跟踪结果Fig.4 Centrifugal pump degradation track by UKF method

图5 UKF方法寿命健康指标跟踪误差Fig.5 Life health index tracking error of UKF method

图6 剩余寿命与衰退期内运行时间的关系Fig.6 Relationship between residual life and running time in degradation period

4 结 论

根据理论分析，考虑离心泵的叶轮和轴承磨损，推导了离心泵扬程的衰退模型，并利用支持向量机和无迹卡尔曼滤波方法建立了离心泵剩余寿命预测模型，得出如下结论。

(1)考虑叶轮和轴承磨损失效时，推导了离心泵寿命与其扬程、功率及效率的关系，该关系可用二次函数的形式来表示。根据离心泵加速退化试验数据，并对试验数据进行拟合，拟合结果比较符合二次函数规律。

(2)根据在线监测的部分离心泵性能退化健康指标数据，支持向量机能够较为准确地预测离心泵性能退化过程，可以为卡尔曼滤波提供较为准确的观测值。

(3)MATLAB仿真结果表明，建立的模型在生命健康指标跟踪初期误差较大，但当跟踪一段时间后，随着可用数据的增加，跟踪误差逐渐减小，可以较为准确地描述离心泵性能退化过程，并对剩余寿命进行预测。