但召江，樓洪梁，李興林，郭明月，陳炳順

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

隨著軸承質(zhì)量的提高，在可靠性定時(shí)截尾試驗(yàn)過(guò)程中會(huì)出現(xiàn)大量的無(wú)失效數(shù)據(jù)。由于缺乏足夠的試驗(yàn)信息，傳統(tǒng)的估計(jì)方法無(wú)法對(duì)其做出適當(dāng)?shù)目煽啃怨烙?jì)。因此有必要對(duì)無(wú)失效數(shù)據(jù)下軸承的可靠性估計(jì)方法展開(kāi)研究。軸承壽命服從兩參數(shù)Weibull分布，兩分布參數(shù)估計(jì)的準(zhǔn)確性決定軸承可靠性估計(jì)的準(zhǔn)確性。在無(wú)失效數(shù)據(jù)這種乏信息的試驗(yàn)過(guò)程中，Bayes估計(jì)方法能充分利用驗(yàn)前信息和試驗(yàn)信息，顯示了其獨(dú)特的優(yōu)越性。無(wú)失效數(shù)據(jù)的產(chǎn)品可靠性研究起源于文獻(xiàn)[1]，自該文獻(xiàn)發(fā)表以來(lái)，對(duì)無(wú)失效數(shù)據(jù)的研究逐漸受到重視，并取得一定的成果。

目前，在無(wú)失效數(shù)據(jù)情況下，Weibull分布參數(shù)的估計(jì)中，由于參數(shù)的先驗(yàn)信息難以獲取，大多數(shù)都從失效概率入手。在無(wú)失效情況下，失效概率大的可能性小，小的可能性大，根據(jù)這一特性利用共軛分布方法確定失效概率的先驗(yàn)分布，然后通過(guò)Bayes估計(jì)方法得出一組不同時(shí)刻的失效概率估計(jì)值，再通過(guò)最小二乘法求得分布的2個(gè)參數(shù)[2-3]。事實(shí)上，形狀參數(shù)是Weibull分布中極其重要的參數(shù)，其數(shù)值決定了Weibull分布曲線的形狀，如果能獲得形狀參數(shù)的先驗(yàn)信息，對(duì)提高參數(shù)估計(jì)的準(zhǔn)確度具有較大意義。文獻(xiàn)[4]用數(shù)理統(tǒng)計(jì)方法對(duì)軸承形狀參數(shù)進(jìn)行了深入研究，提出了軸承形狀參數(shù)的具體分布，為軸承可靠性估計(jì)帶來(lái)方便。

下文將Weibull分布轉(zhuǎn)化為指數(shù)分布，以指數(shù)分布的失效率為切入點(diǎn)，在無(wú)失效數(shù)據(jù)下，利用共軛分布方法確定失效率的先驗(yàn)分布。Weibull分布形狀參數(shù)的先驗(yàn)信息用2種方法取得：(1)通過(guò)文獻(xiàn)[4]的研究成果，擬合一分布作為形狀參數(shù)的先驗(yàn)分布；(2)根據(jù)生產(chǎn)經(jīng)驗(yàn)以某一區(qū)間上均勻分布作為先驗(yàn)分布，然后利用無(wú)失效試驗(yàn)數(shù)據(jù)，得出失效率和形狀參數(shù)的Bayes估計(jì)，進(jìn)而計(jì)算得到軸承壽命的估計(jì)值，并在先驗(yàn)信息和截尾時(shí)間改變的情況下討論其穩(wěn)定性。

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

軸承壽命T服從形狀參數(shù)為m,特征壽命參數(shù)為η的Weibull分布，其分布函數(shù)F(t)為

(1)

其中，m,η未知?，F(xiàn)從中隨機(jī)抽取n個(gè)樣品分成k組，各組樣品數(shù)分別為n1,n2,n3,…，nk，且n1+n2+n3+…+nk=n，各組樣品分別獨(dú)立地進(jìn)行定時(shí)截尾試驗(yàn)，截尾時(shí)間分別為t1,t2,t3,…，tk，結(jié)果無(wú)一失效，則獲得一組無(wú)失效數(shù)據(jù)(ti,ni),i=1,2,3,…，k。

若令λ=η-m,t*=tm，則Weibull分布形狀參數(shù)和特征壽命的估計(jì)可轉(zhuǎn)化為討論指數(shù)分布的失效率估計(jì)[5]。

設(shè)產(chǎn)品壽命T*服從參數(shù)為λ的指數(shù)分布，其分布函數(shù)為

F(t*)=1-exp(-λt*),λ>0,t*>0。

(2)

影響B(tài)ayes估計(jì)準(zhǔn)確性的關(guān)鍵因素是先驗(yàn)分布。將Weibull分布轉(zhuǎn)化為指數(shù)分布有利于先驗(yàn)信息的獲取。文獻(xiàn)[5]提出先驗(yàn)分布應(yīng)取共軛分布(Gamma分布)。

2 先驗(yàn)信息的獲取

針對(duì)以上模型，在無(wú)失效數(shù)據(jù)下指數(shù)壽命型產(chǎn)品的失效率λ小的可能性大，大的可能性小，因此應(yīng)選取減函數(shù)作為失效率λ的先驗(yàn)分布，由于模型中關(guān)于λ的信息極少，故采用共軛的方法確定先驗(yàn)分布。

在給定形狀參數(shù)m情況下取Ga(a,b)為λ的先驗(yàn)分布密度，其中a,b為分布的超參數(shù)，由文獻(xiàn)[6]可知，當(dāng)00時(shí)，概率密度遞減。則在給定m條件下λ的先驗(yàn)分布密度函數(shù)為

(3)

其中,a,b由2種方法確定[7]：(1) 由專家經(jīng)驗(yàn)確定；(2) 用2個(gè)分位數(shù)確定，分位數(shù)可以根據(jù)先驗(yàn)信息和歷史資料確定。例如用上下四分位數(shù)λU和λL，則a,b須滿足

對(duì)于形狀參數(shù)m的先驗(yàn)信息，可以從產(chǎn)品的歷史信息中獲得，例如對(duì)于某一型號(hào)的軸承，可以從歷史試驗(yàn)數(shù)據(jù)中取得m的先驗(yàn)分布函數(shù)，至少根據(jù)生產(chǎn)方的經(jīng)驗(yàn)可以確定形狀參數(shù)的分布范圍，最壞的情況下取形狀參數(shù)m的先驗(yàn)分布為某一區(qū)間的均勻分布。關(guān)于形狀參數(shù)的分布，文獻(xiàn)[4]分別以滾子軸承、球軸承和軸承總體為母體討論了其分布范圍和平均值，并繪出了柱狀圖，這為以后軸承可靠性的評(píng)價(jià)提供了便利。在此根據(jù)文獻(xiàn)[4]提供的數(shù)據(jù)，以球軸承為母體，將擬合得到的形狀參數(shù)的概率分布作為先驗(yàn)信息。

對(duì)于文獻(xiàn)[4]中的數(shù)據(jù)，通過(guò)擬合試驗(yàn)，在95%的置信水平下，形狀參數(shù)m的先驗(yàn)分布符合Weibull分布。用Matlab提供的方法擬合得到形狀參數(shù)m的分布密度函數(shù)π(m)為

(4)

式中：d=1.743 9；c=2.176 0。

對(duì)于某一型號(hào)的軸承，通常無(wú)法確定形狀參數(shù)的具體值，但可以根據(jù)以往的生產(chǎn)經(jīng)驗(yàn)，確定形狀參數(shù)的區(qū)間，因此，也可以用該區(qū)間內(nèi)的均勻分布作為形狀參數(shù)的先驗(yàn)分布。在得到所有的先驗(yàn)信息的基礎(chǔ)上，就可以利用試驗(yàn)數(shù)據(jù)對(duì)參數(shù)λ,m進(jìn)行Bayes估計(jì)。

3 參數(shù)λ,m的Bayes估計(jì)

由(3)和(4)式得λ,m的先驗(yàn)聯(lián)合密度分布π(λ,m)為

π(λ,m)=π(λ|m)π(m)=

(5)

在無(wú)失效數(shù)據(jù)下，相應(yīng)的似然函數(shù)L(0|λ,m)為[8-9]

(6)

根據(jù)(5)和(6)式得λ,m的后驗(yàn)聯(lián)合密度函數(shù)h(λ,m|0)為

(7)

由(7)式得到m,λ的后驗(yàn)邊緣分布h(m|0)，h(λ|0)為

(8)

(9)

(10)

(11)

(12)

有了Weibull分布的參數(shù)后，就可以對(duì)壽命和可靠度作出估計(jì)，即

(13)

(14)

對(duì)形狀參數(shù)m服從區(qū)間為[m1,m2]均勻分布的情況，只需將(4)式改為

(15)

將(15)式代替(4)式，進(jìn)行(5)～(14)式計(jì)算，即可得到形狀參數(shù)服從均勻分布情況下軸承壽命的可靠度估計(jì)。

4 實(shí)例分析

利用文獻(xiàn)[10]中6203軸承定時(shí)截尾試驗(yàn)結(jié)果(表1)對(duì)上述方法進(jìn)行驗(yàn)證，得到的一組無(wú)失效數(shù)據(jù)見(jiàn)表2，表2的截尾時(shí)間靠近失效時(shí)間。

文獻(xiàn)[10]用最佳線性不變估計(jì)(BLIE)、最佳線性無(wú)偏估計(jì)(BLUE)和最大似然函數(shù)(MLE)3種方法得到估計(jì)結(jié)果，見(jiàn)表3。根據(jù)方法1和表2的數(shù)據(jù)，用(4)式作為形狀參數(shù)的先驗(yàn)分布，先驗(yàn)信息中超參數(shù)a=0.01,b=1.5時(shí)，得到的估計(jì)結(jié)果見(jiàn)表3；用方法2，在先驗(yàn)信息中超參數(shù)a=0.05,b=5, 形狀參數(shù)先驗(yàn)分布取[0.5,4]區(qū)間的均勻分布，得到的估計(jì)結(jié)果見(jiàn)表3。

表1 6203軸承定時(shí)截尾試驗(yàn)數(shù)據(jù)

表2 無(wú)失效數(shù)據(jù)的截尾時(shí)間

表3 文獻(xiàn)[10]的估計(jì)結(jié)果和文中估計(jì)結(jié)果的對(duì)照

當(dāng)先驗(yàn)信息不變，截尾時(shí)間分別為表2的0.4，0.6和0.8倍時(shí)，2種方法的估計(jì)值見(jiàn)表4。

表4 截尾時(shí)間按比例變化的Bayes估計(jì)結(jié)果

由表4可以看出，隨著截尾時(shí)間延長(zhǎng)，形狀參數(shù)的變化不明顯，但特征壽命在不斷增加。對(duì)于方法2，當(dāng)截尾時(shí)間取表2中的0.5倍，形狀參數(shù)先驗(yàn)信息的均勻分布區(qū)間伸縮與平移變化時(shí)，估計(jì)結(jié)果的變化見(jiàn)表5。

由表5可以看出，當(dāng)形狀參數(shù)的均勻分布區(qū)間伸縮變化時(shí)，形狀參數(shù)變化不大，特征壽命的變化較大。當(dāng)形狀參數(shù)的均勻分布區(qū)間平移變化時(shí)，特征壽命變化不大，形狀參數(shù)的變化較大。這說(shuō)明先驗(yàn)信息的改變對(duì)估計(jì)結(jié)果的影響比截尾時(shí)間改變帶來(lái)的影響更大些。

表5 形狀參數(shù)的均勻分布區(qū)間變化時(shí)估計(jì)結(jié)果

5 結(jié)論

(1)在無(wú)失效數(shù)據(jù)Bayes估計(jì)方法中，先驗(yàn)信息的選取對(duì)結(jié)果的影響很大。在先驗(yàn)信息準(zhǔn)確的情況下，可以得到相對(duì)準(zhǔn)確的估計(jì)結(jié)果。

(2)在Bayes估計(jì)中，先驗(yàn)信息的改變對(duì)估計(jì)結(jié)果的影響更明顯，截尾時(shí)間的影響相對(duì)較小。

(3)從理論上說(shuō)，方法1更可靠些，但需要較多的形狀參數(shù)歷史數(shù)據(jù)；方法2則相對(duì)簡(jiǎn)單易行。因此，具體使用哪種方法需要根據(jù)實(shí)際情況而定。

此外，文中提出的方法能否在軸承可靠性估計(jì)中推廣使用，還需要更進(jìn)一步的實(shí)踐檢驗(yàn)，畢竟Bayes法先驗(yàn)信息的取得在很大程度上是主觀的。