文章快速检索     高级检索
  浙江大学学报(理学版)  2018, Vol. 45 Issue (6): 679-684, 693  DOI:10.3785/j.issn.1008-9497.2018.06.006
0

引用本文 [复制中英文]

冯海林, 李秀秀. 基于随机效应Wiener退化模型的剩余寿命预测[J]. 浙江大学学报(理学版), 2018, 45(6): 679-684, 693. DOI: 10.3785/j.issn.1008-9497.2018.06.006.
[复制中文]
FENG Hailin, LI Xiuxiu. Residual life prediction based on a stochastic effect Wiener degradation model[J]. Journal of Zhejiang University(Science Edition), 2018, 45(6): 679-684, 693. DOI: 10.3785/j.issn.1008-9497.2018.06.006.
[复制英文]

基金项目

国家自然科学基金资助项目(201300001)

作者简介

冯海林(1966-), ORCID:http://orcid.org/0000-0002-2280-3248, 女, 博士, 教授, 主要从事系统可靠性预测等研究

通信作者

李秀秀, ORCID:http://orcid.org/0000-0001-6492-151X, E-mail:xxli-1@stu.xidian.edu.cn

文章历史

收稿日期:2017-12-04
基于随机效应Wiener退化模型的剩余寿命预测
冯海林 , 李秀秀     
西安电子科技大学 数学与统计学院, 陕西 西安 710126
摘要: 针对退化率较高的产品具有不稳定的退化路径以及产品个体差异对退化过程的影响,建立了一种新的随机效应退化模型,即漂移参数和扩散参数均为随机变量且两者之间呈线性关系的Wiener退化过程模型.基于该模型获得了产品剩余寿命分布与可靠度函数,同时设计了估计模型参数的EM(expectation maximization)算法.最后,通过分析钛合金疲劳裂纹数据以及与现有模型结果的比较,验证了所建模型的有效性和准确性.
关键词: Wiener过程    随机效应    剩余寿命    EM算法    
Residual life prediction based on a stochastic effect Wiener degradation model
FENG Hailin, LI Xiuxiu     
School of Mathematics and Statistics, Xidian University, Xi'an 710126, China
Abstract: Since the products with high degradation rate don't conform to a stable degradation path, a new random effect degradation model namely Wiener degradation process model is established accounting for the influence of individual differences of products on the degradation process. Both drift parameter and diffusion parameter of the model are random variables, and the relationship between them is linear. Then, the residual life distribution and reliability function of the products are obtained based on the model, simultaneously, the EM algorithm to estimate parameters of the model is designed. Finally, the validity and accuracy of the model are verified by analyzing the fatigue crack data of a kind of titanium alloy and compared with the existing model.
Key Words: Wiener process    random effects    residual life    the EM(expectation maximization) algorithm    
0 引言

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

$ X\left( t \right) = v\mathit{\Lambda }\left( t \right) + \sigma B\left( {\mathit{\Lambda }\left( t \right)} \right), $ (1)

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

实际应用中,产品性能的退化具有个体差异性,为此,PENG等[5]首次将个体差异性视为随机影响引入Wiener退化模型的漂移参数中. WANG[6]将个体差异性建模为漂移参数和扩散参数随机化,并用EM算法估计模型参数. YE等[7]提出了一类新的随机效应Wiener退化模型,漂移参数是随机的,扩散参数与其正相关,符合产品高退化率和较高波动性的高可靠性特点.实际上,受测量工具及使用方法等因素影响,测得的产品性能退化(疲劳裂纹长度、磁头磨损等)量存在一定误差,文献[7]未考虑测量误差对退化率与波动性的影响.本文在文献[7]模型的基础上,建立了新的随机效应Wiener过程模型,通过引入误差变量aa~N(μa, wa2),即考虑了测量误差的影响,使模型更符合实际产品的退化数据.但模型精确化的同时也令计算更困难.本文的另一工作是解决参数估计问题.

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

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

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

$ X\left( t \right) = v\mathit{\Lambda }\left( t \right) + \left( {\zeta v + a} \right)B\left( {\mathit{\Lambda }\left( t \right)} \right), $ (2)

其中,{X(t):t≥0}表示产品性能的退化量, v为漂移参数,ζv+a为扩散参数,a为误差项,表示测量误差,且va都是随机的,即v~N(μ, w2), a~N(μa, wa2),B(·)是标准的Wiener过程,Λ(·)通常表示时间尺度的单调增函数,且Λ(0)=0和B(0)=0.

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

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

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

$ \begin{array}{l} {f_{X\left( t \right)}}\left( t \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right)\mathit{\Lambda }\left( t \right)} }} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left( { - \frac{{{{\left( {x - \mu \mathit{\Lambda }\left( t \right)} \right)}^2}}}{{2\left[ {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right]\mathit{\Lambda }\left( t \right)}}} \right). \end{array} $ (3)
2 基于Wiener退化过程的剩余寿命预测

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

2.1 Wiener过程首达时的分布

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

$ T = \inf \left( {t\left| {X\left( t \right) \ge \xi ,t > 0} \right.} \right). $ (4)

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

$ \begin{array}{l} {f_T}\left( t \right) = \frac{\xi }{{\sqrt {2{\rm{ \mathsf{ π} }}\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right)\mathit{\Lambda }\left( t \right)} }} \times \\ \;\;\;\;\;\;\;\;\;\;\;\exp \left( { - \frac{{{{\left( {\xi - \mu \mathit{\Lambda }\left( t \right)} \right)}^2}}}{{2\left[ {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right]\mathit{\Lambda }\left( t \right)}}} \right), \end{array} $ (5)

其相应的概率分布函数为

$ \begin{array}{l} F\left( t \right) = \mathit{\Phi }\left( { - \frac{{\mu \mathit{\Lambda }\left( t \right) - \xi }}{{\sqrt {\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right){\mathit{\Lambda }^3}\left( t \right)} }}} \right) + \\ \;\;\exp \left( {\frac{{2\mu \left( {{\zeta ^2}{w^2} + w_a^2} \right)\xi + 2{w^2}{\xi ^2}}}{{{\zeta ^4}{w^4} + w_a^4}}} \right) \times \\ \;\;\mathit{\Phi }\left( {\frac{{2{w^2}\xi t + \left( {{\zeta ^2}{w^2} + w_a^2} \right)\left( {\mu t + \xi } \right)}}{{\left( {{\zeta ^2}{w^2} + w_a^2} \right)\sqrt {\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right)\mathit{\Lambda }\left( t \right)} }}} \right). \end{array} $ (6)

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

$ \begin{array}{l} R\left( t \right) = P\left( {T > t} \right) = 1 - \int_0^t {{f_T}\left( t \right){\rm{d}}t} = \\ \;\;\;\;\;\;1 - \mathit{\Phi }\left( { - \frac{{\mu \mathit{\Lambda }\left( t \right) - \xi }}{{\sqrt {\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right){\mathit{\Lambda }^3}\left( t \right)} }}} \right) - \\ \;\;\;\;\;\;\exp \left( {\frac{{2\mu \left( {{\zeta ^2}{w^2} + w_a^2} \right)\xi + 2{w^2}{\xi ^2}}}{{{\zeta ^4}{w^4} + w_a^4}}} \right) \times \\ \;\;\;\;\;\;\mathit{\Phi }\left( {\frac{{2{w^2}\xi t + \left( {{\zeta ^2}{w^2} + w_a^2} \right)\left( {\mu t + \xi } \right)}}{{\left( {{\zeta ^2}{w^2} + w_a^2} \right)\sqrt {\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right)\mathit{\Lambda }\left( t \right)} }}} \right). \end{array} $ (7)
2.2 产品的剩余寿命预测

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

$ X\left( {t + {t_k}} \right) = X\left( {{t_k}} \right) + v\mathit{\Lambda }\left( t \right) + \left( {\zeta v + a} \right)B\left( {\mathit{\Lambda }\left( t \right)} \right). $ (8)

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

$ M = \inf \left( {t\left| {X\left( {t + {t_k}} \right) \ge \xi ,t > 0} \right.} \right). $ (9)

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

$ \begin{array}{l} M = \inf \left( {t\left| {X\left( {t + {t_k}} \right) \ge \xi ,t > 0} \right.} \right) = \\ \;\;\;\;\;\;\;\inf \left( {t\left| {X\left( {t + {t_k}} \right) - X\left( {{t_k}} \right) \ge \xi - X\left( {{t_k}} \right),t > 0} \right.} \right) = \\ \;\;\;\;\;\;\;\inf \left( {t\left| {X\left( t \right) \ge \xi - X\left( {{t_k}} \right),t > 0} \right.} \right), \end{array} $ (10)

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

$ \begin{array}{l} {f_M}\left( t \right) = \frac{{\xi - X\left( {{t_k}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right){\mathit{\Lambda }^3}\left( t \right)} }} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\exp \left( { - \frac{{{{\left( {\xi - X\left( {{t_k}} \right) - \mu \mathit{\Lambda }\left( t \right)} \right)}^2}}}{{2\left[ {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right]\mathit{\Lambda }\left( t \right)}}} \right). \end{array} $ (11)
3 Wiener退化模型的参数估计

由退化数据估计参数

$ \theta = {\left( {\mu ,{w^2},{\mu _a},w_a^2,\zeta ,\beta '} \right)^\prime }, $ (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 < jn.基于贝叶斯理论,v分布可由式(13)更新为式(12):

$ \begin{array}{l} f\left( {v\left| {\Delta x} \right.} \right) \propto {v^n}\exp \left( { - \frac{1}{{2\xi }}\sum\limits_{j = 1}^n {\frac{{\left( {v\Delta {x_j} - {\lambda _j}} \right)}}{{{\lambda _j}}}} - \frac{{\left( {v - \mu } \right)}}{{2{w^2}}}} \right) \propto \\ \;\;\;\;\;\;\;\;\;\;\;\;{v^n}\exp \left( { - \frac{{{{\left( {v - {{\tilde \mu }_n}} \right)}^2}}}{{2\tilde w_n^2}}} \right), \end{array} $ (13)

其中,

$ {{\tilde w}_n} = {\left( {\frac{1}{{{w^2}}} + \sum\limits_{j = 1}^n {\left( {\frac{{{{\left( {\Delta {x_j}} \right)}^2}}}{{{\zeta ^2}{\lambda _j}}} + \frac{{{{\left( {\Delta {x_j}} \right)}^2}}}{{w_a^2{\lambda _j}}}} \right)} } \right)^{ - 1/2}}, $
$ {{\tilde \mu }_n} = \tilde w_n^2\left( {\frac{\mu }{{{w^2}}} + \left( {\frac{1}{{{\zeta ^2}}} + \frac{1}{{w_a^2}}} \right)\sum\limits_{j = 1}^n {\Delta {x_j}} } \right). $

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).

$ \begin{array}{l} l\left( \theta \right) = \sum\limits_{i = 1}^N {\left[ { - \frac{1}{2}\sum\limits_{j = 1}^{{n_i}} {\ln \left( {2{\rm{ \mathsf{ π} }}\left( {{w^2}\mathit{\Lambda }\left( t \right) + {\zeta ^2}{w^2} + w_a^2} \right){\lambda _{ij}}} \right)} } \right. + } \\ \;\;\;\;\;\;\;\;\;\left. {\frac{{\hat \mu _{i:{n_i}}^2}}{{2\hat w_{i:{n_i}}^2}} - \sum\limits_{j = 1}^{{n_i}} {\frac{{{\lambda _{ij}}}}{{2{\zeta ^2}{w^2}}}} - \sum\limits_{j = 1}^{{n_i}} {\frac{{{\lambda _{ij}}}}{{2w_a^2}}} - \frac{{{\mu ^2}}}{{2{w^2}}} - \frac{{{\mu ^2}}}{{2w_a^2}}} \right]. \end{array} $ (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θ的条件分布为

$ \begin{array}{l} f\left( {{v_i}\left| {\Delta {x_i}} \right.} \right) = {\left[ {\sqrt {2{\rm{ \mathsf{ π} }}} {w_{\bar i:{n_i}}}\sum\limits_{j = 0}^{{n_i}} {\left( {C_{{n_i}}^j\bar \mu _{i;{n_i}}^{{n_i}}\bar w_{i;{n_i}}^j} \right)} } \right]^{ - 1}} \times \\ \;\;\;\;\;v_{\rm{i}}^{{n_i}}\exp \left( {\frac{{{{\left( {{v_i} - {{\bar \mu }_{i;{n_i}}}} \right)}^2}}}{{2\bar w_{i;{n_i}}^2}}} \right). \end{array} $ (15)

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

$ E\left[ {{v_i}\left| {\Delta {x_i}} \right.} \right] = \frac{{\sum\limits_{j = 0}^{{n_i} + 1} {\left( {C_{{n_i}}^j\bar \mu _{i;{n_i} + 1}^{{n_i} + 1}\bar w_{i;{n_i} + 1}^j} \right)} }}{{\sum\limits_{j = 0}^{{n_i}} {\left( {C_{{n_i}}^j\bar \mu _{i;{n_i}}^{{n_i}}\bar w_{i;{n_i}}^j} \right)} }}, $
$ E\left[ {v_i^2\left| {\Delta {x_i},{\theta ^{\left( m \right)}}} \right.} \right] = \frac{{\sum\limits_{j = 0}^{{n_i} + 2} {\left( {C_{{n_i}}^j\bar \mu _{i;{n_i} + 2}^{{n_i} + 2}\bar w_{i;{n_i} + 1}^j} \right)} }}{{\sum\limits_{j = 0}^{{n_i}} {\left( {C_{{n_i}}^j\bar \mu _{i;{n_i}}^{{n_i}}\bar w_{i;{n_i}}^j} \right)} }}. $

然后用uivi分别表示E[vixi, θ(m)]和E[vi2xi, θ(m)]. θ(m)中相应的参数通过依次迭代参数$ \widetilde \mu $i; ni$ \widetilde w $i; ni获得ui(m)vi(m).

M-step  最大化期望值Q(θ|θ(m)),即找到一个θ(m+1)θ(m+1)=$ {\rm{argma}}\mathop {\rm{x}}\limits_\theta Q(\theta |\theta ^{(m)}) $,由式(14)知,Q函数可分解成两部分.第1部分为

$ \begin{array}{l} {Q_1}\left( {\theta \left| {{\theta ^{\left( m \right)}}} \right.} \right) = - \frac{1}{2}\sum\limits_{i = 1}^N {\left[ {\ln \left( {{w^2} + w_a^2} \right) + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{{v_i^{\left( m \right)} - 2\left( {\mu + {\mu _a}} \right)u_i^{\left( m \right)} + {{\left( {\mu + {\mu _a}} \right)}^2}}}{{{w^2} + w_a^2}}} \right], \end{array} $ (16)

式(16)可改写为式(15),可求得μawa2μw2的最大值:

$ {\mu ^{\left( {m + 1} \right)}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {u_i^{\left( m \right)} - {{\left( {{\mu _a}} \right)}^{\left( {m + 1} \right)}}} \right)} , $
$ {\left( {{\mu _a}} \right)^{\left( {m + 1} \right)}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {u_i^{\left( m \right)} - {{\left( \mu \right)}^{\left( {m + 1} \right)}}} \right)} , $
$ \begin{array}{l} {\left( {{w^2}} \right)^{\left( {m + 1} \right)}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {v_i^{\left( m \right)} - {{\left( {w_a^2} \right)}^{\left( {m + 1} \right)}}} \right)} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {\frac{1}{N}\sum\limits_{i = 1}^N {\left( {u_i^{\left( m \right)} - {{\left( {w_a^2} \right)}^{\left( {m + 1} \right)}}} \right)} } \right)^2}, \end{array} $
$ \begin{array}{l} {\left( {w_a^2} \right)^{\left( {m + 1} \right)}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {v_i^{\left( m \right)} - {{\left( {{w^2}} \right)}^{\left( {m + 1} \right)}}} \right)} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {\frac{1}{N}\sum\limits_{i = 1}^N {\left( {u_i^{\left( m \right)} - {{\left( {{w^2}} \right)}^{\left( {m + 1} \right)}}} \right)} } \right)^2}; \end{array} $

第2部分为

$ \begin{array}{l} {Q_2}\left( {\theta \left| {{\theta ^{\left( m \right)}}} \right.} \right) = - \frac{1}{2}\sum\limits_{i = 1}^N {\left[ {{n_i}\ln {\zeta ^2} + \sum\limits_{j = 1}^{{n_i}} {\ln {\lambda _{ij}}} + } \right.} \\ \;\;\;\;\;\;\;\left. {\frac{1}{{{\zeta ^2}}}\sum\limits_{j = 1}^{{n_i}} {\frac{{v_i^{\left( m \right)}{{\left( {\Delta {x_{ij}}} \right)}^2} - 2{\lambda _{ij}}u_i^{\left( m \right)}\Delta {x_{ij}} + \lambda _{ij}^2}}{{{\lambda _{ij}}}}} } \right], \end{array} $ (17)

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

$ {\zeta ^2} = {\left( {\sum\limits_{i = 1}^N {{n_i}} } \right)^{ - 1}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^{{n_i}} {\frac{{v_i^{\left( m \right)}{{\left( {\Delta {x_{ij}}} \right)}^2} - 2{\lambda _{ij}}u_i^{\left( m \right)}\Delta {x_{ij}} + \lambda _{ij}^2}}{{{\lambda _{ij}}}}} } . $ (18)

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

$ \begin{array}{l} Q_2^p\left( {\beta \left| {{\theta ^{\left( m \right)}}} \right.} \right) = \\ - \frac{1}{2}\left[ {\ln \left( {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^{{n_i}} {\frac{{v_i^{\left( m \right)}{{\left( {\Delta {x_{ij}}} \right)}^2} - 2{\lambda _{ij}}u_i^{\left( m \right)}\Delta {x_{ij}} + \lambda _{ij}^2}}{{{\lambda _{ij}}}}} } } \right) + } \right.\\ \left. {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^{{n_i}} {\ln {\lambda _{ij}}} } } \right], \end{array} $ (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 Table 1 MLEs and AICs of the five Wiener process models
图 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).

$ {\rm{AIC}} = 2k - 2\ln l, $

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

下面给出NRDV最优的解释.使用简单的Wiener退化模型(1)分别验证30个样本的退化路径,测得的30组数据均有一个共同的指数参数b.分别给出该数据的30对估计值($ {\hat v_i}, {\hat \sigma _i} $),如图 2所示.

图 2 参数$ {\hat v_i} $$ {\hat \sigma _i} $之间的散点图 Fig. 2 Scatter plot of the estimated drift parameter $ {\hat v_i} $ and estimated volatility parameter $ {\hat \sigma _i} $.

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

为进一步评估该模型的拟合优度,注意到(Δxij-viλij)/(ζviλij1/2)是独立的标准正态随机变量,vi是第i个单元的随机效应,因此可以考虑分位数(Q-Q)图.然而, 在标准的正态随机变量中所涉及的所有参数都应由数据的估计值代替,则拟合的差异取决于参数值和样本大小.因此,本文是在相同情景和相同样本下模拟获得的极大似然估计值,并以该估计值为参数进行Q-Q图评估,见图 3.由图 3知,该模型拟合良好.

图 3 参数vi的标准正态Q-Q图 Fig. 3 Standard normal Q-Q plot for parameter $ {\hat v_i}$
5 结论

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

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

参考文献
[1] MEEKER W Q, HAMADA M. Statistical tools for the rapid development and evaluation of high-reliability products[J]. IEEE Transactions on Reliability, 1995, 44(2): 187–198. DOI:10.1109/24.387370
[2] LU C J, MEEKER W Q, ESCOBAR L A. A comparison of degradation and failure-time analysis methods for estimating a time-to-failure distribution[J]. Statistica Sinica, 1996, 6(3): 531–546.
[3] 金光. 基于退化的可靠性技术[M]. 北京: 国防工业出版社, 2014.
JIN G. Reliability Technology Based on Degradation[M]. Beijing: National Defense Industry Press, 2014.
[4] 訾佼佼, 刘宏昭, 蒋喜, 等. 基于退化量分布的电主轴可靠性评估[J]. 中国机械工程, 2014, 25(6): 807–812.
ZI J J, LI H Z, JIANG X, et al. Reliability assessment of electric spindle based on degradation values distribution[J]. China Mechanical Engineering, 2014, 25(6): 807–812. DOI:10.3969/j.issn.1004-132X.2014.06.020
[5] PENG C Y, TSENG S T. Mis-specification analysis of linear degradation models[J]. IEEE Transactions on Reliability, 2009, 58(3): 444–455. DOI:10.1109/TR.2009.2026784
[6] WANG X. Wiener processes with random effects for degradation data[J]. Journal of Multivariate Analysis, 2010, 101(2): 340–351. DOI:10.1016/j.jmva.2008.12.007
[7] YE Z S, CHEN N, SHEN Y. A new class of Wiener process models for degradation analysis[J]. Reliability Engineering & System Safety, 2015, 139: 58–67.
[8] WU W F, NI C C. A study of stochastic fatigue crack growth modeling through experimental data[J]. Probabilistic Engineering Mechanics, 2003, 18(2): 107–118. DOI:10.1016/S0266-8920(02)00053-X