加速寿命试验(accelrated life testing, ALT)通过增大试验件所受到的应力, 如温度、电压、压力等, 以便在更短的时间内获得相应的故障数据, 是获得高可靠性、长寿命产品失效数据的一种有效途径。茆诗松^[1]与Voiculescu等^[2]学者已经对加速寿命试验的数学模型进行了大量研究。在现有研究中, 通常仅考虑单一应力的影响, 但在工程实际中, 受试产品的寿命会受到多种环境应力综合效果影响。另外, 加速寿命试验要求产品的失效机理在试验过程中不发生改变^[3], 当某一应力水平超过产品承受上限时, 产品失效机理就会发生变化, 所以单一应力水平不宜增加过高, 这样试验时间的缩短就会受到限制。因此, 为了增加试验效率, 提高试验结果的准确性, 有必要针对多应力加速寿命试验及其统计分析进行研究。张春华等^[4-5]提出一种步降应力加速寿命试验方法, 该方法在损耗型长寿命产品试验中可显著提高试验效率。孙天宇等^[6]在双应力交叉步加试验下研究了寿命服从威布尔分布部件可靠性指标的估计, 给出参数的Bayes估计。寇海霞等^[7]提出一种双应力同步步降加速寿命试验方法来提高试验效率, 并针对威布尔分布产品进行仿真分析。Chen等^[8]以气缸为对象开展双应力交叉步降加速寿命试验, 获得了气缸的可靠性指标并验证了加速模型的正确性。

以上研究均假设产品只有一种失效模式, 然而在工程实际中, 产品的失效往往是由多种失效模式导致, 因此有必要引入竞争失效模型, 在该模型中, 所观测到的失效数据由失效时间和失效机理标志两部分组成。竞争失效模型在医学、工程、计量经济学等很多领域得到了广泛应用, 目前, 国内外针对竞争失效模型在加速寿命试验中的应用已有很多研究。Han等^[9]研究了指数分布下步进应力竞争失效模型, 给出参数的极大似然估计(maximum likelihood estimation, MLE)并构建了置信区间。张春芳等^[10]研究了指数-威布尔分布的竞争失效模型, 分别应用极大似然估计方法和Bayes方法得到参数的点估计和区间估计。Alam等^[11]在Ⅰ型和Ⅱ型截尾下讨论了Fréchet分布竞争失效部分加速寿命试验的统计分析问题。文献[12-14]分别分析了寿命服从威布尔分布、Lindley分布以及Kumaraswamy分布时的竞争失效模型。

在引入双应力竞争失效模型的同时, 为了进一步节省试验时间和成本, 还需要考虑截尾方案。最常用的截尾方案是Ⅰ型截尾(定时截尾)和Ⅱ型截尾(定数截尾), Epstein^[15]首先提出综合Ⅰ型和Ⅱ型截尾的混合截尾方案。然而, 这几种截尾方案只在试验截止时移走未失效的试验件, 使得试验缺少灵活性。在之后被提出的逐步截尾方案中, 每次失效发生时, 均会随机移走一定数量的未失效试验件。Kundu等^[16]结合逐步Ⅱ型截尾和混合截尾方案, 提出一种逐步Ⅱ型混合截尾方案, 避免了逐步截尾试验时间过长的问题。然而在这种方案中, 可能会出现失效样本数量过小的问题, 从而对数据的统计分析产生影响。于是, Ng等^[17]提出了自适应Ⅱ型逐步截尾(adaptive type-Ⅱ progressive censoring, AT-Ⅱ PCS), 既能控制试验时间又能确保在试验中获得失效样本数量。近年来, AT-Ⅱ PCS的应用得到了广泛研究, 文献[18-20]分别讨论了广义Pareto分布、广义指数分布和Weibull分布的参数估计问题。Ren等^[21]将AT-Ⅱ PCS和竞争失效模型相结合, 并通过数值模拟与试验数据验证了所提模型的正确性。

AT-Ⅱ PCS能够有效节省试验时间和成本, 但由于模型和计算的复杂性, 针对其在加速寿命试验中的应用研究仍然较少。因此, 本文在AT-Ⅱ PCS下考虑寿命服从Weibull分布的独立竞争失效产品, 建立起双应力恒加试验的统计分析方法。建立似然函数方程, 并采用Newton-Raphson迭代和渐进似然理论得到未知参数的MLE和渐进置信区间。采用Gibbs与M-H相结合的算法得到参数的Bayes估计, 并采用Monte Carlo(MC)方法构建参数的最大后验密度可信区间。最后, 通过数值模拟方法验证模型的正确性, 并比较在不同的试验方案下MLE与Bayes方法的参数估计效果。

1 模型描述 1.1 自适应Ⅱ逐步型截尾下双应力恒加试验

假设有2组应力水平S¹, S²满足

(1)

式中: k和l分别表示2组加速应力的水平数; S₀¹, S₀²分别为2组应力的正常应力水平。第一组应力中的第i个加速应力水平和第二组应力中的第j个加速应力水平的组合(S¹_i, S²_j)称为水平组合, 简记为(i, j), i=1, 2, …, k, j=1, 2, …, l。水平组合共可以有kl个, 将n个参试样品分为kl组, 记每组样品数量为n_ij, 若只选择其中的部分水平组合进行恒加试验, 则样本量相应减少, 为了叙述方便, 假设有kl个水平组合。给定每组的预期失效数m_ij, 预期试验时间τ_ij和逐步移走方案(R_{ij, 1}, R_{ij, 2}, …, R_{ij, mij})。试验过程中, 在第1个样品失效时刻t_{ij, 1}, 从剩余n_ij-1个未失效样品中随机移走R_{ij, 1}个样品, 在第2个样品失效时刻t_{ij, 2}, 从剩余n_ij-R_{ij, 1}-2个未失效样品中随机移走R_{ij, 2}个样品, 以此类推, 直到第m_ij个失效时刻t_{ij, mij}, 移走剩下所有的未失效产品, 即R_{ij, mij}=n_ij- 。试验会出现以下2种情形:

CaseⅠ  若t_{ij, mij} < τ_ij, 试验结束, 移走方案不变, 这种情况等同于普通的逐步Ⅱ型截尾方案。

CaseⅡ  若t_{ij, J} < τ_ij < t_{ij, J+1}, J=0, 1, 2, …, m_ij-1且t_{ij, 0}=0, 令R_J+1=R_J+2=…=R_mij-1=0, R_mij=n_ij-m_ij- , 试验继续, 也就是将移走方案(R_{ij, 1}, R_{ij, 2}, …, R_{ij, mij})变为(R_{ij, 1}, R_{ij, 2}, …, R_{ij, J}, 0, …, 0, R^*_{ij, mij}=n_ij-m_ij- )。

自适应Ⅱ型逐步截尾方案能保证每组最终的失效样本数为m_ij, 并且当试验时间超过τ_ij后, 试验将尽可能快结束。

1.2 基本假设

1) 共n个参试样品分为kl组, 在水平组合(i, j)下有n_ij个样品进行试验且 。

2) 每个样品的失效由m种独立失效机理的其中一种造成, 记水平组合(i, j)下m种失效机理的发生时间分别为T_{ij, 1}, T_{ij, 2}, …, T_{ij, m}, 则样品的失效时间为min{T_{ij, 1}, T_{ij, 2}, …, T_{ij, m}}。

3) 水平组合(i, j)下, T_{ij, 1}, T_{ij, 2}, …, T_{ij, m}均服从两参数Weibull分布W(η_{ij, u}, β_u), 其中β_u为形状参数, η_{ij, u}为尺度参数, 分布函数(cumulative distribution function, CDF)和密度函数(probability density function, PDF)分别如(2)~(3)式所示。

(2)

(3)

4) 在各水平组合下产品失效机理一致, 即β_11u=β_12u=β_21u=…=β_u(u=1, 2, …, m)。

5) 第u个失效机理的加速模型为

(4)

式中: c_0u, c_1u, c_2u, u=1, 2, …, m为待估参数; φ₁, φ₂为已知函数。

1.3 竞争失效模型

基于1.2节中的假设, 在水平组合(i, j)下进行AT-Ⅱ PCS双应力恒加试验, 观测到的失效数据包含失效时间和失效机理标志, 即(t_{ij, 1}, α_{ij, 1}), (t_{ij, 2}, α_{ij, 2}), …, (t_{ij, mij}, α_{ij, mij}), 其中t_{ij, 1}, t_{ij, 2}, …, t_{ij, mij}为次序统计量, α_{ij, s}∈{1, 2, …, m}。令

则水平组合(i, j)下由第u个失效机理造成的失效数为 。本文所提出的方法很容易推广到m>2的情形, 为方便表述, 考虑失效机理个数m=2。在2种独立失效模式下所观测到的失效时间t_ij的CDF和PDF分别为

(5)

(6)

产品的失效由失效模式1造成的概率为

(7)

由(7)式可得 。

2 极大似然估计 2.1 似然函数和点估计

假设每个水平组合下至少有一个失效由第u个失效机理造成, 其似然函数为

(8)

式中

将(2)~(3)式代入(8)式得

(9)

第u个失效机理的完整似然函数和全部失效数据的完整似然函数分别为

(10)

(11)

式中

参数θ _u=(β_u, η_{11, u}, η_{12, u}, …, η_{kl, u})的MLE记为 , 可直接通过(10)式求得。将(10)式取对数得

(12)

(12) 式分别对参数β_u, η_{ij, u}求偏导并令偏导数为0, 可得以下方程组

(13)

(14)

采用Newton-Raphson迭代法求解(13)和(14)式即可得到参数的MLE。

定理1   存在且唯一。

证明  由(13)式可得

对任意i, j, 均满足t_{ij, mij}≤η_{ij, u}时

对任意i,j , 存在t_{ij, mij}>η_{ij, u}时

故 有解。

(12) 式对β_u求二阶导得

故 单调递减, (13)式有唯一解, 证毕。

由(14)式得η_{ij, u}的极大似然估计 为

(15)

η_{ij, u}, β_u>0, 因此对于确定的β_u, 存在且唯一。

2.2 渐进置信区间

由于参数极大似然估计的精确分布推导复杂, 计算困难, 难以获得参数的精确区间估计, 本节采用渐进似然理论构建参数的渐进置信区间。通过对参数的观测Fisher信息矩阵求逆获得参数的渐进方差-协方差矩阵。

参数 的观测Fisher信息矩阵为

(16)

参数MLE的方差-协方差矩阵为

矩阵 对角线上的元素为对应参数的方差, 即

对给定的置信水平100(1-γ)%, 参数 的渐进置信区间为

z_α/2是标准正态分布的α/2分位数。

3 Bayes估计 3.1 后验分布

Bayes方法能结合先验信息与样本信息对未知参数进行估计, 是传统估计方法的一种有效替代方法。本节采用Bayes方法, 在平方损失函数下获取参数 的点估计和区间估计。根据文献[21]所述, 令λ_{ij, u}=1/η_{ij, u}^β_u, u=1, 2, 选用无信息先验作为参数λ_{11, u}, λ_12u, …, λ_{kl, u}的先验分布, 选取Gamma分布作为参数β_u的先验分布, 形式为：

其中, 超参数a_u, b_u均为非负常数。于是参数 ω_u=(β_u, λ_{11, u}, λ_{12, u}, …, λ_{kl, u})联合先验分布表示为

(17)

根据(10)和(17)式得到参数联合后验密度函数

(18)

式中

由于后验分布的形式十分复杂, 难以求得Bayes估计的显式表达式, 故而采用Gibbs抽样算法来解决这一问题。

3.2 Gibbs抽样算法

由(18)式可以得到(β_u, λ_{11, u}, λ_{12, u}, …, λ_{kl, u})等参数的全条件密度

(19)

(20)

式中

β_u的全条件密度是对数凹函数, 证明详见文献[22], 并且由(20)式可知λ_{ij, u}的全条件密度服从参数为N_{ij, u}和T_ij(β_u)的Gamma分布。通过Gibbs抽样得到后验分布随机样本的流程如下

Step1  给定参数的初值(β⁽⁰⁾_u, λ_{11, u}⁽⁰⁾, …, λ_{kl, u}⁽⁰⁾), 可以选择MLE的结果作为初值;

Step2  设置迭代步数q=0;

Step3  采用Metropolis-Hasting(M-H)算法(或自适应拒绝抽样)从 中生成β_u^(q+1);

Step4  从 中生成λ_{ij, u}^(q+1), 并计算 。设置q=q+1;

Step5  重复Step3~4 N_g次, 并且给定预烧期A(burnin period), 最终得到平方损失下参数的Bayes估计

(21)

根据文献[23]提出的MC方法, HPD可信区间构建方法如下:

1) 将Gibbs抽样中得到的 升序排列, 得到 。

2) 计算100(1-α)%的可信区间[β_u^q, β_u^{q+(1-α)·N_g}], 。

3) 从 中选择最短的区间作为 的HPD可信区间。

4 产品特征寿命估计

正常应力下参数η_{00, u}的估计值由(22)式计算得到

(22)

再利用最小二乘法获得参数 的估计值。

令

(23)

。利用最小二乘法得到 c_u的估计为

(24)

5 数值模拟与数据分析

本节采用数值模拟评估所提的MLE与Bayes估计方法的表现。考虑以温度和电压作为加速应力, 每组应力有4个加速应力水平, 分别为

正常应力水平为T₀=338 K, V₀=125 V。应力水平组合为{(1, 2), (3, 1), (2, 4), (4, 3)}, 每个水平组合下投入的样本量n_ij, 预定失效数m_ij, 移走方案R_ij和预期时间τ_ij见表 1, 每组方案分别在预期时间为τ₁, τ₂时进行AT-Ⅱ PCS恒加试验。

假定产品有2种失效模式, 加速方程为

由加速方程得尺度参数分别为

形状参数的值取为β₁=1.5, β₂=1.6。在应力水平组合(i, j)下, 产品的失效数据由以下步骤获得:

1) 根据方案给定的样本量n_ij, 预定失效数m_ij, 移走方案R_ij, 从PDF为(6)式的竞争失效Weibull分布中抽取Ⅱ型逐步截尾样本t_{ij, 1}, t_{ij, 2}, …, t_{ij, mij}, 详细过程见文献[24];

2) 根据预期时间τ, 找出满足t_{ij, J} < τ_ij < t_{ij, J+1}的J, 舍弃t_{ij, J+2}, t_{ij, J+3}, …, t_{ij, mij};

3) 从截断分布f_ij(t_ij)/[1-F_ij(t_{ij, J+1})]中抽取样本量为 , 失效数为(m_ij-J-1)的失效数据作为t_{ij, J+2}, t_{ij, J+3}, …t_{ij, mij};

4) 根据(7)式得N_{ij, 1}~B(m_ij, q_ij), 由此抽取一组失效机理标志(α_{ij, 1}, α_{ij, 2}, …, α_{ij, mij})。

为了克服随机性, 整个模拟过程重复1 000次。计算Bayes估计时, 形状参数β_u的先验取Gamma分布, 尺度参数η_ij取无信息先验, Gibbs算法中迭代次数N_g=2 000预烧期A=200, 提议分布为正态分布。计算得到的点估计平均值, 估计的均方根误差(root mean squared error, RMSE)和显著性水平α=0.05时的平均置信区间长度(mean interval length, MIL)如表 2所示。RMSE按(25)式计算

(25)

式中: 表示待估参数第i次计算的估计值; θ_u表示待估参数的真值。

图 1给出了方案3中当τ=τ₂时Gibbs采样算法中的迭代路径图, 各参数的值在均值附近上下波动, 表明算法中的马尔科夫链已经充分混合, 并收敛到平稳分布。

从表 2中可以得到以下几点结论:

1) 随着样本数量增大, 或者m_ij/n_ij值增大时, 2种估计方法对尺度参数η_ij的点估计均值更加接近真值, 同时RMSE和平均置信区间长度减小, 说明失效样本数量的增加能够提高估计的准确性。

2) 在同样的(n_ij, m_ij)下, 当预期时间从τ₁增加到τ₂时, 2种方法对于尺度参数η_{ij, 1}, η_{ij, 2}的点估计均值更接近真值, 同时RMSE与平均置信区间长度减小。这是由于τ=τ₁时试验结束得快, 使得观测到的失效数据较τ=τ₂时更小, 因此估计误差也相对较大。

3) 2种估计方法对形状参数β₁, β₂的估计效果在各种方案下没有显著差异, 即失效样本的增加不会使形状参数的估计更准确, 在实际应用中应利用设计经验或是先验信息对形状参数的估计做出更多改进与修正。

4) 当m_ij/n_ij较小(如方案1和方案3)且试验预期时间τ较小时, MLE得到的结果波动较大, 特别是试验得到一些极端数据时(例如最后一个失效样本的失效时间很小), 在采用Newton-Raphson迭代求解MLE时会得到异常值甚至出现迭代不收敛的情况, 此时Bayes估计的效果优于MLE。而当试验预期时间τ较大或m_ij/n_ij较大时, 2种估计方法的效果均能够接受, 但整体上MLE对于尺度参数的估计值较Bayes估计更接近真值, MSE和MIL也较小, 此类情况下MLE优于本文所给先验组合下的Bayes估计。

表 3给出了基于上述模拟结果进行拟合得到的加速方程中未知系数的估计值, 根据加速方程外推得到常应力下的特征寿命估计值以及使用时长分别为500, 1 500和3 000 h时产品的可靠度估计值。方案1中τ=τ₁时MLE的估计结果误差较大, 故省略该情形下的估计值。从表 3可以发现各水平组合(i, j)下失效样本数量以及试验预期时间τ对于估计准确性有显著影响。在各方案下MLE与Bayes估计对可靠度的估计值均呈现略微偏大的趋势, 随着m_ij/n_ij和τ的增加, 估计的结果逐渐接近于真值。

6 结论

本文提出了寿命服从Weibull分布的竞争失效产品在自适应Ⅱ型逐步截尾下双应力恒加试验的统计分析方法。首先建立似然函数方程, 通过Newton-Raphson迭代得到Weibull分布中未知参数的MLE, 基于渐进似然理论构建参数的渐进置信区间。其次, 在形状参数的先验分布为Gamma分布, 尺度参数取无信息先验的条件下获得参数的后验密度函数, 采用Gibbs采样与M-H抽样结合的算法在平方损失下计算参数的Bayes估计, 并基于MC方法建立参数的HPD可信区间。最后通过数值模拟验证了所提方法的可行性, 结果表明: ①预期试验时间τ对于参数估计效果有显著影响, 特别是当失效样本数m_ij占投入样本数量n_ij的比例较小时, 过小的τ会使MLE产生较大波动, 此时建议优先考虑使用Bayes估计。②当失效样本数量较大时, 2种估计均有良好的表现且MLE的估计效果略优于本文所给先验组合下的Bayes估计。