1. 引言
二参数威布尔分布模型广泛应用于机械零部件比如滚动轴承的故障诊断或寿命预测。考虑到产量与成本的因素，如何利用小样本的产品失效试验数据尽可能精确地估计出威布尔参数，成为当前可靠性研究的热点之一[1] 。
常用的威布尔分布参数估计方法包括极大似然估计和线性化最小二乘法，其中后者由于计算简便直观而广泛应用于工程实践中，但估计精度不高。为提高线性化最小二乘法估计威布尔参数的精度，可以变换威布尔参数估计的方向[2] [3] ，这主要是由于单方向的参数估计只考虑了单个方向的坐标误差。另一思路认为在非线性回归模型转化为线性模型时，模型中随机变量的分布特征会改变，不再满足Gauss-Markov假定，从而影响线性最小二乘法估计的精度[4] 。因此，该思路通常对普通最小二乘法进行加权处理，如经典的Bergman加权最小二乘法[5] 。然而，加权最小二乘法需要可靠度的估计值与实际值相差足够小，这样才能保证误差模型的转换成立。因此加权最小二乘法在小样本的情况下并非最佳。
本文从误差模型的角度切入，提出运用泰勒级数展开—最小二乘法(Taylor Series Expansion-Ordinary Least Squares Estimation, TSE-OLS) [6] 来估计威布尔分布参数，利用Monte-Carlo方法考察完全失效数据的拟合效果，并通过实例计算分析TSE-OLS对比普通最小二乘法估计(Ordinary Least Squares Estimation, OLSE)及加权最小二乘法估计(Weighted Least Squares Estimation, WLSE) [5] 的优点。
2. 精度问题的提出
一般地，威布尔分布最小二乘法估计采取反变换法。威布尔分布的累积分布函数有如下形式：
(1)
如果F(t)代表失效概率，则可靠度R(t)可表示为
(2)
对等式两边取两次对数将其线性化为
(3)
令
，
，
，
，从而将(3)式表示成一元线性回归模型的形式：
(4)
其中，e为随机误差。将故障间隔时间或失效时间数据从小到大排序
(n为可靠性试验的样本量)，因此新的回归模型的自变量
也按从小到大排序。
对于无删失数据的试验样本，可以利用近似中位秩公式等方法初步估算出各个失效时间的失效概率
(5)
其中，i = 1 ~ n。根据式(3)可获得因变量的序列
，再由最小二乘原理估计出回归系数
和
：
(6)
最后通过反变换获得形状参数和尺寸参数的估计值：
，
。
上述方法避免了非线性拟合和极大似然估计的复杂数值积分过程，同时解决了原方法可能无法得到最优解的难题，简化了计算，但却使参数估计的精度下降。下面从理论分析的角度，证明通过对数线性化最小后得到的回归系数，不再是最小二乘意义下的最优回归系数。
由最小二乘原理可知[7] ，尽管可以使式(4)中的残差平方和取到最小，如式(7)所示，但无法保证原始回归模型的残差平方和最小。
(7)
由于
(8)
为了使
尽可能小，不妨取理想状态
，那么有
，即所有点都在拟合直线上。通过反变换，有
，代入式(8)则
，这与前述的理想状态是矛盾的。因此，对数线性化后的回归模型不能保证原始的非线性回归模型的残差平方和最小。
3. 对数线性化最小二乘法的加权处理
针对上述问题，Bergman加权处理从原始的误差模型出发得到[8]
(9)
(10)
由式(3)可推导出权函数
，那么根据最小二乘原理，我们推得：
(11)
同样地，通过第1节中的反变换得到形状参数和尺寸参数的估计值。要使加权函数成立，必须满足随机误差e的绝对误差足够小[6] 。但在小样本(n ≤ 20)的情况时，由于近似中位秩误差较大，很难满足该条件，此时加权处理并非符合严格意义上的最小二乘原理。
4. 泰勒级数展开–最小二乘法估计
文[7] 提出利用泰勒级数展开结合最小二乘法的思想来提高对数线性化，他们的工作在幂函数的参数估计中取得了良好的效果。与之不同的是，这里我们主要结合威布尔函数介绍其改进算法。
设
，以OLSE得到的估计参数β0和η0作为初值，将其展开成一阶泰勒级数：
(12)
上式可以写成
(13)
其中：
该方法将非线性回归问题转化为多元线性回归，同时保持了非线性问题的最小二乘原理，而且回归系数没有发生变化，因此为复杂的非线性函数的参数估计带来了便利。这里，利用多元最小二乘法求出形状参数和尺度参数的最佳估计值如下：
(14)
对于任意小的
，若
，
(15)
则
为较理想的估计值。若式(15)不成立，则令
，作为新的初始估计值重复以上步骤，直到满足式(15)，从而保证了参数估计的精度。
5. 数值模拟与实例计算
利用Monte-Carlo法[9] 模拟出若干组不同样本量(每组N = 5000次)的标准威布尔分布(β = η = 1)的随机变量，以均方误差的期望[10] 
来评估三种方法的拟合效果。其模拟结果如图1所示。
从图1可见，WLSE和TSE-OLSE方法的残差平方和都比OLSE小，随着样本量的增加，二者的MSE大小趋于一致。然而在小样本的情况下(n = 5~20)，TSE-OLSE比WLSE下降得更明显，尤其在n = 5时，WLSE的均方误差期望仅比OLSE低4 × 10−4 ，而TSE-OLSE的均方误差期望降低更显著(约1.6 × 10−3)。因此，若可靠性试验限于成本等因素只有少量样本时，TSE-OLSE方法更合适。
接下来通过3个实例数据(如表1~3所示)计算，进一步对比三种参数估计方法。计算结果如表4所
示，其中O代OLSE，W代表WLSE，T代表TSE-OLSE。从表中的均方误差对比可以看出，无论是小样本(n = 7)，中样本(n = 33)，还是大样(n = 105)，TSE-OLSE都优于其它两种方法。
Figure 1. The expectation of mean square errors in different sample sizes
图1. 不同样本量下的均方误差期望
Table 1. Time between failures of complex product A/h
表1. 某复杂产品A的故障间隔时间/h [2]
Table 2. Time between failures of numerical control system B/h
表2. 某数控系统B的故障间隔时间/h [11]
Table 3. Time between failures of complex product C/h
表3. 某复杂产品C的故障间隔时间/h [2]
Table 4. Calculation results of instance A/B/C
表4. 实例A/B/C 计算结果
6. 结论
本文介绍了威布尔分布参数估计的对数线性化最小二乘法及其加权处理，提出了利用泰勒级数展开结合最小二乘法的迭代思想以提高参数估计的精度。通过对标准威布尔分布的完全样本数据进行数值模拟和实例计算，对比分析了泰勒级数展开–最小二乘法、普通最小二乘法、加权最小二乘法等三种方法的拟合效果。结果表明泰勒级数展开–最小二乘法比其它两种方法的均方误差更小，尤其在小样本的情况下更能体现其优势。基于本文工作，可进一步探讨删失样本、不同的形状参数和尺寸参数组合的拟合问题，为工程可靠性试验及零部件寿命预测提供参考。
致谢
感谢国家自然科学基金(11202145, 11202144)和江苏省自然科学基金(BK2012175, BK20130303)对本工作开展提供的支持。