变点模型下Weibull分布场合步加寿命试验的极大似然估计

Maximum Likelihood Estimation for Step-Stress Accelerated Life Testing with a Weibull Change-Point Model

摘要: 加速寿命试验是一种极为有效而经济的寿命试验方法，其中步加寿命试验又是比恒加寿命试验具有失效更快与参试样品数更少的优点，因此被广泛应用。此文中解决了失效机理改变时步加寿命试验失效数据的时间折算问题，给出了失效机理产生突变的Weibull分布定时和定数场合步进应力加速寿命试验的极大似然估计。

Abstract: Accelerated life testing is an extremely effective and economical method for life testing, among which step-stress accelerated life testing offers the advantages of faster failure occurrence and fewer test samples compared to constant-stress accelerated life testing, thus being widely adopted. This paper addresses the time transformation issue for failure data in step-stress accelerated life testing when the failure mechanism changes, and provides maximum likelihood estimation for step-stress accelerated life testing under Type-I and Type-II censoring, where the failure mechanism undergoes an abrupt change and is modeled by a Weibull distribution.

1. 引言

加速寿命试验是一种极为有效而经济的寿命试验方法，试验的目的是得到常应力下各种可靠性指标的数据。在进行加速寿命试验时，一般要确定失效机理不变的一个大概范围。然后在失效机理不变的范围内进行加速寿命试验，得到失效数据后进行参数估计，并利用加速方程外推到常应力水平，这类问题的理论已日趋成熟，并在实践中开始得到应用和推行，详见文[1]。然而，本文考虑一种更复杂的情况：在步进应力加速寿命试验过程中当应力水平达到某一点(称为变点)后产品的失效机理产生了突变(这在试验中极有可能发生)，而且突变点是不知道的，我们希望利用有突变点的这类失效数据来进行参数估计，以得到常应力下的各种可靠性指标。变点问题是近几十年统计界中的一个热门话题，在工业质量控制、经济、金融、医学、计算机及可靠性等领域有大量的应用。研究方法有极大似然法、Bayes方法、最小二乘法等，见文[2]。在可靠性统计方面，文[3]讨论到Weibull分布变点问题的Bayes估计，文[4]提出了加速退化试验变点模型的统计分析。文[5]给出了有一个突变点的Weibull分布场合恒加寿命试验的Bayes分析。本文中我们将续文[5]研究有一个突变点的Weibull分布场合步加寿命试验的参数估计。

2. Weibull分布场合步加试验安排与基本假定

加速寿命试验有恒定应力加速寿命试验、步进应力加速寿命试验和序进应力加速寿命试验三种。步进应力加速寿命试验，简称步加试验。它是先选定一组加速应力水平$S_1<S_2<\cdots <S_l$ ，它们都高于正常应力水平$S_0$ 。试验开始时是把一定数量的样品都置于应力水平$S_1$ 下进行寿命试验。经过一段时间，如${\tau }_1\left(h\right)$ 后，把应力提高到$S_2$ ，将未失效的样品在$S_2$ 下继续进行寿命试验。如此继续下去，直到有一定数量的样品发生失效为止。

在步加试验中，一个样品可能会遭遇若干个加速应力水平的考验。因此，相比恒加寿命试验，步加试验可使样品失效更快一些，并可以减少参试样品个数，且比序加试验更容易实施，正是因为有这些优点，步加试验被广泛应用，在组织与实施步加试验时应注意的事项请参考文献[1]。

2.1. Weibull分布场合步加试验的安排如下

1) 确定正常应力水平$S_0$ 和$l$ 个加速应力水平$S_1,S_2,\cdots ,S_l$ ，这些应力水平一般应满足如下关系式：$S_1<S_2<\cdots <S_l$ 。

2) 从一批产品中随机抽取$n$ 个样品进行步加试验，每步对未失效产品继续在下一级应力水平下进行试验。应力水平转换可以是定时转换，也可以是定数转换。

3) 对定时转换步加试验，事先要确定$l$ 个应力水平的持续时间${\tau }_1,{\tau }_2,\cdots ,{\tau }_l$ ，未失效产品在加速应力水平$S_i$ 下工作${\tau }_i$ 时间，不论失效多少，及时把应力水平提高到$S_{i+1}$ ，然后继续试验，直到在$S_l$ 下工作${\tau }_l$ 时间后才停止试验。

4) 对定数转换步加试验，事先要确定$l$ 个失效数：$r_1,r_2,\cdots ,r_l$ ，且$r_1+r_2+\cdots +r_l\le n$ ，然后要求在$S_i$ 下有$r_i$ 个失效发生时把应力水平提高到$S_{i+1}$ ，然后继续试验，直到在$S_l$ 下再有$r_l$ 个产品发生失效时就停止试验。

5) 设$n$ 个产品在$l$ 个加速应力水平$S_1,S_2,\cdots ,S_l$ 下分别失效$r_1,r_2,\cdots ,r_l$ 个，而在$S_i$ 下$r_i$ 个失效时间为

$t_{i1}\le t_{i2}\le \cdots \le t_{i{r}_i}\le {\tau }_i,i=1,2,\cdots ,l$ 。

当$t_{i{r}_i}={\tau }_i$ 时，就是定数转换步加试验数据；当$t_{i{r}_i}<{\tau }_i$ 时，就是定时转换步加试验数据。这里失效时间都是从应力水平提高到$S_i$ 时开始算起。

2.2. Weibull分布场合步加试验的基本假定

A1：在正常应力水平$S_0$ 和$l$ 个加速应力水平$S_1,S_2,\cdots ,S_l$ 下产品的寿命分布都服从威布尔分布$Wei\left(m_i,{\eta }_i\right)$ ，其分布函数为

$F_i\left(t\right)=1-\exp \left\{-{\left(\frac{t}{{\eta }_i}\right)}^{m_i}\right\},t>0,i=1,2,\cdots ,l$

其中诸$m_i>0$ 为形状参数，诸${\eta }_i>0$ 为特征寿命。

A2：假定只存在一个突变点$S^{\prime }$ ，且$S_k<S^{\prime }<S_{k+1}$ (此处$k$ 未知)，在应力水平$S_0$ 和$S_1<S_2<\cdots <S_k$ 下产品的失效机理不变，在$S_{k+1}<S_{k+2}<\cdots <S_l$ 下失效机理也相同，而在$S_k$ 和$S_{k+1}$ 之间失效机理产生了突变。由于威布尔分布的形状参数反映了失效机理，因此假定等价于$m_0=m_1=\cdots =m_k\ne m_{k+1}=\cdots =m_l$ 。

A3：产品的特征寿命${\eta }_i$ 与所施加速应力水平$S_i$ 之间满足加速模型

$\ln {\eta }_i=a+b\phi \left(S_i\right),i=0,1,2,\cdots ,l$ 。

其中$a$ 与$b$ 是待估参数，$\phi \left(S\right)$ 是$S$ 的已知函数，若应力为电压电流等，则取逆幂率模型，若应力为温度，则取阿伦尼斯模型。

A4：产品的残余寿命仅依赖于当时已累积失效部分和当时应力水平，而与累积方式无关。

2.3. 时间折算公式

步加试验统计分析的难点在于观察得到的失效数据不是寿命数据，因此如何把失效数据折算成寿命数据是解决此问题的关键。

由假定A4，产品在应力水平$S_i$ 下工作${\tau }_i$ 时间的累积失效概率$F_i\left({\tau }_i\right)$ 等于此产品在应力水平$S_i$ 下工作某一段时间${\tau }_{ij}$ 的累积失效概率$F_i\left({\tau }_{ij}\right)$ ，即

$F_i\left({\tau }_i\right)=F_i\left({\tau }_{ij}\right),i\ne j.$

把Weibull分布的分布函数代入上式可得

$1-\exp \left\{-{\left(\frac{{\tau }_i}{{\eta }_i}\right)}^{m_i}\right\}=1-\exp \left\{-{\left(\frac{{\tau }_{ij}}{{\eta }_j}\right)}^{m_j}\right\}$

根据假定A2，$m_0=m_1=\cdots =m_k\ne m_{k+1}=\cdots =m_l$ ，

当$m_i=m_j$ 时，有$\frac{{\tau }_i}{{\eta }_i}=\frac{{\tau }_{ij}}{{\eta }_j}$ ，

$\Rightarrow {\tau }_i=\frac{{\eta }_i}{{\eta }_j}{\tau }_{ij}={\tau }_{ij}{\text{e}}^{b\left[\phi \left(S_i\right)-\phi \left(S_j\right)\right]}$ 。 (1)

当$m_i\ne m_j$ 时，不妨设$m_j=q{m}_i$ ，有${\left(\frac{{\tau }_i}{{\eta }_i}\right)}^{m_i}={\left(\frac{{\tau }_{ij}}{{\eta }_j}\right)}^{m_j}$ ，

$\Rightarrow {\tau }_i=\frac{{\eta }_i}{{\eta }_j^q}{\tau }_{ij}^q={\tau }_{ij}^q{\text{e}}^{a\left(1-q\right)+b\left[\phi \left(S_i\right)-q\phi \left(S_j\right)\right]}$ (2)

(1) (2)两式即是两应力之间的时间折算公式。

公式(1)适用于失效机理未发生变化的阶段。它表明，在较低应力$S_{i-1}$ 下累积工作$t$ 时间所造成的损伤，等同于在较高应力$S_i$ 下工作一个更短的时间$t^{\text{*}}$ 。这一折算系数${\eta }_{i-1}/{\eta }_i$ 完全由加速模型决定。公式(2)则适用于失效机理发生突变后的阶段，此时不仅特征寿命比例变化，失效模式的变化(由形状参数$m_1$ 变为$m_2$ )也参与了时间折算，使得累积损伤的等效关系更为复杂，折算公式中引入了形状参数的幂次运算。

由前面的试验安排(5)，我们知道只有在应力$S_1$ 下的失效数据是寿命数据，而在应力$S_2,\cdots ,S_l$ 下得到的失效数据都不是寿命数据，因此我们先把应力$S_2,\cdots ,S_l$ 下的失效数据全部转换为应力$S_1$ 下的寿命数据。

由时间折算公式，当$i\le k$ 时，由于失效机理不变，采用公式(1)进行转换，从$S_i$ 的失效数据折算到$S_1$ 的寿命数据的折算公式为

$\begin{array}{c}{t}_{ij}^{*}\left(b,k\right)=\frac{{\eta }_1}{{\eta }_j}{t}_{ij}+\frac{{\eta }_1}{{\eta }_1}{\tau }_1+\frac{{\eta }_1}{{\eta }_2}{\tau }_2+\cdots +\frac{{\eta }_1}{{\eta }_{i-1}}{\tau }_{i-1}\\ =\frac{{\eta }_1}{{\eta }_j}{t}_{ij}+{\displaystyle \sum _{h=1}^{i-1}\frac{{\eta }_1}{{\eta }_h}{\tau }_h}\\ =t_{ij}{\text{e}}^{b\left[\phi \left(S_1\right)-\phi \left(S_i\right)\right]}+{\displaystyle \sum _{h=1}^{i-1}{\tau }_h{\text{e}}^{b\left[\phi \left(S_1\right)-\phi \left(S_h\right)\right]}},i=2,\cdots ,k;j=1,2,\cdots ,r_i.\end{array}$

当$i\ge k+1$ 时，由于存在突变点，在时间折算上要特别注意，在突变之前的折算用公式(1)，在突变之后折算时要用公式(2)，因此从$S_i$ 的失效数据折算到$S_1$ 的寿命数据的折算公式为

$\begin{array}{l}{t}_{ij}^{*}\left(a,b,q,k\right)=\frac{{\eta }_1}{{\eta }_i^q}{t}_{ij}^q+\frac{{\eta }_1}{{\eta }_1}{\tau }_1+\frac{{\eta }_1}{{\eta }_k}{\tau }_k+\frac{{\eta }_1}{{\eta }_{k+1}^q}{\tau }_{k+1}^q+\cdots \frac{{\eta }_1}{{\eta }_{i-1}^q}{\tau }_{i-1}^q\\ =\frac{{\eta }_1}{{\eta }_i^q}{t}_{ij}^q+{\displaystyle \sum _{h=1}^k\frac{{\eta }_1}{{\eta }_h}{\tau }_h}+{\displaystyle \sum _{h=k+1}^{i-1}\frac{{\eta }_1}{{\eta }_h^q}{\tau }_h^q}\\ =t_{ij}^q{\text{e}}^{a\left(1-q\right)+b\left[\phi \left(S_1\right)-\phi \left(S_i\right)\right]}+{\displaystyle \sum _{h=1}^{k-1}{\tau }_h{\text{e}}^{b\left[\phi \left(S_1\right)-\phi \left(S_h\right)\right]}}\\ +{\displaystyle \sum _{h=k+1}^{i-1}{\tau }_h^q{\text{e}}^{a\left(1-q\right)+b\left[\phi \left(S_1\right)-\phi \left(S_h\right)\right]}},i=k\text{+}1,\cdots ,l;j=1,2,\cdots ,r_i.\end{array}$

这样我们得到了一个容量为$n$ ，取自Weibull分布$Wei\left(m_1,{\eta }_1\right)$ 的定数截尾“样本”，截尾数为$r_1+r_2+\cdots +r_l=r\left(\le n\right)$ 。“样本”及其失效次序如下所示

$t_{11}\le t_{12}\le \cdots \le t_{1r_1}\le t_{21}^{*}\le \cdots \le t_{2r_2}^{*}\le \cdots \le t_{l1}^{*}\le \cdots \le t_{l{r}_l}^{*}.$

它实际上并不是真正的样本，因为它们含有未知参数$a,b,q,k$ 。

以下为了书写的方便，我们把此寿命数据简单地表示为

$t_1\le t_2\le \cdots \le t_r$ (3)

3. Weibull分布步加试验数据的极大似然估计

极大似然估计的数值求解算法可概述如下：

步骤1：对于给定的变点候选值$k$ (即假设突变发生在第$k$ 个应力步)，利用2.3节中的时间折算公式将所有失效数据折算至应力水平$S_1$ ，形成“样本”序列(3)。

步骤2：将折算后的序列视为来自Weibull分布$Wei\left(m_1,{\eta }_1\right)$ 的一个前r个次序统计量的观测值。令$m=m_1$ ，$\eta ={\eta }_1$ 。基于次序统计量的联合概率密度函数，可写出包含所有未知参数$\theta =\left(m,a,b,q,k\right)$ 的似然函数形式如(4)所示

$L\left(m,a,b,q,k\right)=\frac{n!}{\left(n-r\right)!}{\displaystyle \prod _{h=1}^r\frac{m}{{\eta }^m}{t}_h^{m-1}}\exp \left\{-{\left(\frac{t_h}{\eta }\right)}^m\right\}\cdot {\left\{\exp \left[-{\left(\frac{t_r}{\eta }\right)}^m\right]\right\}}^{n-r}$ (4)

对数似然函数为

$\begin{array}{c}f\left(m,a,b,q,k\right)=\ln L\left(m,a,b,q,k\right)\\ =\ln \frac{n!}{\left(n-r\right)!}+r\ln m-mr\ln \eta +\left(m-1\right){\displaystyle \sum _{h=1}^r\ln t_h}-{\displaystyle \sum _{h=1}^r{\left(\frac{t_h}{\eta }\right)}^m}-\left(n-r\right){\left(\frac{t_r}{\eta }\right)}^m\end{array}$ (5)

注意到：

$\ln \eta =\ln {\eta }_1=a+b\phi \left(S_1\right)$ $\Rightarrow \eta ={\text{e}}^{a+b\phi \left(S_1\right)}$

所以有$\frac{\partial \eta }{\partial a}={\text{e}}^{a+b\phi \left(S_1\right)}=\eta$ ，$\frac{\partial \eta }{\partial b}=\phi \left(S_1\right){\text{e}}^{a+b\phi \left(S_1\right)}=\eta \phi \left(S_1\right)$ 。

步骤3：固定$k$ ，通过对数似然函数(5)对其余参数$m,a,b,q$ 求偏导得到的似然方程组(6)如下：

$\{\begin{array}{l}\frac{\partial f}{\partial m}=\frac{r}{m}-r\ln \eta +{\displaystyle \sum _{h=1}^r\ln t_h}-{\displaystyle \sum _{h=1}^r{\left(\frac{t_h}{\eta }\right)}^m\ln \frac{t_h}{\eta }}-\left(n-r\right){\left(\frac{t_r}{\eta }\right)}^m\ln \frac{t_r}{\eta }=0\\ \frac{\partial f}{\partial a}=-mr+\left(m-1\right){\displaystyle \sum _{h=1}^r\left(\frac{1}{t_h}\frac{\partial t_h}{\partial a}\right)}-m{\displaystyle \sum _{h=1}^r{\left(\frac{t_h}{\eta }\right)}^{m-1}\frac{\frac{\partial t_h}{\partial a}-t_h}{\eta }}-\left(n-r\right)m{\left(\frac{t_r}{\eta }\right)}^{m-1}\frac{\frac{\partial t_r}{\partial a}-t_r}{\eta }=0\\ \frac{\partial f}{\partial b}=-mr\phi \left(S_1\right)+\left(m-1\right){\displaystyle \sum _{h=1}^r\left(\frac{1}{t_h}\frac{\partial t_h}{\partial b}\right)}-m{\displaystyle \sum _{h=1}^r{\left(\frac{t_h}{\eta }\right)}^{m-1}\frac{\frac{\partial t_h}{\partial b}-t_h}{\eta }}-\left(n-r\right)m{\left(\frac{t_r}{\eta }\right)}^{m-1}\frac{\frac{\partial t_r}{\partial b}-t_r}{\eta }=0\\ \frac{\partial f}{\partial q}=\left(m-1\right){\displaystyle \sum _{h=1}^r\left(\frac{1}{t_h}\frac{\partial t_h}{\partial q}\right)}-m{\displaystyle \sum _{h=1}^r{\left(\frac{t_h}{\eta }\right)}^{m-1}\frac{\frac{\partial t_h}{\partial q}}{\eta }}-\left(n-r\right)m{\left(\frac{t_r}{\eta }\right)}^{m-1}\frac{\frac{\partial t_r}{\partial q}}{\eta }=0\end{array}$ (6)

(6)是超越方程组，需通过如Newton-Raphson法等数值迭代求解，得到给定$k$ 下的各参数估计记为(7)

$\widehat{m}=m\left(k\right),\widehat{a}=a\left(k\right),\widehat{b}=b\left(k\right),\widehat{q}=q\left(k\right)$ (7)

步骤4：在变点参数$k$ 的合理整数取值范围内(通常为$2\le k\le l-1$ ，其中$l$ 为应力步数)进行一维搜索，寻找使

$f\left(k\right)=\ln \frac{n!}{\left(n-r\right)!}+r\ln \widehat{m}-\widehat{m}r\ln \widehat{\eta }+\left(\widehat{m}-1\right){\displaystyle \sum _{h=1}^r\ln {\widehat{t}}_h}-{\displaystyle \sum _{h=1}^r{\left(\frac{{\widehat{t}}_h}{\widehat{\eta }}\right)}^{\widehat{m}}}-\left(n-r\right){\left(\frac{{\widehat{t}}_r}{\widehat{\eta }}\right)}^{\widehat{m}}$

达到最大的$\widehat{k}$ ，即为变点位置的估计。

步骤5：将$\widehat{k}$ 代回(7)式，即可得所有参数的极大似然估计：$\widehat{k},\widehat{m}=m\left(k\right),\widehat{a}=a\left(k\right),\widehat{b}=b\left(k\right),\widehat{q}=q\left(k\right)$ 。

进而利用加速方程外推可得到常应力$S_0$ 下的各种可靠性指标。

4. 仿真例子

为验证本文提出的变点模型下Weibull分布场合步加寿命试验极大似然估计方法的有效性与可靠性，为此设计并实施了蒙特卡洛仿真研究。通过模拟不同加速条件下的步加寿命试验数据，评估本文方法在变点检测、参数估计及正常应力下可靠性指标预测方面的性能。

4.1. 基本试验方案

仿真模拟生成40个寿命服从Weibull分布的产品，进行4步定数截尾步加试验。具体设置如下：

1) 样本量：$n=40$ ；

2) 应力步数：$l=4$ ；

3) 截尾方案：定数截尾，每步失效数$r_i=5$ (即$r_1=r_2=r_3=r_4=5$ )，总失效数$r=20$ ；

4) 加速模型：采用阿伦尼斯模型$\ln {\eta }_i=a+b\frac{1}{T_i}$ ，其中取：$a=-17.636$ ；$b=8000$ ；

5) 温度应力水平设定(绝对温度，单位：K)：

正常应力：$T_0=298\text{K}\left(S_0=25˚\text{C}\right)$ ；

加速应力1：$T_1=333\text{K}\left(S_1=60˚\text{C}\right)$ ；

加速应力2：$T_2=353\text{K}\left(S_2=80˚\text{C}\right)$ ；

加速应力3：$T_3=373\text{K}\left(S_3=100˚\text{C}\right)$ ；

加速应力4：$T_4=393\text{K}\left(S_4=120˚\text{C}\right)$ ；

6) 变点位置：失效机理在应力水平$S_2$ 和$S_3$ 之间发生改变，即变点$k=3$ ；

7) 形状参数设定：变点前形状参数：$m_1=1.5$ ，变点后形状参数：$m_2=3.0$ 。

4.2. 仿真结果与分析

Table 1. Parameter estimation results for step-stress accelerated life testing under the Weibull distribution

表1. Weibull分布步加寿命试验各参数估计结果

参数	真实值	估计均值	标准差	偏差	RMSE	相对误差
$a$	−17.636	−17.521	0.854	0.115	0.861	4.88
$b$	8000	7963.2	312.6	−36.8	314.8	3.93
$m_1$	1.5	1.482	0.126	−0.018	0.127	8.47
$m_2$	3.0	2.974	0.214	−0.026	0.216	7.20
${\eta }_0$	10000.0	9876.4	385.2	−123.6	405.8	4.06
$k$	3	2.96	0.197	−0.04	0.198	-

从表1可以看出，所有参数的估计偏差均接近0，表明本文提出的极大似然估计方法具有良好的无偏性，形状参数的估计精度相对较低，这与Weibull分布形状参数估计的固有难度相符，但仍处于可接受范围。变点准确率达到93.0%，表明本文方法能有效识别失效机理的突变点。

5. 结论–展望

本文针对失效机理可能发生突变的实际情况，研究了Weibull分布下步进应力加速寿命试验的统计分析问题。通过引入变点模型，导出了分段形式的时间折算公式，并构建了相应的联合似然函数，提出了基于对数似然和数值优化的极大似然估计方法。模拟研究表明，该方法能够有效估计变点位置及各项可靠性参数，为处理复杂失效数据的加速寿命试验提供了实用的解决方案。

本研究仍存在一些值得进一步探讨的方向：首先，本文方法依赖于Weibull分布的假设，未来可研究其他分布族(如对数正态分布、广义Gamma分布)下的变点模型。其次，极大似然估计在小样本下的性质可能不佳，可考虑采用贝叶斯方法，引入先验信息以改善估计的稳定性。最后，本文假定变点位置与应力水平挂钩且只有一个，对于多个变点或变点位置与累积损伤量相关的情形，将是更具挑战性的后续研究课题。

基金项目

广西民族师范学院2022年度高层次人才项目：大规模Minimax优化问题有效算法研究(2022GCC004)。