基于MATLAB的威布爾分布參數(shù)估計(jì)的圖形界面設(shè)計(jì)

唐軍軍，姜年朝，宋軍，徐艷楠，劉達(dá)

（總參第六十研究所，江蘇 南京 210016）

0 引言

三參數(shù)威布爾分布模型能夠較好地?cái)M合各類試驗(yàn)數(shù)據(jù)，在元件失效概率、材料疲勞壽命和強(qiáng)度分布等領(lǐng)域有著極廣的應(yīng)用，是可靠性學(xué)科中最為重要和常用的一種概率分布形式之一[1]。由于三參數(shù)威布爾概率密度模型中包含尺度參數(shù)、位置參數(shù)和形狀參數(shù)3個(gè)待定參數(shù)，相對(duì)其它的概率模型，如正態(tài)分布模型，參數(shù)的準(zhǔn)確估計(jì)難度較大。早期做法是用圖估法等，但存在著主觀性強(qiáng)、小樣本誤差大等缺陷；研究人員先后提出了最小二乘法、相關(guān)系數(shù)法[1-2]、概率權(quán)重法[3]和參數(shù)擬合割線法[4]等，但參數(shù)的整個(gè)求解過程都比較復(fù)雜，往往需要借助計(jì)算機(jī)編程才能完成。

本文在MATLAB平臺(tái)上，利用其強(qiáng)大的數(shù)值計(jì)算功能，設(shè)計(jì)了一款針對(duì)威布爾分布參數(shù)估計(jì)的圖形化交互界面（GUI），該工具集成了幾種常用的求解方法，包括最小二乘法、相關(guān)系數(shù)法及超概率權(quán)重法，并可給出相應(yīng)的故障累積概率、故障率等隨時(shí)間的變化曲線，還可根據(jù)北航傅惠民的理論[5]得到特定置信度下的置信區(qū)間，操作簡(jiǎn)便，計(jì)算準(zhǔn)確。同時(shí)，將該工具應(yīng)用至某型發(fā)動(dòng)機(jī)的MTBF分析評(píng)估上，得到了較為可信的計(jì)算結(jié)果。

1 威布爾分布

三參數(shù)威布爾分布對(duì)于較復(fù)雜的數(shù)據(jù)，具有較強(qiáng)的適應(yīng)性。設(shè)隨機(jī)變量t服從三參數(shù)威布爾分布，則t的累積分布函數(shù)，即壽命分布函數(shù)為：

式（1）中：β——形狀參數(shù)；

η——尺度參數(shù)；

t0——位置參數(shù)（當(dāng)t0=0時(shí)，即為兩參數(shù)威布爾分布，η也被稱為特征壽命）。

在實(shí)際運(yùn)用中，只知道時(shí)間點(diǎn)t，如發(fā)動(dòng)機(jī)故障發(fā)生的時(shí)間點(diǎn)，再由概率統(tǒng)計(jì)理論推算故障累積概率F（t），也就是要由一個(gè)直接已知量t和一個(gè)間接已知量F（t）根據(jù)式（1）求解β、η和t0這3個(gè)未知量，要求得到準(zhǔn)確度較高的解比較困難，由此可以看出，影響求解結(jié)果精度的因素有兩個(gè)：

a）間接已知量F（t）的準(zhǔn)確性

對(duì)于故障累積概率F（t）的推算，常用的計(jì)算公式有近似中位秩公式、海森公式和平均秩（即數(shù)學(xué)期望）公式等，但在實(shí)際過程中，常常會(huì)遇到參試樣本中途退出的情形，即不完全壽命數(shù)據(jù)，此時(shí)應(yīng)根據(jù)故障樣本和中止樣本的總順序號(hào)，估計(jì)出新順序號(hào)，如式（2），再代入式（1），這樣才能得到較為準(zhǔn)確的故障累積概率F（t）。

式（2）中：Ak——故障樣本的平均順序號(hào)；

i——所有樣本的排列順序號(hào)。

b）求解方法本身的優(yōu)越性

對(duì)威布爾分布的3個(gè)參數(shù)求解的方法較多，在處理不同樣本容量的故障數(shù)據(jù)時(shí)，各方法可能差異較大，工程選用時(shí)存在一定的困難，往往是選取相關(guān)性最好的一組解，使求得的威布爾分布表達(dá)式能夠更好地?cái)M合原始數(shù)據(jù)。

2 MATLAB的實(shí)現(xiàn)

MATLAB由于其強(qiáng)大的數(shù)值計(jì)算功能，在數(shù)值解析、建模仿真等方面有著重要的應(yīng)用，利用其GUIDE工具開發(fā)圖形化交互界面，避免了命令行的繁瑣操作，使求解過程簡(jiǎn)潔、直觀、明了。

根據(jù)上述對(duì)威布爾參數(shù)求解過程的分析，本文所設(shè)計(jì)的GUI計(jì)算工具考慮到了故障累積概率F（t）的計(jì)算、參數(shù)估計(jì)方法的選擇這兩大關(guān)鍵點(diǎn)，給出了多種公式、多種方法供選擇，工程應(yīng)用人員可將自己認(rèn)為合適的F（t）值直接輸入并參與計(jì)算，擴(kuò)大了其適用性，也提高了計(jì)算的準(zhǔn)確性。整個(gè)界面分為3個(gè)區(qū)域：輸入?yún)^(qū)、選擇區(qū)和顯示區(qū)，如圖1所示；其內(nèi)在的邏輯判斷圖如圖2所示。

圖1 威布爾分布參數(shù)估計(jì)GUI

圖2 威布爾分布參數(shù)估計(jì)GUI邏輯判斷圖

將所需要分析的試驗(yàn)數(shù)據(jù)（如發(fā)動(dòng)機(jī)的故障數(shù)據(jù)、材料性能的測(cè)試數(shù)據(jù)等）輸入其中，通過在計(jì)算界面上的簡(jiǎn)單操作，結(jié)果可直觀地顯示至設(shè)計(jì)界面上，從而得到能描述該組數(shù)據(jù)的威布爾分布模型的參數(shù)，進(jìn)而為開展產(chǎn)品可靠性水平、材料疲勞壽命等方面的研究提供幫助。

3 發(fā)動(dòng)機(jī)MTBF的計(jì)算

本文以某型發(fā)動(dòng)機(jī)在某試驗(yàn)期間的故障數(shù)據(jù)為例，簡(jiǎn)要說明計(jì)算過程。將所收集到的該型發(fā)動(dòng)機(jī)的故障，按GJB 451A等相關(guān)規(guī)定進(jìn)行分類和統(tǒng)計(jì)，剔除諸如人為因素等引起的非相關(guān)故障，在0～1000 h的總試驗(yàn)時(shí)間內(nèi)共收集到了21個(gè)相關(guān)故障數(shù)據(jù)和4個(gè)試驗(yàn)中止數(shù)據(jù)，因此，可將這25個(gè)數(shù)據(jù)看作是不完全壽命數(shù)據(jù)，具體的故障數(shù)據(jù)如表1所示，并給出其故障比率隨工作時(shí)間的分布柱狀圖和故障率曲線，如圖3、4所示。

表1 某型發(fā)動(dòng)機(jī)試驗(yàn)期間責(zé)任故障統(tǒng)計(jì)表

圖3 某型發(fā)動(dòng)機(jī)故障隨工作時(shí)間的分布

圖4 某型發(fā)動(dòng)機(jī)故障率曲線

3.1 威布爾分布的F檢驗(yàn)

確定產(chǎn)品壽命、強(qiáng)度及其他特性參數(shù)服從何種分布，是可靠性數(shù)據(jù)分析中的一項(xiàng)重要的前提工作，往往是利用已知的特性數(shù)據(jù)作分布擬合檢驗(yàn)。本文利用F檢驗(yàn)法[6]對(duì)上述數(shù)據(jù)進(jìn)行威布爾分布假設(shè)檢驗(yàn)，取原假設(shè)H0：t1＜t2＜…＜tn來自威布爾分布，F(xiàn)檢驗(yàn)的統(tǒng)計(jì)量為：

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

其中，E（Zi）是標(biāo)準(zhǔn)極值分布順序統(tǒng)計(jì)量之均值。

當(dāng) Fα/2（2（r-r1-1），2r1）＜F

對(duì)表1中的數(shù)據(jù)進(jìn)行計(jì)算，取顯著性水平α=0.1，得到的計(jì)算結(jié)果如表2所示，由于F0.05（20，20）

表2 F檢驗(yàn)計(jì)算結(jié)果

3.2 MTBF計(jì)算

將表1中的數(shù)據(jù)及通過式（2）計(jì)算得到的中位秩輸入至自行開發(fā)的GUI計(jì)算工具中，選擇相關(guān)系數(shù)法對(duì)符合該型發(fā)動(dòng)機(jī)的威布爾分布參數(shù)及其MTBF進(jìn)行求解，并給出在置信度γ=95%的MTBF單側(cè)置信上下限，計(jì)算如圖5所示。

圖5 某型發(fā)動(dòng)機(jī)MTBF計(jì)算過程

由圖5得知，該型發(fā)動(dòng)機(jī)MTBF分布均值為51.58 h，此外，得知置信度γ=95%的MTBF單側(cè)下限為28.64 h，單側(cè)上限為61.99 h。因此，在置信度為2 γ-1=90%的MTBF置信區(qū)間為[28.64 h，61.99 h]，即該型發(fā)動(dòng)機(jī)的MTBF真值有90%的概率落在此區(qū)間內(nèi)，計(jì)算結(jié)果如表3所示。

根據(jù)以上計(jì)算得知的威布爾分布參數(shù)，得到該型發(fā)動(dòng)機(jī)的故障累積概率等特性曲線，如圖6所示，其故障累積概率曲線（圖6（a））與實(shí)際故障點(diǎn)的吻合度較好，相關(guān)系數(shù)達(dá)0.99146，表明求得的威布爾分布表達(dá)式能很好地描述該型發(fā)動(dòng)機(jī)的故障發(fā)生規(guī)律；同時(shí)其故障率曲線（圖6（d））的形狀和走勢(shì)與圖4相吻合，也反映了計(jì)算結(jié)果的正確性。

表3 某型發(fā)動(dòng)機(jī)故障數(shù)據(jù)威布爾分布模型計(jì)算結(jié)果

圖6 某型發(fā)動(dòng)機(jī)威布爾分布特性曲線

4 結(jié)束語

本文利用MATLAB強(qiáng)大的數(shù)值運(yùn)算功能，運(yùn)用其GUIDE工具開發(fā)了威布爾分布參數(shù)估計(jì)的計(jì)算工具，將原先復(fù)雜、繁瑣的求解過程封裝至GUI界面的后臺(tái)程序中，可通過簡(jiǎn)單、快捷的前臺(tái)操作得到符合可靠性評(píng)估、材料分析等工程實(shí)際的計(jì)算結(jié)果，并集成了多種求解算法，對(duì)不同數(shù)據(jù)的處理具有較好的適應(yīng)性。

此外，將該GUI計(jì)算界面運(yùn)用至某型發(fā)動(dòng)機(jī)25個(gè)故障數(shù)據(jù)的計(jì)算分析中，得到了其MTBF分布均值、給定置信度下的置信區(qū)間和特性曲線。結(jié)果表明，通過該界面求得的結(jié)果與故障數(shù)據(jù)的吻合性較好，相關(guān)系數(shù)達(dá)0.99146，一方面，表明威布爾分布模型能較好地描述該發(fā)動(dòng)機(jī)的實(shí)際故障分布規(guī)律；另一方面，也說明了該界面能夠滿足一般的實(shí)際工程需要。

[1]傅惠民，高鎮(zhèn)同.確定威布爾分布三參數(shù)的相關(guān)系數(shù)法[J] .航空學(xué)報(bào)，1990，11（7）： 323-327.

[2]孔瑞蓮.航空發(fā)動(dòng)機(jī)可靠性工程[M].北京：北京航空航天大學(xué)出版社，1992.

[3]張秀之.概率權(quán)重矩法及其在Weibull分布參數(shù)估計(jì)中的應(yīng)用海洋預(yù)報(bào)[J].1994，11（3）：56-61.

[4]郭必柱，鄧建.可靠性分析威布爾三參數(shù)估計(jì)方法比較分析[J].科學(xué)技術(shù)與工程，2010，10（25）：6117-6122.

[5]傅惠民，高鎮(zhèn)同，徐人平.三參數(shù)威布爾分布的置信限[J] .航空學(xué)報(bào)，1992，13（3）： 153-162.

[6]趙宇，楊軍，馬小兵.可靠性數(shù)據(jù)分析教程[M].北京：北京航空航天大學(xué)出版社，2009.