基于MATLAB的威布尔分布参数估计的图形界面设计

唐军军，姜年朝，宋军，徐艳楠，刘达

（总参第六十研究所，江苏 南京 210016）

0 引言

三参数威布尔分布模型能够较好地拟合各类试验数据，在元件失效概率、材料疲劳寿命和强度分布等领域有着极广的应用，是可靠性学科中最为重要和常用的一种概率分布形式之一[1]。由于三参数威布尔概率密度模型中包含尺度参数、位置参数和形状参数3个待定参数，相对其它的概率模型，如正态分布模型，参数的准确估计难度较大。早期做法是用图估法等，但存在着主观性强、小样本误差大等缺陷；研究人员先后提出了最小二乘法、相关系数法[1-2]、概率权重法[3]和参数拟合割线法[4]等，但参数的整个求解过程都比较复杂，往往需要借助计算机编程才能完成。

本文在MATLAB平台上，利用其强大的数值计算功能，设计了一款针对威布尔分布参数估计的图形化交互界面（GUI），该工具集成了几种常用的求解方法，包括最小二乘法、相关系数法及超概率权重法，并可给出相应的故障累积概率、故障率等随时间的变化曲线，还可根据北航傅惠民的理论[5]得到特定置信度下的置信区间，操作简便，计算准确。同时，将该工具应用至某型发动机的MTBF分析评估上，得到了较为可信的计算结果。

1 威布尔分布

三参数威布尔分布对于较复杂的数据，具有较强的适应性。设随机变量t服从三参数威布尔分布，则t的累积分布函数，即寿命分布函数为：

式（1）中：β——形状参数；

η——尺度参数；

t0——位置参数（当t0=0时，即为两参数威布尔分布，η也被称为特征寿命）。

在实际运用中，只知道时间点t，如发动机故障发生的时间点，再由概率统计理论推算故障累积概率F（t），也就是要由一个直接已知量t和一个间接已知量F（t）根据式（1）求解β、η和t0这3个未知量，要求得到准确度较高的解比较困难，由此可以看出，影响求解结果精度的因素有两个：

a）间接已知量F（t）的准确性

对于故障累积概率F（t）的推算，常用的计算公式有近似中位秩公式、海森公式和平均秩（即数学期望）公式等，但在实际过程中，常常会遇到参试样本中途退出的情形，即不完全寿命数据，此时应根据故障样本和中止样本的总顺序号，估计出新顺序号，如式（2），再代入式（1），这样才能得到较为准确的故障累积概率F（t）。

式（2）中：Ak——故障样本的平均顺序号；

i——所有样本的排列顺序号。

b）求解方法本身的优越性

对威布尔分布的3个参数求解的方法较多，在处理不同样本容量的故障数据时，各方法可能差异较大，工程选用时存在一定的困难，往往是选取相关性最好的一组解，使求得的威布尔分布表达式能够更好地拟合原始数据。

2 MATLAB的实现

MATLAB由于其强大的数值计算功能，在数值解析、建模仿真等方面有着重要的应用，利用其GUIDE工具开发图形化交互界面，避免了命令行的繁琐操作，使求解过程简洁、直观、明了。

根据上述对威布尔参数求解过程的分析，本文所设计的GUI计算工具考虑到了故障累积概率F（t）的计算、参数估计方法的选择这两大关键点，给出了多种公式、多种方法供选择，工程应用人员可将自己认为合适的F（t）值直接输入并参与计算，扩大了其适用性，也提高了计算的准确性。整个界面分为3个区域：输入区、选择区和显示区，如图1所示；其内在的逻辑判断图如图2所示。

图1 威布尔分布参数估计GUI

图2 威布尔分布参数估计GUI逻辑判断图

将所需要分析的试验数据（如发动机的故障数据、材料性能的测试数据等）输入其中，通过在计算界面上的简单操作，结果可直观地显示至设计界面上，从而得到能描述该组数据的威布尔分布模型的参数，进而为开展产品可靠性水平、材料疲劳寿命等方面的研究提供帮助。

3 发动机MTBF的计算

本文以某型发动机在某试验期间的故障数据为例，简要说明计算过程。将所收集到的该型发动机的故障，按GJB 451A等相关规定进行分类和统计，剔除诸如人为因素等引起的非相关故障，在0～1000 h的总试验时间内共收集到了21个相关故障数据和4个试验中止数据，因此，可将这25个数据看作是不完全寿命数据，具体的故障数据如表1所示，并给出其故障比率随工作时间的分布柱状图和故障率曲线，如图3、4所示。

表1 某型发动机试验期间责任故障统计表

图3 某型发动机故障随工作时间的分布

图4 某型发动机故障率曲线

3.1 威布尔分布的F检验

确定产品寿命、强度及其他特性参数服从何种分布，是可靠性数据分析中的一项重要的前提工作，往往是利用已知的特性数据作分布拟合检验。本文利用F检验法[6]对上述数据进行威布尔分布假设检验，取原假设H0：t1＜t2＜…＜tn来自威布尔分布，F检验的统计量为：

式中：r=21，r1=[r/2]=10，而

其中，E（Zi）是标准极值分布顺序统计量之均值。

当 Fα/2（2（r-r1-1），2r1）＜F

对表1中的数据进行计算，取显著性水平α=0.1，得到的计算结果如表2所示，由于F0.05（20，20）

表2 F检验计算结果

3.2 MTBF计算

将表1中的数据及通过式（2）计算得到的中位秩输入至自行开发的GUI计算工具中，选择相关系数法对符合该型发动机的威布尔分布参数及其MTBF进行求解，并给出在置信度γ=95%的MTBF单侧置信上下限，计算如图5所示。

图5 某型发动机MTBF计算过程

由图5得知，该型发动机MTBF分布均值为51.58 h，此外，得知置信度γ=95%的MTBF单侧下限为28.64 h，单侧上限为61.99 h。因此，在置信度为2 γ-1=90%的MTBF置信区间为[28.64 h，61.99 h]，即该型发动机的MTBF真值有90%的概率落在此区间内，计算结果如表3所示。

根据以上计算得知的威布尔分布参数，得到该型发动机的故障累积概率等特性曲线，如图6所示，其故障累积概率曲线（图6（a））与实际故障点的吻合度较好，相关系数达0.99146，表明求得的威布尔分布表达式能很好地描述该型发动机的故障发生规律；同时其故障率曲线（图6（d））的形状和走势与图4相吻合，也反映了计算结果的正确性。

表3 某型发动机故障数据威布尔分布模型计算结果

图6 某型发动机威布尔分布特性曲线

4 结束语

本文利用MATLAB强大的数值运算功能，运用其GUIDE工具开发了威布尔分布参数估计的计算工具，将原先复杂、繁琐的求解过程封装至GUI界面的后台程序中，可通过简单、快捷的前台操作得到符合可靠性评估、材料分析等工程实际的计算结果，并集成了多种求解算法，对不同数据的处理具有较好的适应性。

此外，将该GUI计算界面运用至某型发动机25个故障数据的计算分析中，得到了其MTBF分布均值、给定置信度下的置信区间和特性曲线。结果表明，通过该界面求得的结果与故障数据的吻合性较好，相关系数达0.99146，一方面，表明威布尔分布模型能较好地描述该发动机的实际故障分布规律；另一方面，也说明了该界面能够满足一般的实际工程需要。

[1]傅惠民，高镇同.确定威布尔分布三参数的相关系数法[J] .航空学报，1990，11（7）： 323-327.

[2]孔瑞莲.航空发动机可靠性工程[M].北京：北京航空航天大学出版社，1992.

[3]张秀之.概率权重矩法及其在Weibull分布参数估计中的应用海洋预报[J].1994，11（3）：56-61.

[4]郭必柱，邓建.可靠性分析威布尔三参数估计方法比较分析[J].科学技术与工程，2010，10（25）：6117-6122.

[5]傅惠民，高镇同，徐人平.三参数威布尔分布的置信限[J] .航空学报，1992，13（3）： 153-162.

[6]赵宇，杨军，马小兵.可靠性数据分析教程[M].北京：北京航空航天大学出版社，2009.