利用性能退化数据进行产品的可靠性评估是当前较为通用的做法。Nelson^[1]、Lu与Meeker等^[2-3]在基于性能退化数据的可靠性评估方面做了很多工作；其中使用较多的数学模型有回归模型和Wiener过程模型。对于退化过程存在随机波动的产品，由于回归模型的方差是常量，无法对退化过程中退化量的波动性进行建模，只能对大样本下不同样本的回归方程系数进行建模来评估产品总体的退化失效随机性和可靠性；而Wiener过程数学模型是由Brown运动和线性漂移项合成的，对产品的性能退化过程波动性具有很好的建模效果，因此常用于具有波动性的性能退化过程的建模。

彭宝华等^[4]采用Wiener过程对电容器的电容值耗损量建立了数学模型，并对其可靠性进行了评估；朱德馨等^[5]将Wiener过程模型中的轴承退化量初值和漂移系数随机化，结合Bootstrap方法对极小样本下高速列车轴承的可靠性进行了建模；王明磊等^[6]采用二元Wiener过程对小样本的电主轴可靠性进行了分析；Wang^[7]基于Wiener过程对桥横梁的可靠性进行了研究；Tang等^[8]采用Wiener过程对发光二极管的亮度建立了模型，并评估其可靠性；Nicolai等^[9]利用Wiener过程分析了钢构表面锈迹的扩散；Whitmore^[10]利用Wiener过程建立了晶体管增益的退化轨迹模型。

基于Wiener过程的可靠性评估有多种不同的实现方法。如果假设未知参数均为常数，根据所有试验数据构造似然函数来估计未知参数，虽然计算简单但是无法体现不同样本之间的差异；但若假设所有参数均服从正态分布，然后使用所有数据构造似然函数估计未知参数，将会使得评估过程异常复杂。本文先根据样本数据估计退化模型参数，然后根据参数的估计结果假设其恒定或服从正态分布，根据以上假设估计模型参数，最后进行转台可靠性估计。

1 Wiener过程退化模型

Wiener退化过程的数学表达式为

$ y\left( {t;v,{\delta ^2}} \right) = a + vt + \delta B\left( t \right). $

(1)

式（1）中：y(t; v，δ²)是时刻t退化量的测量值；a是退化量的初始值；v是漂移系数；δ是扩散系数。B(t)是标准Brown运动，满足以下性质：

1）独立增量性。每个时刻的性能参数的观测值都服从独立同分布。

2）稳定增量和正态性。研究样本属于退化失效型产品，即某个性能参数随着时间的推移表现为单调上升或者下降的趋势，且性能增量服从正态分布。

3）B(0)=0，从而B(t_i)~N(0，t_i)。

产品的某个性能指标达到某个特定水平时失效，则其失效时间T可以定义为产品的性能退化量X首次达到失效阈值D_f所用的时间，即

$ T = \inf \left\{ {t\left| {X\left( t \right) = {D_{\rm{f}}},t \ge 0} \right.} \right\}. $

(2)

寿命T服从逆Gauss分布，其概率密度函数为^[11]

$ f\left( t \right) = \frac{{{D_{\rm{f}}} - a}}{{\sqrt {2{\rm{ \mathsf{ π} }}{\delta ^2}{t^3}} }} \cdot \exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - a - vt} \right)}^2}}}{{2{\delta ^2}t}}} \right). $

(3)

产品性能退化量呈递增趋势时，即D_f>a，产品的可靠度函数为^[12]

$ \begin{array}{*{20}{c}} {R\left( t \right) = \mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - a - vt}}{{\delta \sqrt t }}} \right) - }\\ {\exp \left( {\frac{{2v\left( {{D_{\rm{f}}} - a} \right)}}{{{\delta ^2}}}} \right)\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - a} \right) - vt}}{{\delta \sqrt t }}} \right).} \end{array} $

(4)

1.1 模型参数估计

由Wiener过程模型（1）有：

$ \Delta {y_i} = v\Delta {t_i} + \delta \Delta {B_i},i = 1,2, \cdots ,N. $

(5)

式（5）中：Δy_i=y_i-y_i-1，Δt_i=t_i-t_i-1，ΔB_i=B_i-B_i-1，i=1，2，…，N。由于ΔB_i~N(0，Δt_i)，则：

$ \Delta {y_i} \sim N\left( {v\Delta {t_i},{\delta ^2}\Delta {t_i}} \right),i = 1,2, \cdots ,N. $

(6)

根据式（6），且由于Wiener过程的退化增量具有独立性，则其似然函数为

$ L\left( {v,\delta } \right) = \sum\limits_{i = 1}^N {\frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}{\delta ^2}\Delta {t_i}} }}{{\rm{e}}^{ - \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{2{\delta ^2}\Delta {t_i}}}}}} . $

(7)

对式（7）取对数，得到

$ \begin{array}{*{20}{c}} {\ln L\left( {v,\delta } \right) = }\\ {\sum\limits_{i = 1}^N {\left[ { - \ln \sqrt {2{\rm{ \mathsf{ π} }}\Delta {t_i}} - \ln \delta - \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{2{\delta ^2}\Delta {t_i}}}} \right]} .} \end{array} $

(8)

分别对v和δ求偏导，得到：

$ \frac{{\partial \ln L}}{{\partial v}} = \sum\limits_{i = 1}^N {\frac{{\Delta {y_i} - v\Delta {t_i}}}{{{\delta ^2}\Delta {t_i}}}} = 0, $

(9)

$ \frac{{\partial \ln L}}{{\partial \delta }} = \sum\limits_{i = 1}^N {\left( { - \frac{1}{\delta } + \frac{{{{\left( {\Delta {y_i} - v\Delta {t_i}} \right)}^2}}}{{{\delta ^3}\Delta {t_i}}}} \right)} = 0. $

(10)

由此可知:

$ \hat v = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{\Delta {y_i}}}{{\Delta {t_i}}}} , $

(11)

$ \hat \delta = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{{{\left( {\Delta {y_i} - \hat v\Delta {t_i}} \right)}^2}}}{{\Delta {t_i}}}} } . $

(12)

如果将模型中的每个未知参数都假设为服从正态分布，将导致计算异常困难。虽然不同样本的性能退化量初值、漂移系数和扩散系数均不相同，但是根据不同参数的差异性，有些参数可以近似认为其保持不变，因此只需要将差异较大的未知参数假设为服从正态分布的随机变量即可，从而降低计算的复杂度。进行参数性质假定时，可以参考其变异系数${C_{\rm{v}}} = \left( X \right) = \frac{{\sqrt {{\rm{Var}}\left( X \right)} }}{{E\left( X \right)}}$^[11]。其中：Var(X)表示参数的标准差, E(X)表示参数的平均值。将变异系数较小的参数假定为恒定值，将变异系数较大的参数假定为服从正态分布。

1.2 模型参数分布假设

假设样本初值服从正态分布，漂移系数和扩散系数为恒定值，即a~N(a₀，σ²_a)，$\widehat {{a_0}} = \frac{1}{N}\sum\limits_{i = 1}^N {{a_i}} $，$\widehat \sigma _a^2 = \frac{1}{{N-1}} \cdot {\sum\limits_{i = 1}^N {\left( {{a_i}-\widehat {{a_0}}} \right)} ^2}$，$\hat v = \frac{1}{N}\sum\limits_{i = 1}^N {{v_i}} $，$\hat \delta = \frac{1}{N}\sum\limits_{i = 1}^N {{\delta _i}} $。

当漂移系数和扩散系数恒定、初值a~N(a₀，σ²_a)时，结合式（4），产品的可靠度函数为

$ \begin{array}{*{20}{c}} {R\left( t \right) = \int {\left( {\mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt}}{{\hat \delta \sqrt t }}} \right) - } \right.} }\\ {\exp \left( {\frac{{2\hat v\left( {{D_{\rm{f}}} - \widehat {{a_0}}} \right)}}{{{{\hat \delta }^2}}}} \right),}\\ {\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - \widehat {{a_0}}} \right) - \hat vt}}{{\hat \delta \sqrt t }}} \right)p\left( a \right){\rm{d}}a.} \end{array} $

(13)

其中：$p\left( a \right) = \frac{1}{{\widehat {{\sigma _a}}}}\varphi \left( {\frac{{a-\widehat {{a_0}}}}{{\widehat {{\sigma _a}}}}} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}\widehat \sigma _a^2} }}{{\rm{e}}^{-\frac{{{{\left( {a-\widehat {{a_0}}} \right)}^2}}}{{2\widehat {\sigma _a^2}}}}}$，积分上下限取为$\left( {\widehat {{a_0}}-3\widehat {{\sigma _a}}, \widehat {{a_0}} + 3\widehat {{\sigma _a}}} \right)$，修正因子$k = \frac{1}{{2\mathit{\Phi }\left( 3 \right)-1}}$。

概率密度函数为

$ \begin{array}{*{20}{c}} {f\left( t \right) = \int {\left\{ {\frac{{{D_{\rm{f}}} - \widehat {{a_0}}}}{{\sqrt {2{\rm{ \mathsf{ π} }}{{\hat \delta }^2}{t^3}} }} \cdot } \right.} }\\ {\left. {\exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt} \right)}^2}}}{{2{{\hat \delta }^2}t}}} \right)} \right)p\left( a \right){\rm{d}}a.} \end{array} $

(14)

类似地，还可以假设漂移系数或扩散系数符从正态分布，推导过程同上。

2 试验验证

采用某型号数控转台进行可靠性试验，在不同转台角度θ下得到的转台回转定位精度误差退化量见表 1。回转定位精度误差退化量的退化趋势见图 1。

从图 1可以看出，数控转台的回转精度误差退化量随时间增加不仅有上升趋势也有波动趋势，反映了其性能逐渐退化，且由于试验环境、载荷以及测量误差等因素的影响，在退化过程中性能存在一定的随机波动。这种情况适合采用Wiener过程对其进行建模评估。

由1.2节的计算过程并结合表 1中的数据，得出Wiener模型的各参数，见表 2。

根据表 2数据，计算各个参数的平均值、标准差以及变异系数，见表 3。

由表 3可以看出，漂移系数和扩散系数比较稳定，可以将它们设为恒定值，而退化量初值的变异系数较大，可假设它服从正态分布，即a~N(6.125，3.275 541²)，其概率密度函数$p\left( a \right) = \frac{1}{{\widehat {{\sigma _a}}}} \cdot \varphi \left( {\frac{{a-6.125}}{{3.275\;541}}} \right)$。相应地，漂移系数和扩散系数取均值，即${\hat v}$=0.039 766，${\hat \delta }$=0.825 169。a的积分上下限取(6.125±3×3.275 541)，补偿系数k=$\frac{1}{{2\mathit{\Phi }\left( 3 \right)-1}}$=1.002 71, 则可靠度函数为

$ \begin{array}{*{20}{c}} {R\left( t \right) = 1.002\;71\int_{ - 3.701\;623}^{15.951\;623} {\left( {\mathit{\Phi }\left( {\frac{{{D_{\rm{f}}} - \hat a - \hat vt}}{{\hat \delta \sqrt t }}} \right) - } \right.} }\\ {\left. {\exp \left( {\frac{{2\hat v\left( {{D_{\rm{f}}} - \hat a} \right)}}{{{{\hat \delta }^2}}}} \right)\mathit{\Phi }\left( {\frac{{ - \left( {{D_{\rm{f}}} - \hat a} \right) - \hat vt}}{{\hat \delta \sqrt t }}} \right)} \right) \cdot }\\ {\frac{1}{{\widehat {{\sigma _a}}}}\varphi \left( {\frac{{a - 6.125}}{{3.275\;541}}} \right){\rm{d}}a.} \end{array} $

(15)

令失效阈值D_f=2′，可靠度曲线如图 2所示，并将基于伪寿命分布法的可靠度曲线同时绘于图 2中。

根据式（14），概率密度函数为

$ \begin{array}{*{20}{c}} {f\left( t \right) = 1.002\;71\int_{ - 3.701\;623}^{15.951\;623} {\left( {\frac{{{D_{\rm{f}}} - \widehat {{a_0}}}}{{\sqrt {2{\rm{ \mathsf{ π} }}{{\hat \delta }^2}{t^3}} }} \cdot } \right.} }\\ {\left. {\exp \left( { - \frac{{{{\left( {{D_{\rm{f}}} - \widehat {{a_0}} - \hat vt} \right)}^2}}}{{2{{\hat \delta }^2}t}}} \right)} \right) \cdot }\\ {\frac{1}{{\widehat {{\sigma _a}}}}\varphi \left( {\frac{{a - 6.125}}{{3.275\;541}}} \right){\rm{d}}a.} \end{array} $

(16)

概率密度函数曲线见图 3，同时将基于伪寿命分布法的结果绘于图 3中。

从图 2和3可以看出，基于Wiener过程的数控转台可靠性估计结果要高于基于伪寿命分布法的可靠性估计结果。在开始的一段时间内，基于Wiener过程的可靠性估计结果显示产品出现失效的概率非常低；而基于伪失效寿命分布法的结果表明产品出现失效的概率较大。从试验过程以及厂家反馈的实际生产使用情况来看，转台在刚开始使用时通常不会出现失效的情况，因此从实际使用结果的角度而言，基于Wiener过程的可靠性评估结果更符合实际情况。从理论上看，数控转台的性能在使用过程中由于环境、试验条件、使用方法以及测量等因素的影响具有随机波动性。对于退化过程存在随机波动的产品，由于回归模型的方差是常量，无法对退化过程中退化量的波动性进行建模，且本文的例子由于属于极小子样情况，也无法对大样本下不同样本的回归方程系数进行建模来评估产品总体的退化失效随机性。因此，基于伪寿命分布法得出的产品可靠性必然小于实际值。Wiener过程数学模型是由Brown运动和线性漂移项合成的，对产品的性能退化过程波动性具有很好的建模效果。因此，从理论分析的角度而言，基于Wiener过程的可靠性评估结果也更符合实际情况。

为了使上一段的论述更具说服力，增加一组试验数据（见表 4），并得到其可靠度函数曲线（见图 4）和失效概率密度函数曲线（见图 5）。图 4和5所示的结果同样表明：基于Wiener过程的数控转台可靠性估计结果要高于基于伪寿命分布法的可靠性估计结果，与上一段的结论一致。可见，对于数控转台的可靠性分析而言，基于Wiener过程的分析结果比基于伪寿命分布法的分析结果更加符合实际情况。

3 结论

本文根据数控转台性能退化的特点，采用Wiener过程对数控转台的性能退化数据进行可靠性建模。利用极大似然法进行模型参数估计，结合试验数据分析结果将退化模型中的性能初值假设为服从正态分布的随机变量而不是定值，充分考虑了数控转台由于各种外界原因产生的个体差异性问题，提高了可靠性评估的精度。最后，将用本文方法得到的可靠性评估结果与基于伪寿命分布法的评估结果进行比较分析，结果表明基于Wiener退化过程的分析结果更加符合实际情况。