基于生存分析與極大似然法的可靠性評估

戴 怡，劉 琴

(1.天津職業(yè)技術師范大學 機械工程學院，天津 300222；2.天津職業(yè)技術師范大學 理學院，天津 300222)

0 引言

國家“高檔數(shù)控機床與基礎制造裝備”科技重大專項2016年提出，我國數(shù)控機床平均無故障工作時間(Mean Time Between Failures，MTBF)要達到2 000 h[1]，然而如何評估MTBF值是一個難題，另外可靠性評估還是撬動市場機制的重要手段。由于評估時間和評估成本的限制，可靠性評估試驗一般都是截尾的，即數(shù)學家們定義的“刪失數(shù)據(jù)”。 嚴格來說，MTBF是數(shù)控機床無故障工作時間的數(shù)學期望，其與數(shù)控機床失效概率密度密切相關，文獻[2-6]指出數(shù)控機床的失效概率密度多為Weibull分布；文獻[7-8]對失效概率密度為指數(shù)分布的截尾情況進行了系統(tǒng)地理論闡述；文獻[9]為基于指數(shù)分布且為定時截尾情況下的國家標準。然而，以上文獻對于復雜分布的截尾情況(一般是非指數(shù)分布)缺少相應的理論闡述，直到出現(xiàn)生存分析理論。生存分析是近三、四十年發(fā)展非常迅速的數(shù)學學科，1986年美國國家科學委員會提出的數(shù)學發(fā)展概況中將生存分析列為六大發(fā)展方向之一。在生存分析方面，乘積限估計(1958年)[10]、弱收斂理論(1974年)[11]、大樣本理論(1980年)[12]、漸進正態(tài)收斂理論(1983年)[13]和均值型泛函估計理論(2002年)[14]等，使生存分析理論體系更加堅實和完善，并顯示出強大的生命力。因此，對于復雜分布下截尾情況的數(shù)控機床可靠性，研究其基于生存分析的可靠性評估方法非常有意義。

極大似然法是一個對統(tǒng)計模型普遍適用的估計方法，似然估計的思想雖然簡單，但是十分深刻。文獻[15]給出了截尾情況下的似然函數(shù)；文獻[16]指出，隨機截尾情形下的Weibull分布參數(shù)在非常大的樣本下具有相合性。由于生存分析方法與極大似然法是目前處理截尾數(shù)據(jù)屈指可數(shù)的幾種有效方法中的兩種，為了對比兩種方法的優(yōu)劣，并選擇適合有替換定時截尾試驗的評估方法，本文將分別進行基于生存分析與極大似然法的可靠性評估，并進行對比研究。

1 截尾試驗模型

截尾壽命試驗有隨機截尾、有替換定時截尾、無替換定時截尾、有替換定數(shù)截尾和無替換定數(shù)截尾5種典型形式，圖1和圖2所示分別為隨機截尾模型和有替換定時截尾模型，圖中：t0為定時截尾時間；字母A，B，C，D，E，F(xiàn)，G表示樣本(如數(shù)控機床)；符號“×”表示失效。圖1中樣本進入試驗的時間是隨機的(一般不同)，樣本發(fā)生失效的時間也是隨機的；圖2試驗過程中，樣本失效后馬上替換新樣本進行試驗，試驗開始時間與結束時間為已知和共同的。

在5種可靠性截尾試驗模型中，有替換定時截尾試驗模型最適合描述現(xiàn)場跟蹤試驗。因為現(xiàn)場可靠性試驗范圍廣、時間長，基層工作人員對完全把握故障(或失效)的判斷情況存在難度，因此定數(shù)截尾試驗風險相對較大，而且考慮到可靠性試驗應盡量與企業(yè)生產(chǎn)合拍，以降低試驗成本，有替換試驗模型更合適，所以有替換定時截尾試驗模型是本文研究的重點。本文研究的數(shù)控機床有替換定時截尾壽命試驗模型如圖3所示。圖中假設有n臺機床進行可靠性試驗，當某臺機床發(fā)生故障時立刻進行維修，且修舊如新；每臺機床的工作時間(如x1,1，x1,2)和維修時間(如d1,1，d1,2)依次發(fā)生；n臺機床中某一臺機床發(fā)生故障的次數(shù)是隨機的，一般情況下不同；第1臺至第n臺機床發(fā)生故障的次數(shù)分別為J1,J2,…,Jn。根據(jù)文獻[17-18]的研究，平均維修時間是一個較小的數(shù)值，簡便起見，本文在有替換定時截尾試驗模型中忽略維修時間d1,1，d1,2，…。

2 隨機截尾試驗模型與有替換定時截尾試驗模型之間的關系

生存分析與極大似然法主要研究隨機刪失模型。如上所述，本文認為有替換定時截尾模型適合現(xiàn)場跟蹤可靠性試驗。為了應用生存分析與極大似然法，需要找出兩個模型之間的關系，即有替換定時截尾模型是否屬于隨機截尾這一大類模型。

圖4所示為n臺機床的有替換定時截尾試驗模型，圖中為了簡化問題忽略了維修時間。因為總截尾時間t0在試驗前設定，所以截尾變量Ci,1=Ci,2=…=Ci,j=Cn,1=t0；第i臺機床的工作時間為Ti,1,…,Ti,j,…,Ti,Ji，其相應的截尾變量為Ci,1，…，Ci,j，…，Ci,Ji；用T1，J1，…，Ti,Ji，…，Tn,Jn分別表示第1，…，i，…，n臺機床最后一個工作時間(壽命)的數(shù)據(jù)，且為截尾數(shù)據(jù)；前面的數(shù)據(jù)，如T1,J1-1,…,Ti,Ji-1,…,Tn,Ji-1等為完全失效數(shù)據(jù)。

第i臺機床壽命變量與截尾變量之間的關系為：

(1)

經(jīng)變換得到：

(2)

當這些數(shù)控機床壽命數(shù)據(jù)T1,T2，…,Tn相互獨立且服從某分布時，其截尾變量C1，C2，…，Cn也相互獨立并服從與Ti分布有關的某分布，Ti與Ci的關系類同于生存分析中定義的Ti與Ci之間的關系(詳見3.1節(jié)和式(3))。因此，有替換定時截尾模型可以認為隸屬于隨機截尾模型之大類，有替換定時截尾數(shù)據(jù)可以應用生存分析或極大似然法計算其平均壽命。

3 生存分析模型及其平均壽命積分

3.1 乘積限估計

設T1,T2,…,Tn是非負獨立同分布且表示壽命的隨機變量，其分布函數(shù)記為F；C1,C2…,Cn為非負獨立同分布且表示刪失的隨機變量，其分布函數(shù)記為G；在隨機刪失模型中，不能完全觀察Ti，僅觀察到

Xi=min(Ti,Ci),δ=I[Ti≤Ci];

i=1,2,…,n。

(3)

式中：I[·]表示示性函數(shù)，δ包含了刪失信息。觀察到的數(shù)據(jù)對為(Xi,δi),i=1,2,…,n，設X(1)

(4)

3.2 平均壽命積分

若數(shù)控機床失效分布服從分布函數(shù)F, 則其生存函數(shù)S(x)=1-F，數(shù)控機床平均壽命

(5)

(6)

(7)

式中Mn小于但趨于。Susarla等[12]證明，在一定條件下，正態(tài)漸進收斂于

由于“小散弱”的行業(yè)基本狀況沒有根本改變，這也衍生了很多問題。李殿平表示，這主要集中在三個方面，即信息化發(fā)展水平相對較低、綜合服務能力較弱、假冒偽劣農(nóng)資依然泛濫。

在乘積限定義中，X(1)

4 基于生存分析與極大似然法的可靠性評估

4.1 數(shù)控機床可靠性模擬實驗及截尾數(shù)據(jù)

因為數(shù)控機床可靠性試驗成本高、耗時長，可靠性試驗數(shù)據(jù)全部來自現(xiàn)場試驗不現(xiàn)實，而文獻[2-6]已給出數(shù)控機床大致的失效分布規(guī)律，所以可以進行數(shù)控機床模擬試驗。另外，模擬試驗還可以得到MTBF真值，有助于分析評估方法的優(yōu)劣。試驗樣本無限增多可以反映評估方法最本質的特性，而這正是以概率統(tǒng)計理論為指導的一類分析方法——蒙特卡洛方法的優(yōu)點。

本文模擬試驗假設數(shù)控機床失效服從Weibull分布，用MATLAB軟件產(chǎn)生服從Weibull分布的n行m列隨機數(shù)(Xi,j)n×m，其中每一個數(shù)表示數(shù)控機床一個完整的無故障工作時間，即一個完全失效數(shù)據(jù)。然后根據(jù)有替換定時截尾運行機制，按流程圖(如圖5)產(chǎn)生有替換定時截尾模擬試驗數(shù)據(jù)[19]。需要說明的是，該流程圖產(chǎn)生的n行失效數(shù)據(jù)，每行最后一個為截尾數(shù)據(jù)。

4.2 基于生存分析與極大似然法的數(shù)控機床MTBF值的區(qū)間估計

假設數(shù)控機床失效的故障數(shù)據(jù)服從尺度參數(shù)α=1 500 h、形狀參數(shù)β=1.2的Weibull分布，即MTBF的真值為1 411 h，參加試驗的機床數(shù)為50臺，進行如下模擬仿真試驗和計算，結果如表1所示：

(1)確定試驗總截止時間t0=1 000，1 500，2 000，2 500，3 000，3 500，4 000，4 500，共計8種截尾方案，其中單位為h。

(2)對每一種截尾方案產(chǎn)生200組服從Weibull(1.2，1 500)的完全失效數(shù)據(jù)，然后按照圖5流程進行相應地截尾，產(chǎn)生有替換定時截尾試驗數(shù)據(jù)，每組含50個截尾數(shù)據(jù)和數(shù)量不等的完全失效數(shù)據(jù)。

(3)應用生存分析平均壽命積分公式，計算每一組數(shù)據(jù)的MTBF值，并得到其點估計。對計算結果進行統(tǒng)計分析，統(tǒng)計出點估計值在區(qū)間(1 211，1 611)上覆蓋MTBF真值的概率p，列于表1第2行，p即為一般可靠性評估含義的置信概率。

(4)應用極大似然法計算每一組數(shù)據(jù)的MTBF值，并得到其點估計。對計算結果進行統(tǒng)計分析，統(tǒng)計出點估計值在區(qū)間(1 211，1 611)上覆蓋MTBF真值的概率p，列于表1第3行，p即為一般可靠性評估含義的置信概率。

(5)將每種截尾方案截尾數(shù)據(jù)所占的比例列于表1的第4行。

由表1可見，隨著總截尾時間t0從1 000 h遞增到4 500 h，基于生存分析評估方法，其置信概率從0%提高到97.5%，截尾比例從50%～75%降低到22%～29%，但在截尾時間為1 000 h～3 000 h時，置信概率比較低，不超過68.5%。

表1 截尾時間為1 500 h～4 500 h下兩種方法覆蓋真值的概率

另外，隨著截尾時間t0從1 000 h遞增到4 500 h，應用極大似然評估方法，其置信概率從43.5%提高到94%，在截尾時間為1 000 h～30 00 h，置信概率從43.5%增加到88%，但在3 000 h后，置信概率增加得不明顯。

5 結束語

有替換定時截尾試驗是產(chǎn)生刪失數(shù)據(jù)(或稱截尾數(shù)據(jù))的重要方式之一。生存分析是近幾十年專門研究刪失數(shù)據(jù)的理論，然而本文的仿真結果卻顯示，在通常的試驗時間區(qū)間，極大似然法的評估結果優(yōu)于生存分析，這樣意外的結果提示需審慎對待生存分析的適用場合。

本文得出如下結論：

(1)當MTBF真值為1 411 h，有替換定時截尾試驗的試驗時間為MTBF真值的1.4～2.5倍時，基于極大似然法的可靠性評估置信概率為82.5%～91.5%，明顯高于生存分析方法的0.5%～79.5%。提示在可靠性評估的工程實踐中，當試驗時間是平均壽命的1.5～2.5倍的情況下，極大似然評估結果優(yōu)于生存分析方法。

(2)當試驗時間為MTBF值1.5～2.5倍時極大似然法優(yōu)于生存分析的結論，是通過仿真方法得到的，尚缺乏系統(tǒng)的理論分析。但由于傳統(tǒng)的統(tǒng)計理論基本上是基于大樣本理論的，新的隸屬于智能方法面向有限樣本的統(tǒng)計學習理論，或可以有所作為。

(3)有替換定時截尾壽命試驗適宜這樣一些可靠性評估：①教科書式試驗難以組織，需要結合工程實際抽象模型;②故障判斷有一定難度，且對評估結果關系重大;③試驗對象的壽命不是特別大，以至于難以評估，尤其是試驗時間為平均壽命1.5～2.5倍的情況。