多臺數控機床的時間截尾可靠性評估

王智明，楊建國，王國強，張根保

(1.上海交通大學機械與動力工程學院，200240上海，wangzm301@sjtu.edu.cn; 2.重慶大學機械工程學院先進制造技術實驗室，400030重慶)

在數控機床的可靠性評估中，平均故障間隔時間(tMTBF)、給定時間的可靠度R(t)和一定可靠度下的工作時間t(R)是評估機床可靠性的重要指標之一.可靠性試驗是可靠性評估的基礎，一般常采用多臺定時截尾試驗［1-2］.

分布模型的選擇是評估的難點和關鍵.常見的分布模型有指數分布和Weibull分布，其選擇的方法是先假設故障間隔時間(tTBF)服從某一分布，隨后用K-S或χ2檢驗進行檢驗［3］，但有時會出現2種分布均被通過的情況［4］;對于可線性化的分布模型，可用最小二乘法求解模型，并用相關系數進行檢驗［5］，但該方法所得模型并非最優;另外，還有其他選擇模型的方法，如:綜合法［3］、誤差面積比［4］、灰關聯［6］、神經網絡［7］等方法.但對于少樣本數據，上述方法所得模型精度一般較低.在機床可靠性評估中，參數估計多局限于參數的點估計，對于參數區間估計的報道則較少.

本文基于AIC、BIC信息準則優選tTBF模型和參數，給出了tTBF模型參數和可靠性指標的點估計，用 Fisher信息矩陣法給出了模型參數及R(tMTBF)、t(0.9)等可靠性指標的區間估計.最后用蒙特卡洛法(MC)提高了區間估計的精度.

1 基于AIC、BIC信息準則的模型和方法選擇

1.1 模型選擇

數控機床時間截尾數據通常為不完全數據，其中包含右刪失數據.設某次試驗投入k臺機床，所得數據為t=(ti，tcj)(i=1，2，…，n;j=1，2，…，k)，其中ti為失效數據，tcj為截尾數據.一般每臺機床有一個截尾數據，且不相等.設截尾數據的概率密度和可靠度分別為f(t;θ)、R(t;θ)，則用極大似然法(MLE)優選模型時，其對數似然函數為

令

則可得參數的估計值^θ.

不同模型對試驗數據信息反映的準確性各異，其中，只有對數據擬合度好、參數少、形式簡單的模型才是最佳模型.AIC信息準則將K-L距離和極大似然方法相結合，并利用了似然估計性質［8-9］，其形式如下:

AIC=-2(maximum log likelihood)+2m. BIC信息準則［10-11］是另一種模型選擇的有效方法之一，其形式如下:

BIC=-2(maximum log likelihood)+m×ln(n).其中m是參數數量，n是觀測數據數.可見，BIC對增加的參數加上更大的懲罰項m×ln(n).因此，BIC更有可能選擇參數少的模型.AIC、BIC作為模型選擇的有效方法之一，其值越小者模型越好.

故障間隔時間分布模型一般為指數分布或Weibull分布，因此，以指數分布、Weibull分布或其他分布作為備選分布模型，用MLE法分別求取最大似然函數值，最后得到AIC、BIC值，該值最小者即為所求最佳模型.

1.2 方法選擇

常見的估計參數的方法有最小二乘法(LSM)和MLE法.其中，LSM根據因變量的不同可分X方向秩回歸(RRX)和Y方向秩回歸(RRY)［12］.

Weibull模型可經對數變換線性化，然后用LSM求解.對于不完整數據可用Johnson法得到不完整數據中失效數據的失效順序號，再按中位秩估計得到相應的可靠度［13］.Weibull模型用MLE法估計參數時，所得極大似然估計方程沒有封閉形式，可用迭代法解之.為了使其快速收斂，初值可取LSM所得估計值.最后用AIC、BIC準則優選最佳參數估計方法.

2 參數及可靠性指標的區間估計

2.1 參數的區間估計

通常在大樣本情況下，極大似然估計量是漸進正態分布的.在小樣本情況下，參數的對數更接近正態分布.所以，根據漸近正態原理有

進而可得參數θ的估計區間值如下:

其中，var(^θ)為估計參數的方差，zα/2為置信度為1-α的正態分布的分位數.因此，若求得參數的方差，即可得其估計區間.參數的方差可用Fisher信息矩陣法(FM)解決［14-15］，即參數方差矩陣為Fisher信息矩陣的逆.為了便于說明，下面以兩參數Weibull分布為例.

設時間截尾Weibull分布的失效數為n，截尾數為k，則由式(1)知其對數似然函數為

其中，β為形狀參數，η為尺度參數，均大于零，所以

由Fisher信息矩陣法知

2.2 可靠性指標的區間估計

在獲得模型的參數估計值后，易得其點估計，平均故障間隔時間(tMTBF)、給定時間的可靠度R(t)和一定可靠度下的工作時間t(R)計算公式如下:

2.2.1 給定時間可靠度R(t)的區間估計

在Weibull分布中，設T=ln t，則T服從極小值分布，其可靠性函數為

μ、σ為該函數參數.設μ=ln η，σ=1/β，u= β(ln t-ln η)，則極小值分布可靠性函數可變換為Weibull分布可靠性函數

所以，若u的置信區間已知，便可得出R(t)的置信區間.u的置信區間為

其中，

至此，由式(7)可得出給定時間可靠度的置信區間估計

2.2.2 一定可靠度下的工作時間t(R)的區間估計

在式(6)中，令

則

所以，v的置信區間為 CBv=［vL，vU］ =，最后得到t的置信區間為

3 案例分析

數據來源于文獻［1］，該文給出了10臺加工中心的不完整數據30個，其中失效數據20個，時間截尾數據10個，屬于小樣本數據.具體為: 176.00，248.00，10.50，472.00，45.00，39.00，209.33， 261.25， 510.00，120.00， 224.00，267.50， 32.00， 50.00， 138.50， 398.00，478.00，353.00，348.00，137.06，332.50*，84.00*， 267.00*， 165.00*， 283.16*，700.00*，562.00*，387.00*，383.00*，1.50*(帶*號者為時間截尾數據).

3.1 模型選擇

以指數分布、Weibull分布作為備選分布模型，用MLE法求解，所得結果如表1所示.

表1 模型選擇結果

從表1可見，2參數Weibull分布的AIC、BIC值均較小，所以選取2參數Weibull分布作為時間截尾TBF可靠性評估模型.

3.2 方法選擇

2參數Weibull分布的可靠性函數如式(5)，線性化后令y=ln［-ln(R(t))］，x=ln(t)，則

R(t)可用如下中位秩公式估計:

對于時間截尾數據，失效數據的失效順序ri可用Johnson法估計，公式如下:

得出ri后代替式(11)中的i即可得截尾數據的可靠度Ri.用LSM方法中的RRY估計參數時，由式(10)知β=1.024 4，η=407.942 6.為了選取最佳參數，還可用RRX和MLE法估計，結果如表2所示.由表2可見，當2參數Weibull分布分別用RRX、RRY、MLE 3種不同方法估計參數時，MLE方法所得AIC、BIC最小，模型精度最高.

表2 方法選擇結果

至此，選擇MLE方法得模型參數、機床的tMTBF及該時刻的可靠度R(tMTBF)以及可靠度為0.9的工作時間t(0.9)的點估計分別如下:β= 1.200 9，η =376.989 3，tMTBF=354.545 9，R(tMTBF)=0.395 0，t(0.90)=57.88.

3.3 參數及可靠性指標的區間估計

由式(3)、(4)知參數的協方差矩陣為

由式(2)得參數β、η的置信度為0.95的置信區間分別為［0.827 7，1.742 3］、［261.515 2，543.451 7］.另外，由式(8)得u的方差后知其估計區間為CBu=［-0.511 5，0.365 1］，其中方差var(u)=0.050 0.同樣，由式(9)得v的方差后可得其估計區間，最后得出R(tMTBF)和t(0.9)的估計區間分別為CBR(t)=［0.236 8，0.549 0］，CBt(0.9)=［26.91，124.47］.

由估計區間結果來看，該區間太大，估計精度較低，所以有必要提高區間估計的精度.事實上，上述區間估計采用漸進正態分布原理，適合于大樣本.對于小樣本數據，其估計精度一般較低.通常將樣本單位數不少于30個的樣本稱為大樣本，不及30個的稱為小樣本.數控機床的故障數據一般為少樣本數據，因此，對于少樣本數據可用模擬抽樣的方法來擴大樣本容量，從而提高區間估計精度.

3.4 MC仿真抽樣

2參數Weibull分布的tTBF的抽樣公式為

r為(0，1)區間上均勻分布的隨機變量.取樣本容量 N分別為 500、1 000、2 000、5 000、8 000、10 000，其參數和可靠性指標抽樣區間估計分別如圖1、2所示.從圖可見，樣本容量為8 000和10 000時的估計結果已經很接近，因此停止抽樣，并取容量為10 000時的抽樣結果為最終估計值.其結果如表3所示.從表3可以看出，估計區間均已大大減少.其中，t(0.9)、R(tMTBF)區間長度分別減少93.21%、95.07%.

圖1 參數β、n的區間仿真

圖2 可靠性指標R(tMTBF)、t(0.9)的區間仿真

表3 樣本容量為10 000時的區間模擬抽樣結果

4 討論

1)在Weibull分布模型中，當形狀參數β=1時，Weibull分布就轉換為指數分布，模型處于浴盆曲線的偶然失效區;當β＜1或β＞1時，則模型對應于浴盆曲線的早期失效或耗損失效.可見本案例中的加工中心處于耗損失效期.因此，Weibull分布模型因其靈活性作為故障間隔時間的分布模型而廣泛應用于NC機床的可靠性分析中.在2-Weibull模型的參數估計中，MLE法較之LSE法有更高的估計精度.這是因為LSE法在對數變形代換過程中改變了因變量的形態，變形后因變量殘差平方和最小并不能保證變形前因變量殘差平方和最小，同時，用Johnson法對不完整數據中的失效數據排序時，只考慮了其失效順序而沒有考慮截尾數據的相對大小.換言之，對于不同大小的截尾數據，當用Johnson法分析時，失效數據可能有相同的失效順序.當截尾數據較多時，只擬合重新排序后的失效數據(本案例只有20個失效數據)的LSE法的評估精度自然會受到影響.MLE適用于完整數據和不完整數據，精度較之LSE法有較好的評估效果.但LSE法的評估結果可作為MLE法中迭代計算模型參數的初值，因此，使用MLE法時可結合LSE法而達到快速收斂的目的.本文所提方法雖然是針對時間截尾的，但仍適用于定數截尾實驗.當采用定數截尾時，各有關公式中tcj=0.

5 結論

1)在NC機床的可靠性評估中，對不完整數據可結合LSE法的結果，優先選用MLE法.

2)AIC、BIC信息準則結合了MLE的計算結果和模型的簡約原則，是機床可靠性評估中優選模型和模型參數的有效方法.

3)評估模型參數和可靠性指標的區間時，對于大樣本數據，可采用基于漸進正態分布的FM法，但當樣本數據較少時，可用MC仿真法擴充樣本容量，提高區間估計精度.

［1］DAI Yi，ZHOU Y F，JIA Y Z.Distribution of time between failures of machining center based on type I censored data［J］.Reliability Engineering and System Safety，2003，79:375-377.

［2］張海波，賈亞洲，周廣文.數控系統故障間隔時間分布模型的研究［J］.哈爾濱工業大學學報，2005，37 (2):198-200.

［3］WANG Yiqiang，YAM R C M，ZUO Mingjian.A multicriterion evaluation approach to selection of the best statistical distribution［J］.Computers＆Industrial Engineering，2004，47:165-180.

［4］張英芝，申桂香，賈亞洲，等.數控車床故障分布規律及可靠性［J］.農業機械學報，2006，37(1):156-159.

［5］ZHANG Haibo，JIA Yazhou，ZHOU Guangwen.Time between failures model and failure analysis of CNC system［J］.Journal of Harbin Institute of Technology (New Series)，2007，14(2):197-201.

［6］張英芝，賈亞洲，申桂香，等.數控機床故障分布的灰關聯分析［J］.農業機械學報，2004，35(6):195-197.

［7］王金武，劉家福，許仲祥.履帶式拖拉機可靠性與維修性的分析［J］.農業機械學報，2004，35(4):81-83.

［8］AKAIKE H.Information theory and an extension of the maximum likelihood principle［C］//Proc Second Int Symp on Information Theory.Budapest:Akademiai Kiado，1973:267-281.

［9］SAYAMA S，SEKINE M.Weibull distribution and K-distribution of sea clutter observed by X-band radar and analyzed by AIC［J］.IEICE Transactions on Communications，2000，E83B(9):1978-1982.

［10］SCHWARZ G.Estimating the dimension of a model［J］.Ann Statist，1978(6):461-464.

［11］ADRIAN E，RAFTERY.Bayes factors and BIC:comment on"a critique of the bayesian information criterion for model selection"［J］.Sociological Methods Research，1999，27(3):411-427.

［12］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.

［13］JOHNSON L G.Theory and technique of variation research［R］.New York:Elsevier Publishing Co，1964.

［14］SANGUN P，BALAKRISHNAN N，GANG Zheng. Fisher information in hybrid censored data［J］.Statistics and Probability Letters，2008，78:2781-2786.

［15］ESCOBAR L A，MEEKER W Q.The asymptotic equivalence of the fisher information matrices for type I and type II censored data from location-scale families［J］. Communications in statistics theory and methods，2001，30(10):2211-2225.