基于Gamma退化过程的剩余寿命预测及维修决策优化模型研究

谷玉波，贾云献，张英波

(军械工程学院 装备指挥与管理系，石家庄 050003)

机械设备在使用过程中，某些部件的运行状态会逐渐退化。状态的退化可以表现为机械部件的磨损、裂纹的增长、腐蚀程度的加深等，这些都是一系列理化效应的结果。产品状态退化的最终结果会导致功能故障的发生，而且在产品的退化过程中能够直接表征退化状态的状态指标(磨损、裂纹、腐蚀等)经常是不可测量或难于直接进行测量的。为获得设备状态退化程度的数据，可以运用CBM状态监控技术对设备的状态数据进行提取分析，间接地反映出设备的健康状态，实现对设备剩余寿命的预测，最终实现维修决策的优化，使维修效果达到最佳。在此，运用Gamma退化过程状态空间模型来描述设备状态性能的退化，通过该模型实现了剩余寿命的预测，并运用剩余寿命预测结果对维修决策模型进行优化。

1 Gamma退化过程的状态空间模型

在常用的分布中，Gamma分布可以用来描述连续累积磨损过程，并且具有退化过程建模所需的所有属性，即非负、平稳增长和独立增量过程[1]。因此，运用Gamma分布建立状态空间模型，可以描述设备的退化过程，并实现剩余寿命的预测。

1.1 状态空间模型

状态空间模型(State Space Model，SSM)是一种用来描述系统演化过程的模型，可以用来建立产品的状态退化模型。状态空间模型假设系统状态随时间的演化可以由一个不可观测或难以观测的状态序列确定，与状态序列伴随的是一个可以观测的测量序列，两者之间的关系可以由状态空间模型来描述，其中应用较为普遍的状态空间模型是典型相关(Canonical Correlation)方法[2]。通过建立状态方程和观测方程，状态空间模型为描述系统的动态演化特性建立了一个有效的模型框架，一些复杂的问题就可以用比较简单的形式表示。状态空间模型的基本方程为

xi+1=Fi(xi,ui,ζ,θ1)，

(1)

yi=Hi(xi,ui,ε,θ2)，

(2)

式中：Fi为状态转移函数，Hi为测量函数，一般而言Fi和Hi假设是已知的；ui为系统的输入变量或控制变量；θ1,θ2分别为状态方程和观测方程的静态参数；ζ和ε分别为系统噪声和观测噪声，一般假设ζ和ε相互独立，且服从某一特定分布。

1.2 Gamma过程

Gamma过程是一个具有独立、非负增量的随机过程。运用Gamma过程对一些退化过程建立模型，能够使建模过程中的数学计算更加明确易懂[3]，适合描述在时域上[4]单调累积微小增量的退化过程，比如磨损、腐蚀、裂纹增长等。

根据Gamma过程，定义非负随机过程{X(t);t≥0}，其概率密度函数为

(3)

2 Gamma-SSM退化模型

2.1 模型的建立

假设设备状态x的退化符合Gamma过程，y表示与x相对应的观测变量，xi表示设备ti时刻所对应状态，则建立的状态空间模型的状态方程和观测方程为

xi-xi-1～Gamma(α(ti)-α(ti-1),λ)，

(4)

yi=Hi(xi)+ε。

(5)

假设当系统状态x达到故障阈值Xf时系统发生故障，如图1所示。则系统状态从0首次达到故障阈值所需的时间为

图1 系统退化过程

Tf=inf{t:x=Xf,t>0}。

(6)

为简化问题，假设观测噪声是可以加的。假设设备状态退化过程为平稳Gamma过程，形状参数α(t)为线性函数，即α(t)=a·t，并且假设观测量y与状态量x之间的关系为y=c·x+ε，则在获得观测序列y0:n的条件下，退化模型可以表示为

xi+1-xi～Gamma(α·(ti+1-ti),λ)，

(7)

yi=c·xi+ε。

(8)

由假设可知，观测噪声ε服从均值为0的正态分布。所以，y～N(c·x,σ)。因此只要获得参数a,λ,c,σ的值，就可以确定该退化模型的形式。

2.2 模型参数估计

在对上述模型中参数进行求解的过程中，会出现新的需要估计的参数，导致运算量增大。为解决此问题，运用文献[5]提到的经验最大化(Experience Maximization，EM)算法与粒子滤波(Particle Filtering，PF)算法相结合的方法求解模型中出现的参数。

EM算法是进行极大似然估计的一种有效方法，简单来说，这一算法的每一次迭代主要包括2个步骤：(1)求期望(Expectation Step)，称为E步；(2)求极大值(Maximization)，称为M步。EM算法通过假设潜在变量的存在，来简化似然函数极大估计[3]。

PF算法是一种基于Monte Carlo方法和递推Bayes估计的统计滤波方法[6]，建立在序贯重要性采样(Sequence Importance Sampling，SIS)和Bayes理论的基础上，对于解决非线性、非Gauss问题有很好的效果[7]。PF算法的基本思想是：首先根据系统状态向量的经验分布在状态空间产生一组随机样本集合(即粒子)，然后根据观测数据不断调整粒子的权重和位置，通过调整后的粒子信息修正最初的经验分布。其实质是用由粒子及权重组成的离散随机测度近似相关的概率分布，并根据算法递推更新离散随机测度。当粒子容量足够大时，就近似于状态变量真实的后验概率分布。

由于EM-PF算法兼有两种算法各自的优点，可以更加有效地对模型中的参数进行求解。模型参数的求解过程如图2所示。

图2 模型参数求解过程

3 剩余寿命预测

剩余寿命是指产品从被检测的某一时刻起到该产品发生故障的时间长度[8]，设备剩余寿命的预测是维修管理中的一个重要环节，由于剩余寿命预测对产品使用过程中的安全性、经济性和任务性具有重要的影响，因此它是制定维修决策的重要依据。

通过确定状态空间模型参数，就能得到Gamma状态空间模型的具体形式，进而得到产品的剩余寿命分布函数。令Ns为粒子个数；Tf为设备故障时间；wi为第i个粒子的权重；Xf为故障阈值；y0:c为直至当前时刻所获得的观测序列；y0:i为当前获得的观测值序列；ti为第i次状态检测时间；xi为ti时刻的状态值；x1:n为n次检测所获得的状态序列，即x1:n={x1,x1,x2,…,xn} ；yi为ti时刻的退化量观测值；y0:n为n次检测所获得的观测序列，即y0:n={y0,y1,y2,…,yn},则ti时刻的剩余寿命τi分布函数为

(9)

p(τi+ti≥Tf|xi)=p(x(τi+ti)≥

(10)

(11)

另外，预测tk时刻的状态概率密度函数可表示为

(12)

(13)

(14)

(15)

4 维修决策优化模型

构建剩余寿命概率密度函数的最终目的是根据剩余寿命与实际需求做出维修决策，确定合理的维修或更换时间。平时对设备进行维修管理的过程中，人们总是希望设备的维修管理费用尽可能达到最低。以费用最小为目标的维修决策模型可以根据设备当前的运行状态确定在下一间隔期内是否进行维修以及何时进行维修，可最大程度地避免故障的发生。

在维修策略中，产品有2种维修方式：(1)修复性维修，即产品在达到预定维修时刻之前发生故障而进行的维修或更换活动，也称故障维修或故障更换；(2)预防性维修，即产品在达到预定的维修时刻进行的维修或更换活动。因此按照维修决策的目标，并根据剩余寿命概率密度函数，建立以最小费用为目标的维修决策模型，其决策过程如图3所示。

图3 以费用最小为目标的维修决策示意图

4.1 假设和参数

(1)每隔Δt进行一次状态检测并获取当前时刻的状态数据。

(2)如果在检测间隔期内发生故障，则立刻进行修复性维修。

(3)设备的预防性维修费用小于故障后的修复性维修费用，即Tp

(4)预防性维修和修复性维修可以使设备恢复初始状态，即修复如新。

(5)参数：E(C)为更新周期内的期望总费用；E(T)为期望更新周期长度；Δt为状态检测间隔期，Δt=ti-ti-1；f(τi|y0:i)为ti时刻状态信息为y0:i时产品的剩余寿命概率密度分布函数；Pp为部件进行预防性维修的概率；Pf为部件进行修复性维修的概率；cp为预防性维修费用；cf为修复性维修费用；TR为最佳维修更换时间；C(TR)为部件维修更换时间为TR时，长期使用下的单位时间费用；ti为当前时刻的状态监控点，t0=0；Tp为预防性维修所需的平均时间；Tc为修复性维修所需的平均时间；ci为每次实施状态检测的费用。

4.2 以费用为目标的决策模型

设备部件在长期使用下的单位时间费用可表示为

(16)

E(C)=cfPf+cpPp+nci=cff(τi

ti|y0:i)+cp[1-f(τi

cp+(cf-cp)f(τi

(17)

(18)

(19)

将(17)～(18)式代入(16)式中，则费用模型为

(20)

其中，nci为更新周期内监控所产生的费用，由于状态检测中检测费用大多为设备投入费用，nci常常可以忽略不计，因此在实际决策中，近似取nci=0。如果Tp和Tc相对于TR来说很小，在费用模型中也可以将其忽略。

利用数值方法最小化(20)式，得出设备进行基于状态的维修时可达到的最小单位费用及所对应的维修时刻TR。当获取新的状态信息值y0:i+1时，需要将数据代入(20)式重新计算更新结果。

4.3 最佳维修时间确定方法

通过费用决策模型可以得到各检测时刻对应的维修时间TR，使得单位时间维修费用最小。当TR-ti>Δt时，不对部件进行维修，直到下次正常检测时刻；当首次出现TR-ti≤Δt时，则到达TR时对部件进行维修或更换，此时的TR即为基于状态维修的最佳维修时间。

5 实例分析

试验采用3套Rexnord ZA-2155双列滚子轴承，对2套轴承进行全寿命试验，用于模型参数估计和模型验证，对第3套轴承进行截尾寿命试验，用于维修决策，分别记为轴承1、轴承2和轴承3。试验得到2组全寿命试验数据和1组截尾寿命数据。

全寿命试验过程中，轴承的转速和载荷保持在2 000 r/min和10 000 LB(约44 500 N)，并采用PCB加速度传感器对轴承的振动信号进行采集，每1 h采集1次信号，采样频率为20 kHz，长度为1 s。

试验过程中，轴承1和轴承2分别在966 h,982 h出现强烈的噪声，此时认为轴承故障。用振动信号的能量均值作为间接状态观测数据y，磨损量为状态数据x，2组全寿命数据的振动信号能量均值在整个寿命周期上的变化如图4所示。

图4 试验轴承平均能量的变化过程

应用轴承1的间接状态观测数据对模型参数进行估计，并采用EM-PF算法对剩余寿命模型中的参数进行计算，结果见表1。

表1 模型参数计算结果

确定模型参数之后就能得到基于Gamma退化过程的状态空间模型，通过计算求得轴承剩余寿命概率密度函数f(τi|y0:i)。对于轴承2的观测数据，运用Matlab 7.0对剩余寿命预测方法进行实现，不同状态检测时刻(500～980 h)的剩余寿命概率密度函数如图5所示，其中*表示剩余寿命的估计值，◁表示剩余寿命真实值。剩余寿命预测值与实际值见表2。

图5 不同检测时刻剩余寿命概率密度

表2 剩余寿命预测值与实际值比较(轴承2) h

由上述结果可以看出，该模型对于此类试验背景下的轴承剩余寿命预测具有一定的实用性，下面利用该模型对轴承3进行维修决策，确定使得单位时间费用最小的最佳维修时间。假设cp=600元，cf=1 200元，在得到轴承的剩余寿命概率密度函数之后，将f(τi|y0:i)代入到(2)式中，计算各维修时刻的单位时间费用，得到每次检测时费用最低的更换时间TR，如图6所示。图中可以看出，在某一检测时刻，随着更换时间的增加，单位时间费用先减少后增加，在极值点处可以取得该次检测的最佳维修时间。图中◁表示各检测时刻使单位时间费用最低的维修时间。

图6 不同检测时刻所确定的维修更换时间

结合最佳维修时间确定的方法，通过Matlab计算，可以得出首次满足TR-ti≤Δt条件时，TR=701 h。所以在701 h对轴承进行更换，可以使单位时间费用C(TR)达到最低,约为289.6元。

6 结束语

为解决维修决策的实际问题，重点研究了剩余寿命预测方法和维修决策优化模型。首先根据产品状态的退化过程的特征建立了基于Gamma 退化过程的状态空间模型，并对模型参数的估计方法进行了说明，进而通过该模型实现剩余寿命的预测，得到产品的剩余寿命概率密度函数。最后以滚子轴承寿命试验得到的数据为例，验证了剩余寿命预测模型，并以费用最小为目标建立了维修决策模型，对另一轴承进行了维修决策的优化，实现了剩余寿命预测与维修决策的有机结合，实例证明了所建立模型的有效性和可行性。