徐 微，胡偉明，孫 鵬

(1.北京天泰志遠(yuǎn)科技有限公司，北京 100029;2.中國(guó)石油塔里木油田塔西南勘探開(kāi)發(fā)公司柯克亞作業(yè)區(qū)輪南項(xiàng)目部，新疆庫(kù)爾勒 841000)

現(xiàn)代機(jī)械設(shè)備正不斷向自動(dòng)化、高速化、復(fù)雜化等方向發(fā)展，一旦其中某部件發(fā)生故障，將直接影響整個(gè)機(jī)械設(shè)備的正常運(yùn)行，所以在故障發(fā)生前進(jìn)行可靠性預(yù)測(cè)，有助于設(shè)備安全性的提高，維修成本的降低，企業(yè)管理水平的提升.

威布爾分布作為一種可靠性模型，已在很多領(lǐng)域內(nèi)得到了不同程度的發(fā)展［1］.但關(guān)于威布爾分布的參數(shù)估計(jì)，盡管前人已有了大量的研究成果，但在工程實(shí)際中仍較難得到應(yīng)用，如溫艷清等人利用Newton-Raphson(牛頓迭代)算法以及完整數(shù)據(jù)的極大化(CM)算法成功實(shí)現(xiàn)了其參數(shù)估計(jì)［2］，但該方法計(jì)算繁瑣且理論性較強(qiáng)，使得工程上應(yīng)用困難，而方華元等人則在極大似然法的基礎(chǔ)上進(jìn)行了改進(jìn)，提出了遺傳算法的優(yōu)化［3］，這種方法雖然提高了精度，但是由于遺傳算法是一種隨機(jī)的選擇，導(dǎo)致了其結(jié)果有一定的波動(dòng)，因此，考慮到這些理論的復(fù)雜性，工程實(shí)施中的難度性，本文提出了比較簡(jiǎn)潔且易于觀察的最小二乘法——水平殘差和垂直殘差平方和最小兩種方法對(duì)模型的參數(shù)估計(jì)，只要在普通坐標(biāo)紙上畫(huà)出經(jīng)變化后的失效時(shí)間的變換x、累積函數(shù)的變換y值或在威布爾概率坐標(biāo)紙上直接畫(huà)出失效時(shí)間t、累積函數(shù)F(t)的原始數(shù)據(jù)，即可完成參數(shù)的估算，從而實(shí)現(xiàn)工程中設(shè)備或部件的可靠性預(yù)測(cè).

1 形狀、尺度參數(shù)威布爾分布模型

威布爾的概率密度函數(shù)為

式中:m為反映不同時(shí)期故障率曲線形狀的參數(shù);n為衡量平均故障間隔時(shí)間的尺度參數(shù);r為壽命極限值或最小值的位置參數(shù)，該值是否取零與模型參數(shù)個(gè)數(shù)有關(guān).

考慮到三參數(shù)威布爾分布模型的復(fù)雜結(jié)構(gòu)，計(jì)算的繁瑣以及在實(shí)際工程中的較少應(yīng)用，因此本文采用簡(jiǎn)化后的兩參數(shù)威布爾分布模型，并利用最小二乘法對(duì)設(shè)備或部件的可靠性進(jìn)行預(yù)測(cè)研究［4］.

二參數(shù)威布爾分布的密度函數(shù)、失效率函數(shù)、累積分布函數(shù)及可靠度函數(shù)分別見(jiàn)圖1—圖4.從圖中可以看出，形狀參數(shù)m值的不同獲得的曲線形狀有很大區(qū)別.

圖1 二參數(shù)威布爾分布的密度函數(shù)Fig.1 Probability density function of two parameters Weibull

圖2 二參數(shù)威布爾分布的失效率函數(shù)Fig.2 Failure rate function of two parameters Weibull

圖3 二參數(shù)威布爾分布的分布函數(shù)Fig.3 Cumulative distribution function of two parameters Weibull

圖4 二參數(shù)威布爾分布的可靠度函數(shù)Fig.4 Reliability function of two parameters Weibull

2 威布爾最小二乘估計(jì)

兩參數(shù)威布爾的累積分布函數(shù)可以表示為

對(duì)等式兩邊進(jìn)行變換后則有

2.1 垂直殘差平方和最小擬合數(shù)據(jù)點(diǎn)的最小二乘法(方法1)

形狀、尺度參數(shù)計(jì)算公式可以表示為

垂直殘差平方和為

2.2 水平殘差平方和最小擬合數(shù)據(jù)點(diǎn)的最小二乘法(方法2)

形狀、尺度參數(shù)計(jì)算公式可以表示為

水平殘差平方和為

結(jié)合式(5)和式(7)可得

綜上，若殘差系數(shù)m＞1時(shí)，則用水平殘差平方和模型進(jìn)行計(jì)算獲得參數(shù)的效果較好，否則，選用垂直殘差平方和模型估算參數(shù)會(huì)更好，同時(shí)說(shuō)明用該組參數(shù)進(jìn)行預(yù)測(cè)效果會(huì)更好.

3 關(guān)于分布函數(shù)的取值

在分布函數(shù)的取值方面，許多學(xué)者常常忽略完整數(shù)據(jù)和截尾數(shù)據(jù)的區(qū)別［7］，而將二者不區(qū)分地加以應(yīng)用，這樣必然會(huì)產(chǎn)生一定的誤差，導(dǎo)致設(shè)備可靠性預(yù)測(cè)的不準(zhǔn)確，管理人員錯(cuò)誤地指導(dǎo)維修管理，造成設(shè)備無(wú)效的維修或維修過(guò)度.

研究發(fā)現(xiàn)，對(duì)威布爾分布函數(shù)取值的選取與要評(píng)估的失效數(shù)據(jù)類型有關(guān)，而試驗(yàn)數(shù)據(jù)和現(xiàn)場(chǎng)數(shù)據(jù)則是失效數(shù)據(jù)的兩種表現(xiàn)類型，若將截尾數(shù)據(jù)按完整數(shù)據(jù)處理或?qū)⑼暾麛?shù)據(jù)按截尾數(shù)據(jù)處理時(shí)均會(huì)產(chǎn)生差異［8］，因此，通過(guò)對(duì)以下分布函數(shù)的分析來(lái)確定各值選取不同對(duì)預(yù)測(cè)結(jié)果的影響.

3.1 完整數(shù)據(jù)的估計(jì)

用平均秩估計(jì)(也稱期望估計(jì))為

式中:F為完整數(shù)據(jù)分布函數(shù)的經(jīng)驗(yàn)取值，F(xiàn)值的具體選擇根據(jù)所需精度而定;i為失效數(shù)據(jù)的序號(hào);k為失效數(shù)據(jù)的個(gè)數(shù).

用中位秩估計(jì)為

對(duì)于F的估計(jì)還有以下兩式:

研究發(fā)現(xiàn)，式(9)的均方誤差最大，式(12)的次之，式(11)的最小，而式(10)卻是最常用的，但當(dāng)k充分大時(shí)，式(9)—(12)給出的F值是非常接近的［9］.

3.2 截尾數(shù)據(jù)的估計(jì)

對(duì)于截尾數(shù)據(jù)的估計(jì)，只要將式(9)—(12)中的i更換為Di=Di－1+(k+1 －Di－1)/(k+2 －j)即可進(jìn)行估算，Di為第i個(gè)失效數(shù)據(jù)的序號(hào)，若令h為截尾數(shù)據(jù)的個(gè)數(shù)，k－h(huán)為失效數(shù)據(jù)的個(gè)數(shù)，以增序排列這組數(shù)據(jù)，則最小數(shù)據(jù)編號(hào)記為1，最大數(shù)據(jù)編號(hào)記為k，最小的失效數(shù)據(jù)編號(hào)記為1，最大的失效數(shù)據(jù)編號(hào)記為k－h(huán)，則j=1，2，…，k，而i=1，2，…，(k－h(huán)).

以某機(jī)械設(shè)備部件的疲勞實(shí)驗(yàn)采集的數(shù)據(jù)(見(jiàn)表1)為例說(shuō)明，其中采用較常用的中位秩作為F的估計(jì)值.

從圖5中可以很清楚地看到，將數(shù)據(jù)歸為完整數(shù)據(jù)和歸為截尾數(shù)據(jù)后所得結(jié)果具有很大不同，因此數(shù)據(jù)判斷的正確性將直接影響到參數(shù)估計(jì)的精度，繼而影響預(yù)測(cè)的準(zhǔn)確度.從最小二乘法關(guān)于垂直殘差最小和水平殘差最小模擬的結(jié)果(見(jiàn)圖6)，可以直觀地看到在此例中二者的差異是存在的，而表2的數(shù)值則更清楚地比較出兩種方法的優(yōu)劣性.

表1 設(shè)備部件的疲勞實(shí)驗(yàn)數(shù)據(jù)Tab.1 Fatigue test data of component of equipment

圖5 完整數(shù)據(jù)和截尾數(shù)據(jù)的比較Fig.5 Comparison between completed data and censored

圖6 垂直殘差和水平殘差平方和最小結(jié)果的比較Fig.6 Comparison of minimizing the sum of squares in horizontal residuals and vertical residuals

表2 兩種方法的參數(shù)值比較Tab.2 Comparison of the parameter value gained from two methods above

綜上可知，方法2的數(shù)據(jù)誤差(Rxy=0.0613)小于方法1的數(shù)據(jù)誤差(Ryx=0.4910)，因此，此例中選擇方法2所得到的參數(shù)值，即n=4328.1和m=2.8713，更能準(zhǔn)確地預(yù)測(cè)該設(shè)備部件的可靠性.

4 可靠性預(yù)測(cè)

將已求得威布爾的兩參數(shù)值分別代入到失效分布函數(shù)公式和可靠度函數(shù)公式中，可以得到［10］:

計(jì)算得到威布爾形狀參數(shù)m=2.8713＞1，則說(shuō)明該部件處于“浴盆曲線”三個(gè)壽命期的耗損階段，若估算5000h后的可靠度則將其代入到式(13)中即可，得到R(5000h)=22.02% ，由此可以預(yù)測(cè)任何時(shí)間點(diǎn)處的可靠性發(fā)展變化情況(如圖7)，根據(jù)本例估算的5000h后所得的可靠性建議采取相應(yīng)措施，如需停機(jī)檢修或更換機(jī)件以提高設(shè)備部件的安全性.

圖7 設(shè)備部件的可靠性預(yù)測(cè)Fig.7 Reliability prediction of component of equipment

5 結(jié)論

本文基于威布爾概率圖的最小二乘法展開(kāi)，在應(yīng)用水平殘差平方和最小和垂直殘差平方和最小分別對(duì)威布爾模型參數(shù)估計(jì)，以確定最佳的參數(shù)值，而在分布函數(shù)值的選取中討論了完整數(shù)據(jù)和截尾數(shù)據(jù)對(duì)預(yù)測(cè)結(jié)果的影響，可以說(shuō)，恰當(dāng)?shù)姆植己瘮?shù)值和最小二乘法的選取能保證可靠性預(yù)測(cè)的準(zhǔn)確性，這對(duì)延長(zhǎng)設(shè)備的壽命，企業(yè)管理人員制定合理的維修決策具有重要的意義.

［1］劉惟信.機(jī)械可靠性設(shè)計(jì)［M］.北京:清華大學(xué)出版社，1996.

LIU Weixin.The mechanical reliability design［M］.Beijing:Tsinghua University Press，1996.

［2］溫艷清，劉寶亮.Weibull分布在完全數(shù)據(jù)條件下的參數(shù)估計(jì)［J］.山西大同大學(xué)學(xué)報(bào)，2009，25(4):17 －19.

WEN Yanqing，LIU Baoliang.The parameter estimation of Weibull distribution under the condition of completed data［J］.Journal of Shanxi Datong University，2009，25(4):17 －19.

［3］方華元，胡昌華，樊紅東，等.基于GA的可靠性壽命分布參數(shù)的極大似然優(yōu)化估計(jì)［J］.上海航天，2006(2):50－53.

FANG Huayuan，HU Changhua，F(xiàn)AN Hongdong，et al.The optimization estimation of maximum likelihood for reliability life distribution based on the GA［J］.Shanghai Aerospace，2006(2):50－53.

［4］于曉紅，張來(lái)斌，王朝暉，等.基于新的威布爾分布參數(shù)估計(jì)法的設(shè)備壽命可靠性分析［J］.機(jī)械強(qiáng)度，2007，29(6):932－936.

YU Xiaohong，ZHANG Laibin，WANG Chaohui，et al.The reliability analysis of equipment life based on the new parameter estimation method of Weibull distribution［J］.Mechanical Strength，2007，29(6):932 －936.

［5］游達(dá)章.最小二乘法在威布爾分布的可靠性評(píng)估［J］.湖北工業(yè)大學(xué)學(xué)報(bào)，2009，24(4):34 －36.

YOU Dazhang.The reliability evaluation about the least square method［J］.Journal of Hubei University of Technology，2009，24(4):34－36.

［6］ZHANG L F，XIE M ，TANG L C.A study of two estimation approaches for parameters of Weibull distribution based on WPP［J］.Reliability Engineering and System Safety，2007，92:360－368.

［7］WU Dongfang， ZHOU Jiancheng， LI Yongdan. Unbiased estimation of Weibull parameters with the linear regression method［J］.Journal of the European Ceramic Society，2006，26:1099–1105.

［8］蔣仁言.威布爾模型族［M］.北京:科學(xué)出版社，1998.

JIANG Renyan. Weibull model-characteristics， parameter estimation and application［M］.Beijing:Science Press，1998.

［9］SAGHAFI A，MIRHABIBI A R，YARI G H.Improved linear regression method forestimatingWeibullparameters［J］.Theoretical and Applied Fracture Mechanics，2009，52:180－182.

［10］黃海，林穗華.基于威布爾分布下一個(gè)失效時(shí)刻的預(yù)測(cè)問(wèn)題［J］.河池學(xué)院月報(bào)，2006，26(5):30 －33.

HUANG Hai，LIN Suihua.The prediction of next failure time based on Weibull distribution［J］.Monthly Report of Stream College，2006，26(5):30－33.