樓洪梁,陳磊,李興林,但召江,陳炳順

(1.中國(guó)計(jì)量學(xué)院 質(zhì)量與安全工程學(xué)院，杭州 310018；2.杭州軸承試驗(yàn)研究中心有限公司 博士后科研工作站，杭州 310022；3.杭州誠(chéng)信汽車軸承有限公司，杭州 310024)

隨著產(chǎn)品質(zhì)量不斷提高，在軸承等機(jī)械產(chǎn)品的可靠性定時(shí)截尾試驗(yàn)中，常會(huì)出現(xiàn)無失效數(shù)據(jù)，特別是在高可靠性、小樣本試驗(yàn)中，滾動(dòng)軸承的壽命可以認(rèn)為近似服從Weibull分布。在壽命試驗(yàn)過程中，如果沒有取得失效數(shù)據(jù)，則使用傳統(tǒng)的統(tǒng)計(jì)方法，如極大似然法、最佳線性不變估計(jì)等，不能對(duì)其可靠性進(jìn)行有效評(píng)價(jià)。所以，對(duì)無失效數(shù)據(jù)情況下軸承的可靠性估計(jì)方法展開研究具有重要的理論意義和實(shí)用價(jià)值。

文獻(xiàn)[1]最先對(duì)無失效數(shù)據(jù)下產(chǎn)品的可靠性展開研究，并取得一定的成果。文獻(xiàn)[2-5]利用修正似然函數(shù)的思想，給出了無失效數(shù)據(jù)下參數(shù)估計(jì)的方法。隨后，Bayes方法逐漸成為處理無失效數(shù)據(jù)的主流方法。文獻(xiàn)[6-8]分別取均勻分布、Beta分布和Gamma先驗(yàn)分布對(duì)ti時(shí)刻的失效概率Pi進(jìn)行Bayes估計(jì)，用加權(quán)最小二乘法配布一條Weibull曲線，從而獲得分布參數(shù)的點(diǎn)估計(jì)，最后獲得可靠性和壽命的估計(jì)。Bayes估計(jì)方法中，選取不同的先驗(yàn)分布會(huì)對(duì)估計(jì)結(jié)果產(chǎn)生嚴(yán)重影響[9]。文獻(xiàn)[10-12]用Beta分布作為失效概率P的先驗(yàn)分布，對(duì)該取值作了解釋，并討論了該取值方法對(duì)估計(jì)結(jié)果的影響。文獻(xiàn)[13]提出了一種綜合估計(jì)方法，在k組無失效數(shù)據(jù)后，再加入一組虛擬的數(shù)據(jù)，并在其中引進(jìn)失效信息對(duì)無失效數(shù)據(jù)下的估計(jì)進(jìn)行修正。文獻(xiàn)[14-15]將文獻(xiàn)[7]中的先驗(yàn)分布由Beta分布改為不完全Beta分布，但沒有給出這種改變的理論依據(jù)，也沒有對(duì)估計(jì)結(jié)果的穩(wěn)定性進(jìn)行討論。

下文在以上研究的基礎(chǔ)上，引入虛擬失效信息，結(jié)合試驗(yàn)信息，通過Bayes方法得出新的估計(jì)結(jié)果。

1 無失效數(shù)據(jù)模型

2 Weibull分布參數(shù)與可靠度的估計(jì)

獲得k組數(shù)據(jù)后，就可通過加權(quán)最小二乘法得出參數(shù)m,η的估計(jì)。令

yi=lnti,

(1)

由此可得

得到Weibull分布參數(shù)后，對(duì)壽命和可靠度作出估計(jì)

(2)

(3)

3 引入虛擬失效信息的方法

3.1 t1時(shí)刻r1的Bayes估計(jì)

對(duì)于可靠度r1，選擇不完全Beta分布作為其先驗(yàn)分布，由于沒有失效數(shù)據(jù)，根據(jù)工程實(shí)踐和專家經(jīng)驗(yàn)給r1取下限μ(0≤μ<1)，取r1的先驗(yàn)分布為(μ,1)的不完全Beta分布為

(4)

式中：a,b為超參數(shù)，a>0,b>0。

由于試驗(yàn)過程中無樣本失效，可知r1大的可能性大，小的可能性小，所以π(r1|a,b)為r1的增函數(shù)。當(dāng)a>1,0

從Bayes估計(jì)的穩(wěn)定性來看，參數(shù)c的取值不宜過大，一般取(2,8)比較適宜[14]。

r1的先驗(yàn)密度函數(shù)π(r1)為

。(5)

由于t1是第1個(gè)截尾時(shí)間點(diǎn)，故可引入r1的估計(jì)無虛擬失效信息，在失效數(shù)服從二項(xiàng)分布的情況下，t1時(shí)刻的似然函數(shù)q(A1|r1) =r1s1(A1為t1時(shí)刻有s1個(gè)樣品未失效事件)。

r1的后驗(yàn)分布為

(6)

其中邊緣分布函數(shù)

在平方損失函數(shù)下r1的Bayes估計(jì)為

(7)

3.2 ti(i>1)時(shí)刻ri的Bayes估計(jì)

第i-1組樣品在ti-1時(shí)刻未出現(xiàn)失效，但在ti時(shí)刻可能出現(xiàn)l個(gè)失效，將這一虛擬失效信息和ti時(shí)刻無失效數(shù)據(jù)綜合起來對(duì)ril做出估計(jì)。

假設(shè)ni-1個(gè)樣本在ti時(shí)刻的失效個(gè)數(shù)服從二項(xiàng)分布，則ni-1個(gè)樣本在ti時(shí)刻出現(xiàn)了l個(gè)失效的概率P(Bil|ril)(Bil為ti-1時(shí)刻的ni-1個(gè)無失效樣本在ti時(shí)刻出現(xiàn)了l個(gè)失效事件)為

，(8)

式中：C,D分別為第i-1組樣品在ti-1時(shí)刻未失效和在ti時(shí)刻未失效的事件；P(C),P(D)分別為第i-1組樣品在ti-1和ti時(shí)刻的可靠度；P(D|C)為第i-1組樣品在ti時(shí)刻沒有失效的概率，顯然P(C|D)=1，因此

(9)

由于在試驗(yàn)中每組試驗(yàn)樣品都是相互獨(dú)立的，即引入的虛擬失效信息Bil和ti時(shí)刻有si個(gè)樣品沒有失效這2個(gè)事件也是相互獨(dú)立的，所以在ti時(shí)刻有si個(gè)樣品沒有失效的概率為P(Ai|ril) =rilsi。引入虛擬失效信息后ril的似然函數(shù)q(AiBil|ril)(Ai為ti時(shí)刻有si個(gè)樣品未失效事件)為

q(AiBil|ril)=P(Bil|ril)P(Ai|ril)=

(10)

。(11)

在Bil虛擬失效信息下，ril的先驗(yàn)分布為

(12)

ril的后驗(yàn)分布為

(13)

其中邊緣分布函數(shù)m(AiBil)為

在平方損失函數(shù)下ril的Bayes估計(jì)為

(14)

由此可得

(15)

4 案例分析

以6006軸承為例，其壽命服從二參數(shù)的Weibull分布，為便于比較，使用文獻(xiàn)[6]中的定時(shí)截尾試驗(yàn)無失效數(shù)據(jù)，見表1。表中，si=∑ni,i=1,2,…，6。

表1 無失效數(shù)據(jù)

由表2可以看出，不同c值下估計(jì)結(jié)果有一定的波動(dòng)。由表3可以看出，在不同的c值下，文中方法得到的特征壽命和形狀參數(shù)估計(jì)值波動(dòng)很小，形狀參數(shù)最大波動(dòng)0.160 4，特征壽命最大波動(dòng)153.6，而方法1與方法2的形狀參數(shù)與特征壽命最大波動(dòng)分別為0.618 4，1 447與0.394 9，334.4。

表2 不同c值下參數(shù)m,η的估計(jì)結(jié)果與可靠度

表3 不同方法得到的參數(shù)m,η的估計(jì)

c=5時(shí)用不同方法得到的Weibull分布失效概率密度圖如圖1所示。由圖可以看出，c=5時(shí)，文中方法所得到的形狀參數(shù)與特征壽命均介于另2種方法之間。綜上所述，文中方法在超參數(shù)c變化時(shí)具有更好的穩(wěn)定性。

圖1 c=5時(shí)不同方法下失效概率密度圖

5 結(jié)束語(yǔ)

提出了一種適合于滾動(dòng)軸承等機(jī)械產(chǎn)品在無失效數(shù)據(jù)條件下，估計(jì)其分布參數(shù)及可靠度與可靠壽命的方法,即在滾動(dòng)軸承截尾試驗(yàn)中出現(xiàn)無失效數(shù)據(jù)時(shí)，在每個(gè)截尾時(shí)間點(diǎn)的可靠度估計(jì)過程中,引入前一個(gè)截尾時(shí)間點(diǎn)的無失效樣本的虛擬失效信息。實(shí)例分析表明，當(dāng)超參數(shù)c變化時(shí)，該方法的估計(jì)結(jié)果波動(dòng)最小，具有更好的穩(wěn)定性與相對(duì)準(zhǔn)確性。Bayes估計(jì)的結(jié)論在很大程度上取決于先驗(yàn)分布的選取，文中提出的估計(jì)方法還需要在生產(chǎn)實(shí)踐中進(jìn)行檢驗(yàn)。