基于随机效应Wiener退化模型的剩余寿命预测

冯海林，李秀秀

(西安电子科技大学 数学与统计学院， 陕西 西安 710126)

0 引 言

随着科学技术的快速发展，航空航天、电子工业、机械等领域出现了大量具有高可靠性与长寿命的产品，即在常应力下很难观测到产品失效,在加速寿命试验中也很少甚至无失效数据. 这给传统的基于寿命数据的可靠性分析带来了极大挑战. 在此背景下，基于性能退化的建模思想，为这些长寿命、高可靠性产品的可靠性分析提供了新的途径[1-2]. 从产品性能退化的角度分析其失效情况，一般需要选择一个与产品寿命或可靠性高度相关且可以测量的或从测量数据中可以提取的变量，称其为性能可靠性特征量[3-4].刻画产品性能退化所提取的特征量称为性能退化量，性能退化量受产品使用环境等因素的影响，且随时间的变化而变化，具有一定的随机性.因此，建立合理的性能退化随机过程模型是准确预测产品剩余寿命的关键环节之一. 基于标准Wiener过程所衍生的性能退化过程模型是近年来广受关注且可描述多种典型产品性能退化的一类随机过程，其一般形式为

X(t)=vΛ(t)+σB(Λ(t)),

(1)

其中，{X(t)：t≥0}表示产品性能的退化量.v是漂移参数，σ是扩散参数，且v和σ均为常数.B(·)是标准的Wiener过程，Λ(·)通常表示时间尺度的单调增函数，且Λ(0)=0和B(0)=0.

根据以上随机效应Wiener过程模型和当前的性能退化量，建立剩余寿命的退化模型，获得剩余寿命的概率密度函数和可靠度函数. 通过分析退化模型参数，可知模型的参数含有随机效应，对数似然函数比较复杂，很难最优化. 因此，选用EM算法进行参数估计，再融合同类产品的退化数据，对该Wiener退化模型中的未知参数进行估计，并与已有模型进行对比，以说明本模型的实用性和精确性.

1 含随机效应的Wiener退化模型

令X(t)为产品在t时刻的性能退化量的测量值，建立随机效应Wiener退化模型为

X(t)=vΛ(t)+(ζv+a)B(Λ(t)),

(2)

由于产品制造工艺、制造材料的差异以及受工作环境的动态影响，同类型不同产品的性能退化过程存在差异. 为描述这种差异性，在退化模型中常采用将相关参数随机化的方法，即将Wiener退化模型(1)中的漂移参数设成正态分布.

模型(2)与模型(1)相比，有σ=ζv+a，表示v和σ呈线性关系. 通过加入误差项a，缩小了性能退化测量值与真实值之间的偏差，从而提高了产品可靠性评估的精确度. 该Wiener模型不仅体现了产品个体在同时具有高退化率与高波动性的性质，而且也解决了产品的退化性能在测量时会产生误差的问题；同时，大大降低了计算难度，更符合实际产品的退化数据，从而能更准确地预测产品的剩余寿命.

根据Wiener退化模型(1)，考虑漂移参数的随机效应，可知模型(2)中退化量X(t)的概率密度函数(PDF)为

(3)

2 基于Wiener退化过程的剩余寿命预测

本文的主要目的是利用新随机效应Wiener退化模型(2)描述产品的退化轨迹，最终预测产品的剩余寿命. 下面将由首达时的概念推导随机效应Wiener退化模型(2)计算的剩余寿命的分布形式. 为此，首先给出Wiener过程首达时的分布.

2.1 Wiener过程首达时的分布

给定产品的失效阈值ξ,产品的寿命T为性能退化量X(t)首次达到失效阈值ξ的时间，则

T=inf(t|X(t)≥ξ,t>0).

(4)

考虑模型(2)中的漂移参数μ和a是随机变量，则产品的寿命T的根率密度函数(PDF)为

(5)

其相应的概率分布函数为

(6)

因此，产品的可靠度函数可表示为

(7)

2.2 产品的剩余寿命预测

为预测某产品的剩余寿命，首先须建立此产品的剩余寿命预测模型. 假定特定产品在时刻tk和t+tk的性能退化量为X(tk)和X(t+tk). 按照Wiener过程独立增量的性质，可得

X(t+tk)=X(tk)+vΛ(t)+(ζv+a)B(Λ(t)).

(8)

产品的剩余寿命是指从当前时刻到产品失效时刻的时间间隔. 如果已知当前的性能退化量为X(tk)，根据寿命的定义可知产品在时刻tk的剩余寿命t为

M=inf(t|X(t+tk)≥ξ,t>0).

(9)

获得产品剩余寿命的关键是找到剩余寿命的概率密度函数(PDF).由Wiener过程独立增量性质知,

M= inf(t|X(t+tk)≥ξ,t>0)=

inf(t|X(t+tk)-X(tk)≥ξ-X(tk),t>0)=

inf(t|X(t)≥ξ-X(tk),t>0),

(10)

则产品的剩余寿命t的PDF为

(11)

3 Wiener退化模型的参数估计

由退化数据估计参数

(12)

其中β为Λ(t;β)中未知参数的集合. 当一个产品的性能退化量可观察时，能够更新v的分布. 假定可获得N个实验产品的退化数据，当时间为(t1,t2,…,tn)′时观察的退化量为x=(x1,x2,…,xn)′，令Δx=(Δx1,Δx2,…,Δxn)′，其中Δx1=x1和Δxj=xj-xj-1. 类似地，令λ1=Λ(t1)和λj=Λ(tj)-Λ(tj-1),1

(13)

其中，

令vi为第i个产品已实现的随机效应. 以vi为条件，第i个产品的退化量{Xi(t);vi,t>0}有独立的增量. ΔXi=(ΔXi1,ΔXi2,…,ΔXini)为退化量增量，其中ΔXij=X(tij)-X(tij-1),j=1,2,…,ni，ΔXi1=X(ti1). Δxi是增量的观察值，N个产品的观察数集D={Δxi,i=1,2,…,N}. 由对数似然函数得到极大似然估计，式(13)可改写为式(14).

(14)

因式中含有随机效应，对数似然函数较复杂，很难最优化. 因此，可用EM算法估计新随机效应Wiener退化模型中的未知参数. 给定参数估计θ(m)，经EM算法反复迭代Eθ(m)[l(θ)]后，用θ更新θ(m).

EM算法的求解分2步：

E-step给定观测数据集D={Δxi,i=1,2,…,N}和初始参数值θ(0)，计算对数似然函数l(θ)关于未知参数θ的期望. 定义对数似然函数的期望Q(θ|θ(m))=E[l(θ;D)|θ(m),D],θ(m)表示m次迭代的估计参数 .vi关于Δxi和θ的条件分布为

(15)

由式 (15)可改写为式(14)，有

(16)

第2部分为

(17)

利用Q2(θ|θ(m))对ζ2求导，且令导数为0，式(18)改写为式(17)：

(18)

将式(18)代入式(17)，得到

(19)

从β开始最大化式(18)，得到β(m+1)，且将结果代入式(17)，得到(ζ2)(m+1).

4 数据分析

基于本文的Wiener退化模型数据分析，以文献[8]中2024-T351铝合金样本的疲劳裂纹扩展数据为例(在同一恒定振幅下获得数据)，所有样品的初始裂纹长度相同；每个样品的裂纹长度增量是在10 000个时间周期后测量的，每间隔5 000个时间周期测量1次，直到40 000个时间周期止. 实验数据见表1[8]，得到的裂纹增长路径图见图1.

表1 5种Wiener过程模型的MLE和AIC

图1 30个样本的裂纹增长路径Fig.1 Fatigue crack growth paths of 30 testing units

WU等[8]通过纵向数据分析技术拟合数据，本文使用Wiener过程分析此组数据. 研究中，简单模型式(1)记为SM，其中参数v和σ为固定常数；由YE等[7]提出的随机效应模型记为RDV，其漂移参数v～N(μ，w2)，扩散参数σ=ζv；由PENG 等[5]提出的随机效应模型记为RD，其v～N(μ，w2)；由WANG[6]提出的随机效应模型记为RDRV，其漂移参数[v|σ]～N(μ,θσ2)，扩散参数σ-2～Ga(r,δ)；本文的随机效应模型记为NRDV,其漂移参数v～N(μ，w2)，扩散参数σ=ζv+a；对这几种模型进行参数估计和AIC值计算.

通常认为裂纹增长遵循幂律，表示为Λ(t)=tb. 表1为SM(固定效应)、RDV(随机效应模型)、RD(单一效应模型)、RDRV(混合效应模型)和NRDV(新随机效应模型)5种模型参数的极大似然估计及对应的AIC(Akaike information criterion).

AIC=2k-2 lnl，

其中，k为参数的数量，l为模型似然函数的最大值. AIC值越小，表示模型越符合实际. 由表1 知，“NRDV”最符合实际数据，因此，由该模型预测的剩余寿命更精确.

图2 参数和之间的散点图Fig.2 Scatter plot of the estimated drift parameter and estimated volatility parameter .

从图2中可以看出，估计的漂移参数和扩散参数呈线性分布，且一次拟合系数为 0.583 6. 说明当产品有更高的退化率时，其退化路径会随退化率的变化而波动. 本文提出的模型可以更好地解释这一现象.

图3 参数vi的标准正态Q-Q图Fig.3 Standard normal Q-Q plot for parameter

5 结 论

高退化率产品也具有高波动性，这是一种常见现象. 由于均值和方差是2个互相独立的参数，以前研究的Wiener退化过程不易描绘上述问题. 为了克服这一缺陷，本文提出了新随机效应的Wiener过程模型，即在简单的Wiener退化模型上，引入随机效应，对v施加一个正态分布，并将扩散参数与漂移参数设定成线性关系， 此模型较现有随机效应维纳过程模型(扩散参数和漂移参数都相互独立)更符合实际情况. 对新的随机效应模型进行统计推断，并用EM算法估计模型参数，验证了本文模型的有效性和准确性.

本模型可用于退化试验，以评估产品的可靠性，也可预测已知当前退化量产品的剩余寿命. 结果表明，本模型更接近实际应用. 因此，利用本模型估计的可靠性和预测的剩余使用寿命更准确. 同时，本模型还具有扩展性，而现有随机效应模型中的假设是有所限制的. 需要指出的是： 除了正态分布外，截断正态分布也常用来表示随机效应分布.