基于二元混合随机过程的轴承剩余寿命预测 ①

金晓航 李建华 郭远晶 贾 虹*

(*特种装备制造与先进加工技术教育部重点实验室(浙江工业大学) 杭州 310023)(**浙江工业大学机械工程学院 杭州 310023)(***宁海县浙工大科学技术研究院 宁海 315600)(****浙江工业大学之江学院 绍兴 312030)

0 引 言

装备(如盾构机、风力发电机、水泵水轮发电机组等)在当今社会建设和经济发展中发挥了重要作用，装备一旦发生故障会造成重大事故和巨额损失[1-5]。为避免灾难事故的发生，同时降低运维成本，需要对装备健康状态进行有效的监测、评估并准确地预测其剩余寿命(remaining useful life， RUL)。目前用于剩余寿命预测的方法有贝叶斯滤波算法、神经网络算法、随机过程模型等[6-7]。基于随机过程的装备剩余寿命预测方法能够较好地用于分析装备在工作过程中受到外围环境、运行条件、制造工艺等不同因素造成的不确定性影响，且具有良好的计算分析特性，成为近年来研究的热点[8-9]。

为了更加全面地反映装备的健康状态和准确地预测其剩余寿命，通常需要同时利用多个性能指标。在利用多个性能指标共同表征装备的健康状态时，需考虑不同性能指标间的相关性问题[10]。如Bai等人[11]在对应力-强度模型的可靠性分析中利用Gumbel Copula函数描述应力和强度变量之间的相关关系，并通过工程案例进行分析，证明了所建立模型的实用性。金晓航等人[12]利用Copula函数分析两个性能指标间的相关关系，提出了一种基于二元维纳过程的轴承剩余寿命预测方法。上述的工作都是基于多个性能指标服从同一种随机过程的假设下开展的。但在一些情况下，所选取的性能指标性质差异较大时，不同指标无法同时服从同一种随机过程。

基于上述分析，考虑一种将维纳过程和伽马过程融合分析的二元混合随机过程，用于预测装备存在两个性质不同的性能指标时的剩余寿命。首先基于维纳过程和伽马过程对不同性能指标的退化过程分别进行建模分析，然后利用赤池信息准则(Akaike information criterion，AIC)选取合适的Copula函数分析不同性能指标间的相关性，采用分步极大似然估计法在线更新模型参数，实现装备的剩余寿命预测。

1 基于二元混合随机过程的退化建模

目前应用较为广泛的随机过程模型有维纳过程和伽马过程等[9]。维纳过程适合分析退化过程连续且呈现随机波动的性能指标；伽马过程适合分析退化过程严格单调递增(或递减)的性能指标[13-14]。因此，需要结合装备性能指标变化的特点选择合适的随机过程建立相应的退化模型。本文针对装备两个不同特点的性能指标，分别利用维纳过程和伽马过程进行建模分析。

1.1 基于维纳过程的退化模型

若装备退化过程的性能指标表现出非严格单调递增(或递减)的过程，可以利用维纳过程进行建模分析。模型的表达式如式(1)所示[12]。

X(t)=X(0)+μt+σB(t)

(1)

式中，X(t)表示为t时刻的性能指标值，μ为漂移系数，σ为扩散系数，B(t)为标准维纳过程。

利用维纳过程分析的退化过程{X(t),t≥0}具有如下性质。

(1)X(0)=0；

(2) {X(t),t≥0}具有平稳独立增量；

(3) 任意时刻间的增量服从正态分布，即：

X(t+Δt)-X(t)～N(μΔt,σ2Δt)

根据维纳过程的定义，增量ΔX(t)服从正态分布，令Δt=1，其概率密度函数可表示为

(2)

1.2 基于伽马过程的退化模型

若装备退化过程的性能指标表现为严格单调递增(或递减)的情况时，可以利用伽马过程进行建模分析[15]。设{G(t),t≥0}是一个形状参数为α，尺度参数为β的伽马过程，利用伽马过程分析的退化过程应满足以下性质。

(1)G(0)=0；

(2) {G(t),t≥0}具有平稳独立增量；

(3) 任意时刻间的增量服从伽马分布，即：

G(t+Δt)-G(t)～Ga(αΔt,β)

根据伽马过程的定义，增量ΔG(t)服从伽马分布Ga(αΔt,β)，令Δt=1，其概率密度函数为

(3)

1.3 相关性分析

利用多个性能指标同时表征装备的健康状态时，需要考虑不同性能指标之间的相关特性。Copula函数能够有效地构建出各个性能指标的联合分布函数和边缘分布函数之间的关系，可被用于分析不同性能指标之间的相关性。因此，引入Copula函数分析装备多个性能指标之间的相关关系，如下式所示。

F(x1,x2,…,xn)=C(u1,u2,…,un;θ)

(4)

式中，F(x1,x2,…,xn)为多个变量的联合分布函数；C(u1,u2,…,un;θ)和θ是分别为Copula函数和Copula函数中的参数；u1,u2,…,un为各个变量的边缘分布函数。

当各变量的边缘分布函数已知时，它们的联合分布函数可根据Sklar定理[15]获得。用Copula函数描述两个退化性能指标之间的相关特性时，联合分布函数可表示为

F(x1,x2)=C(F(x1),F(x2);θ)

(5)

考虑到常用的Copula函数有Gaussian Copula、Frank Copula、Gumbel Copula和Clayton Copula等[16]，因此需要结合装备实际退化过程选择合适的Copula函数。鉴于赤池信息准则(AIC)是一种广泛被用于评价模型分析效果优劣的方法[17]，文章中利用AIC信息准则选择合适Copula函数：

AIC=-2log(A)+2m

(6)

式中，A为模型对应的极大似然函数，m为模型中参数个数，AIC的值越小说明所选择的模型越合适。

2 剩余寿命预测

由于同时利用了装备的两个性能指标(性能指标1服从维纳过程，性能指标2服从伽马过程)进行剩余寿命预测，所以定义当{X1(t),t≥0}和{X2(t),t≥0}中的任意一个性能指标超过其所对应的失效阈值ω1、ω2时，即认为装备发生故障[18]。因此，当利用两个性能指标预测分析装备的剩余寿命时，剩余寿命T定义为

T=inf{t：X1(t)≥ω1或X2(t)≥ω2}

(7)

2.1 剩余寿命分布

装备退化时，若其性能指标变化的情况服从维纳过程，则装备的剩余寿命服从逆高斯分布[19]。确定维纳过程模型中的参数μ、σ后，即可得到剩余寿命TW的概率密度函数：

fTW(t|ωW,μ,σ)=

(8)

若装备的性能指标变化服从伽马过程时(其性能指标初始值为0，失效阈值为ωG)，根据伽马过程的性质可知，装备退化过程中其性能指标到达失效阈值的首达时间TG为

FTG(t)=P{TG≤t}=P{G(t)≥ωG}

(9)

通过式(9)，TG的概率分布函数可以表示为

(10)

式中，Γ(a,b)为不完全伽马函数，可表示为

(11)

装备剩余寿命TG的概率密度函数为

(12)

由于通过式(12)直接计算fTG(t|ωG,α,β)比较困难，可以通过离散化的方式获得装备剩余寿命概率密度函数的近似表示。取tn=n/α,n=0,1,…，并令qn=P{tn

=αqn

(13)

基于不同性能指标获得装备剩余寿命概率密度函数后，利用Copula函数分析性能指标间的相关特性，即可获得装备剩余寿命T的联合概率函数：

fT(t|ωW,ωG,μ,σ,α,β,θ)=c(FTW(t),FTG(t);θ)

·fTW(t|ωW,μ,σ)

·fTG(t|ωG,α,β)

(14)

式中c(·)为C(·)的密度函数。

2.2 参数估计

由式(14)可知剩余寿命预测模型中包含5个未知参数(μ,σ,α,β,θ)，可通过分步极大似然估计法求解各个参数：(1)更新不同性能指标退化模型中的参数(μ,σ,α,β)；(2)更新Copula函数中的参数θ，具体过程如下。

第1步。由维纳过程的性质知道，服从维纳过程的性能指标增量ΔXi～N(μΔti,σ2Δti)，可得到模型参数(μ,σ)的似然函数为

(15)

基于式(15)，分别对漂移系数μ和扩散系数σ求偏导，可得μ、σ的极大似然估计值分别为

(16)

(17)

由伽马过程的性质知道ΔGi～Ga(αΔti,β)，可获得模型参数(α,β)的似然函数为

L(α,β)=

(18)

由极大似然估计法，令

(19)

(20)

3 仿真算例和实例分析

3.1 仿真算例分析

通过仿真算例来验证分析本文所提出的基于二元混合随机过程的剩余寿命预测方法的可行性。分别利用维纳过程和伽马过程的性质(维纳过程的增量服从正态分布，伽马过程的增量服从伽马分布)，在Matlab软件中生成两个性能指标数据，失效阈值分别定义为(ω1,ω2)=(55.18, 31.89)，如图1所示。两个性能指标的增量分别满足正态分布N(0.15, 0.3)和伽马分布Ga(0.2, 0.4)，如图2所示。

图1 两模拟性能指标退化趋势

图2 两个模拟退化增量

利用AIC信息准则选择Copula函数，结果如表1所示。从表1中可以看到，针对该组数据，Frank Copula函数所对应的AIC值最小，因此选择Frank Copula函数分析性能指标间的相关特性。利用式(8)和(13)分别计算基于不同性能指标的装备剩余寿命的边缘概率密度函数。

表1 4种Copula函数的AIC值

根据Sklar定理构建出剩余寿命的联合概率密度函数，利用两步极大似然估计法更新模型参数。首先，基于两个性能指标各个时刻的测量值，结合极大似然估计法在线更新退化模型中的参数(μ,σ,α,β)，进而确定两个不同性能指标剩余寿命的边缘分布函数。然后，将两个边缘分布函数作为Copula函数的输入，再次利用极大似然估计法更新Copula函数中的参数θ，获得各个时刻参数的估计结果，如图3所示。

图3 各个监测时刻点的参数估计值

获得模型参数的估计值后，结合式(14)计算得到装备的剩余寿命概率密度函数，如图4所示。观察不同时刻的概率密度函数可以看出，随着观测数据的增多，剩余寿命预测结果的分布范围逐渐缩小，说明其预测结果的精确度越来越高。将各个时刻概率密度函数最大值所对应的时间作为剩余寿命的预测值，可得到整个过程的剩余寿命预测的结果，如图5所示，从图中可以看出剩余寿命的预测值逐渐逼近真实值。

图4 不同时刻点剩余寿命的概率密度函数

利用均方根误差(root mean square error, RMSE)对不同方法预测结果的误差进行分析，结果如表2所示。RMSE值越小说明预测效果越好[20]。从表2中可以看出，相较于基于一元随机过程的剩余寿命预测方法，基于二元混合随机过程的剩余寿命预测方法的预测效果更好。

图5 剩余寿命预测结果

表2 不同预测方法预测结果的RMSE

3.2 实例分析

利用PRONOSTIA实验平台上的轴承全寿命周期的实验数据对基于二元混合随机过程的剩余寿命预测方法做进一步的验证分析[21]。实验过程中，利用Dytran3035B的加速度传感器采集振动信号，通过NI数据采集系统每隔10 s，以25.6 kHz的采样频率，采集并存储一个时长为0.1 s的振动数据；利用Pt100的温度传感器连续采集温度信号，采样频率为10 Hz。利用轴承振动信号的有效值和温度信号的平均值作为表征其健康状态的性能指标。图6(a)、(b)是轴承性能退化阶段，两个性能指标归一化后的变化过程，分别利用维纳和伽马过程分析振动和温度的数据。对轴承两个性能指标建模分析之前，需要判断它们是否符合随机过程的性质，通过Jarque-Bera验证分析，振动信号有效值的增量服从正态分布N(0.0045, 0.0483)，如图7(a)所示。由于受环境温度以及实验平台其他部件传热等因素的影响，轴承温度信号的实际增量并不一定符合伽马分布，因此需要对温度信号进行预处理分析(即对时序的温度信号进行包络分析，且令后一时刻的温度值大于等于前一时刻的值)。结合伽马过程的性质，预处理后的温度性能指标如图6(c)所示，经检验，该信号的增量服从Ga(1.2645, 0.0015)的伽马分布，如图7(b)所示。图7的结果表明针对轴承所构建的两个性能指标可以分别利用维纳和伽马过程进行建模分析。

图6 轴承两性能指标退化趋势

图7 轴承两个性能指标增量的直方图

通过AIC信息准则选择Copula函数分析不同性能指标间的相关特性，如表1所示。从表1中可知Frank Copula函数所对应的AIC值为-3 856，相较于其他3个Copula函数，Frank Copula函数的AIC值最小，因此选择Frank Copula函数分析轴承两个性能指标间的相关特性。进而构建轴承剩余寿命的联合概率密度函数，利用分步极大似然估计法更新模型中的参数，在轴承性能退化阶段的第100、120、140、160、180时刻的模型参数估计值，如表3所示。

表3 不同监测时刻点所对应的参数估计值

获得各个时刻模型参数的估计值后，结合式(14)即可获得对应时刻的剩余寿命概率密度函数，如图8所示。图中实线表示的是概率密度函数最大值对应的时刻，即为所预测的剩余寿命；虚线表示的是当前时刻剩余寿命的真实值，实线与虚线之间的距离代表预测结果的误差。图8中(a)到(d)依次为轴承退化阶段第100、130、150、170时刻的剩余寿命概率密度函数图，从图中可以看出，随着测量时刻点的增加，实线与虚线之间的距离不断减小，说明预测误差逐渐缩小、预测精度不断提高。

图8 不同时刻剩余寿命概率密度函数

将不同时刻剩余寿命概率密度函数的最大值对应的时间提取出来，获得如图9所示的轴承剩余寿命预测结果。从图中可以看出，相较于单独利用维纳过程或伽马过程的剩余寿命预测方法，基于二元混合随机过程的剩余寿命方法的预测效果更好。

为了客观地衡量不同剩余寿命方法的预测效果，利用均方根误差和误差区间分布对预测结果的精度分别进行了比较分析。不同方法预测结果的RMSE值如表2所示，从表中可以看到基于二元混合随机过程的轴承剩余寿命预测方法的RMSE值为

图9 剩余寿命预测结果

3.6535，小于其他两种方法的RMSE值。同时将所有预测结果的误差值从大到小划分为[-300, 30]、[-100, 20]、[-50, 10]和[-10, 5]4个不同的区间，结果如表4所示。从表中可以看出基于二元混合随机过程的剩余寿命预测方法在各个区间中所占的比例分别为100%、82.61%、54.35%、27.17%，比其他剩余寿命预测方法在各个不同误差区间中所占的比例更高。综合上述分析，说明基于二元混合随机过程的轴承剩余寿命预测方法的预测精度较高。

表4 不同剩余寿命预测方法的预测结果误差区间对比分析

4 结 论

针对装备存在多个性能指标且指标随时间变化的不同特点，提出利用维纳和伽马过程相融合的二元混合随机过程用于预测装备的剩余寿命。

利用AIC信息准则选择合适的Copula函数分析不同性能指标间的相关特性，通过分步参数估计方法更新模型参数并预测装备未来时刻的健康状态，进而实现剩余寿命的预测。

分别利用仿真和轴承实验数据对所提方法进行验证分析，结果表明基于二元混合随机过程的剩余寿命预测方法的预测效果比单独利用维纳过程或伽马过程进行建模分析的剩余寿命预测方法更好。