竞争失效场合恒定应力加速退化试验统计分析

曹宇宇，唐家银，何平，曹玲

(西南交通大学数学学院, 四川 成都 611756)

0 引言

随着科技化水平的不断进步，很多领域都出现了高可靠性、长寿命的产品，如高性能的电子产品、航天产品，因此用传统的可靠性试验会花费大量的时间和成本，且有些传统的可靠性试验方法在工程中较难实现.为了缩短试验周期、降低试验成本，新的可靠性评估方法应运而生.其中，加速退化试验(accelerated degradation testing, ADT)和加速寿命试验(accelerated life testing, ALT)目前已成为可靠性试验领域的两大重要组成部分[1]，受到国内外研究学者的关注.

ADT是指通过提高应力水平来加速产品退化，通过高应力下的数据来估计产品的可靠性并预测常应力下产品的寿命[2].传统的加速退化试验统计分析中通常假设产品仅有单一失效模式，但在实际工程中，产品存在多种失效模式的可能性很大，任何一个失效模式发生，产品则发生失效.从失效模式来看，产品的失效可分为退化失效与突发失效.在实际工作中，产品发生失效是由最早出现的失效模式导致的.因此在竞争失效场合下，传统单一失效模式的统计分析方法并不适用.

关于竞争失效的研究，最早是由Nelson[3]建立竞争失效模型，通过分析给出了对数正态分布下的极大似然估计(MLE).赵建印[4]通过对金属化膜脉冲电容器的研究，得到了竞争失效下的产品可靠性模型.王华伟等[5]利用贝叶斯模型建立航空发动机性能退化模型,实现了在不同退化情况下航空发动机剩余寿命的预测.安宗文等[6]通过Wiener过程对连续退化进行刻画,推导出产品在失效相关模式下的竞争失效可靠性分析模型.龙哲等[7]通过对退化失效采用Wiener过程来描述产品的性能退化轨迹.上述基于伪失效寿命研究中使用的退化模型大多假设为服从简单的一般线性过程或者Wiener过程，并不能很好地拟合退化轨迹，而且基于退化量分布更加关注的是失效过程本身，比起伪寿命结果可信度更高.苏春等[8]则基于性能退化数据对竞争失效模型进行了研究和分析，但研究过于简单，不能够全面地反映竞争失效问题.另外，对竞争失效下恒定应力加速退化试验也鲜少有较为详尽的研究.

为了对突发失效模型进行更好地刻画，引入基于比例危险的突发失效模型，比例危险模型最早由COX[9]提出并应用于寿命数据分析，该模型主要采用非参数统计方法获取产品的可靠度估计，EGHBALI[10]改进了COX的比例危险模型，得到了一种融合加速退化因子的比例危险退化模型；蔡忠义等[11]将监测时刻作为协变量，建立了比例危险退化模型.王智明等[12]采用混合Weibull分布并结合比例危险模型对刀具故障数据进行建模.以上方法均是基于退化数据，而本研究则对比例危险模型进行了改进，使其应用于突发失效数据，并由此获得基于比例危险的突发失效模型.

笔者针对竞争失效最一般的形式(产品具有单一退化失效模式和单一突发失效模式)，对其恒定应力加速退化试验进行建模和分析.首先给出竞争失效场合恒定应力加速试验的突发与退化竞争模型，在此基础上分别给出突发失效模型和加速退化模型，通过对参数进行估计，实现对竞争失效下恒定应力加速退化试验的统计分析.最后给出应用算例.

1 竞争失效场合恒定应力的可靠性评估模型

1.1 模型假设与符号竞争失效场合恒加试验统计分析基于以下假定：

假设1：产品具有单一突发失效模式和单一退化失效模式，其相应的突发失效时间为Tr，退化失效时间为Td，两者相互竞争；

假设2：不同时刻、不同应力下，产品性能退化量X(t)服从相同的分布族，参数为时间t和应力水平S的函数；

假设3：失效阈值D是恒定的；

假设4：性能退化程度可能会影响突发失效，也可能不会.

1.2 突发与退化竞争模型设在应力水平Si(i=1,2,…,k)下，时刻t产品的退化量为X(t)，且当X(t)首次达到失效阈值D时产品发生退化失效.进一步假设产品突发失效后功能立刻完全丧失，因此不能继续测量其退化量.记退化失效时间为Td,突发失效时间为Tr，则产品失效时间T为

T=min{Td,Tr}

(1)

时间t内产品不失效的概率，即可靠度为

R(t)=P{T>t}=P{Td>t,Tr>t}

(2)

在退化失效模式下，产品的可靠度函数为

(3)

式中：Rd(t,S)是仅考虑应力S下产品退化失效模式时的可靠度，gd(t,xt,S)为应力S下退化量在t时刻的分布密度函数，xt表示t时刻X(t)的取值，即X(t)=xt.

为了得到恒定应力下系统竞争失效的可靠度，需要考虑性能退化和突发失效的相关性.若突发与退化相关，则产品突发失效的失效率λr(t,xt,S)是时间t、退化量xt和应力水平S的函数，Tr为突发失效的时间，则突发失效模式下的可靠度函数为

(4)

式中：Rr(t,xt,S)表示应力S下突发失效的可靠度函数，xt表示突发失效发生时刻t对应的退化量值，即X(t)=xt，τ为时间t的积分变量.

因此，竞争失效模式下t时刻产品的可靠度为

(5)

如果突发失效与退化量相互独立，此时突发失效率可以表示为

λr(t,xt,S)=λr0(t,S)

(6)

式中：λr0(t)是不考虑性能退化影响情况下突发失效的失效率.

此时产品的可靠度为

(7)

以下研究仅考虑突发失效与退化量相关的情况.当突发与退化不相关时，这时可采用传统的构造分布函数来刻画突发失效的失效模型，再利用极大似然即可得到突发失效模型中的分布参数.

根据竞争失效模型，在获得应力水平Si下竞争失效的可靠度函数后，可令应力水平为S0，即可得到在常应力S0下产品的各项可靠度估计.

1.3 加速退化方程加速退化方程[13]定义如下：记ζ(t)为产品退化量x(t)的某特征值或退化过程中的某参数，方程ζ(t)=A(t)eb(t)φ(S)或lnζ(t)=a(t)+b(t)φ(S)被称为加速退化方程.其中S为应力水平，φ(S)是应力的函数，a(t)=ln(A(t))和b(t)是与退化量和应力无关的参数.

在实际工程应用中，当温度作为加速应力时，φ(S)=1/S.当电应力作为加速应力时，φ(S)=lnS.

1.4 基于比例危险的突发失效模型为了对恒定应力加速退化当中突发失效进行更好的描述，本研究给出基于比例危险的突发失效模型.在恒定应力加速退化试验当中，比例危险(proportional hazards，PH)族[11]是具有如下性质的模型：同一时刻不同退化程度的两个个体的失效率函数(失效率函数也被称为危险函数)与时间无关，仅与退化程度和当时应力有关.设应力Si下时刻t两个产品的退化量分别为x1和x2，则这两个产品在时刻t的失效率之比λ(t,x1,S)/λ(t,x2,S)与时间t无关.对于突发失效时间T来说，如果其失效率函数具有如上性质，则其条件失效率函数可写成

λ(t,xt,S)=λ0(t)h(xt,S)

(8)

式中：λ0(t)是不考虑性能退化和应力水平影响下突发失效的失效率，也叫做基准失效率函数；h(xt,S)是性能退化水平和应力水平的函数.

在该试验下，可令h(xt,S)=exp(α1+α2xt+α3S)，式中：α1、α2、α3为未知回归系数.

对于上述基于比例危险的突发失效模型可以看到，获得基准失效率函数即λ0(t)是该模型的关键所在，估计λ0(t)之后可以进一步得到突发失效的失效率函数，从而对竞争失效问题得到更清晰的认识.

2 竞争失效场合恒定应力加速退化统计分析方法

2.2 突发与退化模型的参数估计

2.2.1 退化失效模型中的参数估计 在只考虑退化失效模式下，根据性能退化数据的散点图特征，利用拟合优度检验选取最合适的性能退化分布函数.常用的退化量分布函数有对数正态分布、威布尔分布、正态分布等分布类型.

当退化量分布为位置-刻度族模型时，一般地，位置参数可表示为退化时间(或其变换形式)、加速应力(或其变换形式)的线性函数：

μ(t,S;β)=β1+β2t+β3φ(S)

(9)

当退化量分布不为位置-刻度族模型时，则要通过加速退化方程确定退化过程某参数的函数形式.假设ω(t)为该分布的时变参数，由加速退化方程可得lnω(t)=a(t)+b(t)φ(S)，这里令a(t)=β1+β2t，b(t)=β3，则可将加速退化方程改写为：

lnω(t)=β1+β2t+β3φ(S)

(10)

该加速退化方程同样适用于威布尔分布，对数正态分布等位置-刻度族分布类型.

则样本的联合对数似然函数为：

(11)

对上述联合对数似然函数采用极大似然估计可得到各个参数的估计值.

在获得每个参数的估计后，则可得到产品退化失效模式下的可靠度Rd(t,xt,S).

2.2.2 比例危险模型λ0(t)的参数估计 当突发失效与退化量相关时，引入基于比例危险的突发失效模型评估突发失效的失效率函数.也即，在不同应力下，对突发失效时间进行分布拟合，选择最优的分布类型作为基准分布类型，进一步确定基准失效率，从而表示出突发失效的失效率函数.

首先，记录产品在应力Si下发生突发失效时对应的退化量，记第l(l=1,2,…,Ni)个突发失效样品失效时对应的性能退化量为xl，该性能退化量对应的突发失效时间为tl.

选择合适的分布函数来拟合产品的突发失效时间，对于各个应力上突发失效时间分布类型的确定，首先确定其备选分布类型.然后采用分布假设检验，从备选分布中确定各应力下突发失效时间的最优分布类型.

(12)

式中：Fn(x)为经验分布函数；F(x)为假设的累积分布函数.

(13)

在确定各个应力下的最优分布类型后，选取各个应力下的最优分布类型为基准分布，对应的失效率即为基准失效率.

2.2.3 比例危险模型h(xt,S)的参数估计 在选取基准失效率函数后，突发失效的失效率函数可表示为：

λ(t,xt,S)=λ0(t){exp(α1+α2xt+α3φ(S))}

(14)

式中：λ0(t)为基准失效率；α1、α2、α3为待估参数；xt为突发失效发生时刻t对应的退化量.

根据失效率与可靠度的关系，可得到下列表达式：

(15)

当突发失效时刻t已知时，其对应的退化量xt随即已知.观察上式可知exp{α1+α2xt+α3φ(S)}与积分变量t无关，因此可将上式改写为

(16)

将(16)式两边进行变换，可得到如下线性关系：

(17)

yt=α1+α2xt+α3φ(S)

(18)

对(18)式进行线性回归即可得到α1、α2、α3的估计值.

而对(18)式进行线性回归的前提是yt为已知量，但式中R(t,xt,S)却是未知的，因此需要建立突发失效时间与可靠度的关系，下面将结合概率论中概率的基本概念以及古典概型等理论来建立产品突发失效时间与可靠度之间的联系，进而得到yt与t的关系.

2.2.4 突发失效可靠度函数预计 为了简化公式，以下均在应力水平Si下计算突发失效故障率并进行参数估计.

首先，从退化数据中判断发生突发失效的产品并记录其发生突发失效时对应的突发失效时间.记第l(l=1,2,…,Ni)个突发失效样品失效时对应的突发失效时间为tl，并将其进行升序排列，得到{tn}，其中n=1,2,…,Ni.

建立突发失效模式下突发失效故障率的函数关系，即当突发失效时间为tn时，突发失效故障率为：

(19)

由突发失效故障率与可靠度的关系，便可得到突发失效的可靠度.

综上，可得到基于比例危险的突发失效模型的失效率函数，进而结合可得到竞争失效下恒定应力加速退化产品的可靠度函数.

3 算例

现有一批激光器进行退化试验，该试验利用了GaAs激光器的工作电流，在80 ℃温度下随时间变化的百分比数据记录[14].共有22个样本参加试验，每隔500 h测试一次样本，至4 000 h为止.当激光器的工作电流增加至其额定电流的110%时，则认为该激光器发生退化失效，每个激光器监测8次.其中15个样本发生性能退化，另外7个样本发生突发失效.现利用本文中提出的方法对该型激光器进行可靠性评估.并与文献[8]中所给的可靠度曲线进行对比，验证本研究所建模型的正确性.

1)退化失效模型中的参数估计.

由文献[14]可知，1、2、4、6、7、8、10、11、13、14、16、18、19、20、22号样本为性能退化数据.图1显示的是该型激光器性能退化数据的变化趋势.

图1 激光器性能退化数据

通过分析可知各测量时刻的退化量服从威布尔分布故采用威布尔分布进行退化失效建模.由加速退化方程(10)式可得到恒定应力下威布尔分布函数为

(20)

式中：形状参数m为常数，尺度参数ln[η(t,S)]=β1+β2t+β3φ(S)，β1、β2、β3为未知参数.

鉴于该激光器试验仅在额定应力80 ℃下进行，所以可将上述尺度参数进行变换以满足试验，但并不影响该方法在恒定应力加速退化试验当中应用.由于温度应力S为常数，可将β1与β3φ(S)合并，令β0=β1+(1/80)β3，尺度参数则变形为ln[η(t)]=β0+β2t.因此，激光器退化量的分布函数为

(21)

则所有样本的对数似然函数可表示为

(22)

(23)

2)比例危险模型λ0(t)的参数估计.

由文献[14]可知，3、5、9、12、15、17、21号样本为突发失效数据，突发失效时间及其对应的退化量见表1.

表1 突发失效时间及其退化量

确定突发失效时间分布模型，备选分布类型有：正态分布、Weibull分布、对数正态分布、Gamma分布等.鉴于鉴于该激光器试验仅在额定应力80 ℃下进行，所以在应力80 ℃下采用AD检验法确定的最优分布即为基准分布.计算80 ℃的检验统计量AD值(见表2).

表2 80 ℃下检验统计量AD值

(24)

式中：φ(·)与Φ(·)分别指标准正态分布的密度函数和分布函数.

3)突发失效可靠度函数预计与比例危险模型h(xt,S)的参数估计.

将表1中的突发失效数据按突发失效时间排序，并按照1.2.4求出突发失效时间对应的突发失效故障率、可靠度，如表3.

表3 80 ℃下突发失效相关数据

鉴于该激光器试验仅在额定应力80 ℃下进行，所以可将(18)式进行变换以满足试验，但并不影响该方法在恒定应力加速退化试验中的应用.由于温度应力S为常数，可将α1与α3φ(S)合并，令α0=α1+(1/80)α3，则(18)式变形为

yt=α0+α2xt

(25)

由此得到基于比例危险的突发失效模型的失效率函数为

(26)

依据提出的可靠性建模方法，得到竞争失效情况下的可靠度为

(27)

图2 竞争失效可靠度曲线

由图2可见：在监测初期，激光器突发失数据较少，M1所得曲线略高于M2，但总的来说，两条可靠度曲线相差不大，与实际退化情况基本相符.在监测后期，突发失效数据增多，并且退化失效数据也开始增多.以4 000 h为例，文献[8]中给出的可靠度为0.65，本研究给出的可靠度为0.46，而通过频率去估计该时间点的可靠度为0.54，相较于文献[8]，本研究的可靠度值较为接近，更能客观地反映这一现象.相比M2，M1得到的可靠性评估精度略优且略保守，符合工程上的做法.除此之外，采用基于比例危险的突发失效模型，可以更方便地推广到竞争失效下恒定应力场合，适用性较文献[8]中的结果更好.

4 结论

为了解决恒定应力加速退化试验中竞争失效可靠性评估问题，采用加速退化方程和比例危险模型对该问题进行刻画，丰富了该方面的研究.

本研究基于突发与退化竞争失效模型，对突发失效和退化失效分别展开研究.

1)针对伪失效寿命方法的不足，采用退化量分布函数描述产品性能的退化过程.

2)对产品在性能退化中伴有的突发失效，采用基于比例危险的突发失效模型，通过AD检验法获得基准失效率函数λ0(t)，而对比例危险模型中的h(xt,S)，则通过非参数方法和线性回归得到h(xt,S)中的参数估计.

本文中提出的竞争失效场合下恒定应力加速退化可靠性模型，更加客观地反映了退化与突发的竞争失效过程，具有一定的工程应用价值，但对于含有多种退化失效和突发失效的竞争失效问题还需要继续研究.