無失效威布爾情形下基于雙修正多層Bayes的可靠性評估

龍足騰， 鄭 波,2， 甯 洋， 羅金超

（1. 中國民用航空飛行學(xué)院 航空電子電氣學(xué)院， 四川 廣漢 618307；2. 核工業(yè)西南物理研究院， 四川 成都 610225）

隨著工業(yè)化、現(xiàn)代化的發(fā)展，越來越多的產(chǎn)品表現(xiàn)出高質(zhì)量、高可靠性、長壽命及短期內(nèi)無失效的特性。在這類產(chǎn)品的可靠性試驗中，經(jīng)常出現(xiàn)無失效的情形，很難通過定時或定數(shù)截尾試驗獲取其完全失效數(shù)據(jù)來進行可靠性分析，這不僅增加了產(chǎn)品研發(fā)的時間和成本，還會嚴重阻礙產(chǎn)品的后續(xù)定型及應(yīng)用。因此，產(chǎn)品無失效情形下的可靠性評估研究具有很強的實際意義，是當(dāng)下可靠性領(lǐng)域的研究重點。

Martz等[1]是最早開始研究無失效數(shù)據(jù)處理的學(xué)者，提出了一種以失效概率為指標，結(jié)合Bayes理論的可靠性評估方法。隨后，茹詩松等[2]提出了配分布曲線法，通過估計無失效數(shù)據(jù)各個截尾時刻的失效概率，找到失效概率與壽命分布函數(shù)之間的關(guān)系，用最小二乘法配一條分布曲線，得到可靠度函數(shù)。由于配分布曲線法流程清晰，計算相對簡便，估計效果穩(wěn)定，在無失效數(shù)據(jù)的可靠性研究和實際工程中得到廣泛應(yīng)用。

配分布曲線法的關(guān)鍵是估計失效概率，相關(guān)學(xué)者就失效概率的估計進行了深入的研究，一些成熟的Bayes方法被相繼提了出來，如多層Bayes[3]和EBayes 法[4-5]。這些方法在試驗中得到驗證，證明了能夠用于無失效情形的分析[6-7]。隨后，它們在處理不同壽命分布類型的無失效數(shù)據(jù)中成為主流方法，尤其在威布爾分布情形的研究中得到廣泛的認可[8-11]。一些學(xué)者和研究人員通過結(jié)合各種估計原理和算法來構(gòu)造可靠性評估模型，從而解決實際工程問題[12-14]。在實際應(yīng)用中，超參數(shù)和試驗分組等影響可靠性評估精度的問題被引入研究，推動了主流Bayes方法在工程問題中的應(yīng)用[15-16]。

隨著研究的深入，一些學(xué)者發(fā)現(xiàn)基于多層Bayes和E-Bayes的估計效果較于保守，于是通過引入Bootstrap法的區(qū)間估計或融合多源信息數(shù)據(jù)來提高模型估計的精度和可信度。李海洋等[17-18]利用EBayes 法完成點估計，通過引入?yún)?shù)Bootstrap 法完成區(qū)間估計并構(gòu)建融合模型用于機器人軸承可靠性評估，并且在對扭轉(zhuǎn)電機實際運行后的無失效數(shù)據(jù)的分析中取得了良好的估計效果；王瑞祥等[19]利用多層Bayes完成點估計，再結(jié)合Bootstrap構(gòu)建小樣本無失效的可靠性評估模型；Chen等[20]針對航空發(fā)動機壓氣機盤的小樣本無失效數(shù)據(jù)，在Bayes估計的基礎(chǔ)上采用了融合測試數(shù)據(jù)與模擬數(shù)據(jù)的評估思路。極少數(shù)學(xué)者旨在改進Bayes 估計，如：李爽等[21]在威布爾分布失效概率估計時，利用后驗分布調(diào)整先驗分布來修正失效概率的取值上界，提出了一種基于多層Bayes的改進方法，將它用于提高失效概率估計精度。

綜上所述，利用配分布曲線處理無失效情形下可靠度點估計時，E-Bayes和多層Bayes是估計失效概率的主流方法，但其估計效果相對保守。為了提高失效概率的估計精度，本文通過改進多層Bayes，提出雙修正多層Bayes，并結(jié)合加權(quán)最小二乘法和可靠度函數(shù)完成可靠度點估計；通過仿真試驗和工程實例分析，來驗證新模型的可行性和有效性；最后，討論了雙修正多層Bayes的適用性。

1 無失效數(shù)據(jù)與威布爾分布

1.1 無失效數(shù)據(jù)結(jié)構(gòu)

在工業(yè)產(chǎn)品中，隨機抽取n個樣本，將它分成k組，進行定時截尾試驗。各組樣本數(shù)分為n1,n2, …,nk，截尾時間為t1,t2, …,tk(t1＜t2＜…＜tk)。若在試驗中無一樣本失效，則稱該樣本組數(shù)據(jù)為無失效數(shù)據(jù)，有n=n1+n2+…nk(i=1, 2, …,k)。無失效數(shù)據(jù)的結(jié)構(gòu)如表1所示。表中未失效數(shù)si（si=nk+nk-1+…+ni） 即試驗中截尾時間ti時刻的未失效數(shù)。

表1 無失效數(shù)據(jù)結(jié)構(gòu)Table 1 Zero-failure data structure

1.2 威布爾分布

許多工業(yè)產(chǎn)品如機械零部件、電子元器件的壽命能夠擬合成威布爾分布。威布爾分布的形狀參數(shù)能很好地反映產(chǎn)品失效特性，可以根據(jù)其值大小判定產(chǎn)品的失效類型，因此威布爾分布對各類型的試驗數(shù)據(jù)具有較強的適應(yīng)能力，在可靠性研究領(lǐng)域有著廣泛應(yīng)用。

本文選取兩參數(shù)參威布爾分布作為研究對象。其失效概率密度函數(shù)為：

失效分布函數(shù)為：

失效率函數(shù)為：

可靠度函數(shù)為：

式中：β為形狀參數(shù)；η為尺度參數(shù)或特征壽命（即可靠度R=0.367 9時的壽命值），它決定了曲線尺寸比例的大小，能夠縮放坐標尺度，并且受負載影響，負載越大，η越小。

形狀參數(shù)β對威布爾分布函數(shù)的影響如圖1所示。

圖1 形狀參數(shù)β對威布爾函數(shù)的影響Fig.1 Ⅰnfluence of shape parameters on Weibull functions

由圖1 可知，當(dāng)β=3.0～5.0 時，威布爾分布函數(shù)圖像與正態(tài)分布圖像很相似。

2 無失效威布爾可靠度估計

2.1 失效概率pi的多層Bayes估計

無失效數(shù)據(jù)在某一時刻ti處的失效概率pi值較大的概率很小，pi值較小的概率很大，因此取pi的先驗分布為共軛Beta 分布，其概率密度函數(shù)為：

式中：0＜pi＜1；B(a,b)為Beta函數(shù)；a和b為超參數(shù)且相互獨立，一般為了保證其密度函數(shù)為減函數(shù)以及Bayes 估計的穩(wěn)健性；設(shè)0＜a≤1，1＜b＜c，c為常數(shù)，根據(jù)文獻[16]，c通常為整數(shù)，在2～7 中選取。

超參數(shù)a和b的先驗分布為：

由式(5)、式(6)可得pi的多層先驗分布為：

由于所得數(shù)據(jù)是無失效數(shù)據(jù)，pi的似然函數(shù)為：

根據(jù)Bayes定理，由式(7)、式(8)可得pi的多層后驗分布密度函數(shù)為：

在平方損失下，取多層后驗分布的期望作為pi的多層Bayes估計，可以表示為：

為了簡化計算，根據(jù)工程建議[19]，取a=1，則計算得，從而得到pi多層Bayes的解析式為：

2.2 失效概率pi的雙修正多層Bayes估計

由于在無失效場合下，pi始終滿足0＜pi＜1，無法取到0和1，因此考慮通過修正pi的取值范圍來減小估計誤差。根據(jù)Bayes 修正理論[21]，后驗分布是對先驗分布的修正與調(diào)整，將其用于估計pi的基本思路是：利用前一時刻pi的后驗分布替代后一時刻pi的先驗分布。

設(shè)pl、pu分別為pi的下限值和上限值，則0＜pl＜pi＜pu＜1，此時Beta分布就不能很好地描述pi的先驗分布。考慮選用不完全Beta 分布，其相較于Beta分布圖像更加集中，這就意味著pi的取值區(qū)間變小，誤差波動也變小。不完全Beta函數(shù)的表達式為：，式中：x1和x2分別對應(yīng)pl和pu。在多層Bayes 的基礎(chǔ)上，通過改變先驗分布，可以得到雙修正多層Bayes的表達式，如下所示。

1）失效概率pi估計。

根據(jù)不完全Beta分布函數(shù)，得其先驗密度函數(shù)為：

其似然函數(shù)為：

根據(jù)Bayes定理，得其后驗分布密度函數(shù)為：

在平方損失下，對后驗分布取期望得到pi的估計為：

取a=1，簡化式(15)，可得：

2）修正失效概率上限值pu。

根據(jù)截尾時間正序ti(t1,t2, …,tk)，設(shè)初始時刻t1的失效概率符合Beta分布，即失效概率取值范圍為(0, 1)。計算t1時刻失效概率上限值pu1，將其作為下一時刻t2的失效概率下限值，然后得到t2時刻失效概率上限值pu2并作為t3時刻失效概率下限值，依次進行至tk時刻，將得到的失效概率上限值集合pu1,pu2, …,puk作為上限，對應(yīng)式(16)中的pu，即：

式中：pui為第i時刻的失效概率上限值。

3）修正失效概率下限值pl。

式中：ascend[·]為升序排列，pli為第i時刻的失效概率下限值。

2.3 參數(shù)估計及可靠度計算

式中：ωi為權(quán)重。根據(jù)文獻[16]可知，以總試驗時間設(shè)計權(quán)重，估計效果最好，則：。

3 驗證分析

3.1 Monte-Carlo仿真數(shù)據(jù)

為了驗證新方法的估計效果，需要進行大量的對比試驗。利用Monte-Carlo 法可以有效地仿真試驗數(shù)據(jù)，便于評估分析。其詳細仿真方法參考文獻[19]的做法：由計算機生成700個服從威布爾分布的隨機數(shù)（產(chǎn)生的隨機數(shù)要遠多于試驗產(chǎn)品樣本數(shù)，本文取產(chǎn)品樣本數(shù)為28個），并將隨機數(shù)從小到大排列，然后以每10 個數(shù)據(jù)為一組，共分為70 組數(shù)據(jù)，對其中前7組數(shù)據(jù)中的首個數(shù)據(jù)向下取整，可得到7組截尾時間，樣本在截尾時間下均沒有發(fā)生失效，各截尾時刻的樣本數(shù)隨組號依次遞減。

通過前文分析可知，考慮到長壽命機械產(chǎn)品的特性，取β=3，η=2 500，得到一組無失效數(shù)據(jù)，如表2所示。

表2 仿真的無失效數(shù)據(jù)Table 2 Simulation of zero-failure data

3.2 仿真結(jié)果對比分析

分別采用傳統(tǒng)Bayes、E-Bayes、多層Bayes 和雙修正多層Bayes對表2數(shù)據(jù)進行運算，得到失效概率估計值，如表3所示。由表可知，E-Bayes和多層Bayes方法下失效概率估計值相差不大，后者略優(yōu)，但兩者均優(yōu)于傳統(tǒng)Bayes，而本文方法在1 049 h前最接近真值。

表3 不同方法下的失效概率估計值Table 3 Estimation of failure probability under different methods

將表3數(shù)據(jù)代入式(19)，并取超參數(shù)c=5，得到各方法下的β和η，并計算其與真值之間的相對誤差Emr.β和Emr.η，結(jié)果如表4 所示。由表可知：E-Bayes和多層Bayes的估計誤差均在40%左右；傳統(tǒng)Bayes對η的估計效果較好，但對β的估計較差；本文方法對β和η的估計誤差均較低并且較均衡，控制在6.3%左右。

表4 不同方法下的參數(shù)估計及其相對誤差Table 4 Parameter estimation and relative error under different methods

將表4 數(shù)據(jù)代入式(20)，得到各方法下的可靠度，如圖2所示。由圖可知：采用本文方法得到的可靠度更加接近真值，且其魯棒性更好；在1 200 h內(nèi)任意時刻的可靠度均優(yōu)于其他Bayes方法。

圖2 不同方法下的可靠度Fig.2 Reliability with different methods

3.3 工程實例分析

航空陀螺儀是一種非常精密的導(dǎo)航儀器，其中軸承是其關(guān)鍵的機械零部件，軸承失效會直接導(dǎo)致陀螺儀發(fā)生工作故障。因此，對陀螺儀的軸承進行可靠性評估是極為重要的。采用文獻[16]提供的某軸承無失效數(shù)據(jù)，如表5所示，分別用傳統(tǒng)Bayes、E-Bayes、多層Bayes 和本文方法進行可靠性分析。另外，根據(jù)相關(guān)工程部門提供的信息，該軸承在1 300 h內(nèi)的可靠度為100%。

表5 某軸承的無失效數(shù)據(jù)Table 5 Zero-failure data of a certain bearing

估計精度受超參數(shù)c的影響，故討論當(dāng)c=4～7時各方法的估計效果和穩(wěn)健性。c=4～7 時軸承的可靠度如圖3所示。

圖3 在c = 4～7時軸承的可靠度Fig.3 Reliability of bearings when c = 4～7

由圖3可知：

1）當(dāng)c=4～7時，從總體上來看，4種方法按可靠度估計精度從大到小的排序依次是本文方法、E-Bayes、多層Bayes 和傳統(tǒng)Bayes，本文方法估計的結(jié)果更加接近實際值，且其具有良好的穩(wěn)健性。

2）當(dāng)c=4時，本文方法在1 100 h以前的可靠度估計效果好于E-Bayes和多層Bayes，1 100 h以后逐漸出現(xiàn)差于多層Bayes和E-Bayes的趨勢。

3）在試驗組數(shù)不變的情況下，c值越大，各方法的估計效果越好。

4 雙修正多層Bayes的適用性分析

以上驗證結(jié)果表明了本文方法的估計效果較好，但在實驗過程中存在著一定的偶然性和局限性，主要體現(xiàn)在利用Monte-Carlo 法生成呈威布爾分布的無失效數(shù)據(jù)時須提前設(shè)置β和η，因此僅能體現(xiàn)本文方法在某一特定條件下的優(yōu)勢。為了驗證本文算法的適用性，分析當(dāng)β=1.4～5.0，η=1 500～13 500時雙修正多層Bayes結(jié)合加權(quán)最小二乘法的參數(shù)估計效能。利用Monte-Carlo 仿真試驗，循環(huán)10 000次，輸出結(jié)果為相對誤差平均值，具體試驗流程如下圖4所示。

圖4 適用性分析的試驗流程Fig.4 Experimental process for applicability analysis

運用定量分析法討論本文方法隨威布爾雙參數(shù)變化的可靠度估計能力,并與E-Bayes、多層Bayes方法進行對比。設(shè)置η=2 000，得到β變化時各方法的估計效果，結(jié)果如圖5所示。設(shè)置β=3，得到η變化時各方法的估計效果，結(jié)果如圖6所示。

圖5 各方法估計β的相對誤差Fig.5 Relative error of each method to estimate β

圖6 各方法估計η的相對誤差Fig.6 Relative error of each method to estimate η

由圖5和圖6可知：β會影響各Bayes的估計效果，而η幾乎不會產(chǎn)生影響；多層Bayes和E-Bayes的估計誤差隨著β的增大而增大，在β=1.5 左右時估計效果最好，相對誤差控制在20%以下，較本文方法更有優(yōu)勢；本文方法在β＞2.2時的估計效果優(yōu)于E-Bayes和多層Bayes，在β=3左右時估計效果最好，相對誤差在10%以下，在β＞3時Emr-β出現(xiàn)增大的趨勢，而Emr-η基本不變，總體效果仍顯著優(yōu)于E-Bayes 和多層Bayes，且較它們具有更寬的適用區(qū)間。

5 結(jié) 論

為了提高對威布爾分布無失效數(shù)據(jù)的可靠度點估計精度，本文提出了雙修正多層Bayes方法來估計失效概率并完成可靠度評估，得出如下結(jié)論：

1）雙修正多層Bayes方法能有效改善失效概率的估計精度，尤其在結(jié)合加權(quán)最小二乘法后，能將威布爾雙參數(shù)估計的相對誤差控制在10%以下。

2）以雙修正多層Bayes方法為核心，結(jié)合加權(quán)最小二乘法，能夠有效進行參數(shù)估計和可靠度估計，并具有較高的精度和良好的穩(wěn)健性。

3）通過適用性分析發(fā)現(xiàn)，雙修正多層Bayes的估計精度更高，且適用范圍更寬，尤其對于β＞2.2的威布爾分布數(shù)據(jù)，較主流Bayes方法有顯著優(yōu)勢。研究結(jié)果可以為其他壽命分布無失效數(shù)據(jù)的可靠性評估提供參考。