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

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

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

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

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

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

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

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

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

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

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

1.2 威布爾分布

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

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

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

失效率函數(shù)為：

可靠度函數(shù)為：

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

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

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

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

2 無(wú)失效威布爾可靠度估計(jì)

2.1 失效概率pi的多層Bayes估計(jì)

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

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

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

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

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

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

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

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

2.2 失效概率pi的雙修正多層Bayes估計(jì)

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

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

1）失效概率pi估計(jì)。

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

其似然函數(shù)為：

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

在平方損失下，對(duì)后驗(yàn)分布取期望得到pi的估計(jì)為：

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

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

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

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

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

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

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

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

3 驗(yàn)證分析

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

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

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

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

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

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

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

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

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

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

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

3.3 工程實(shí)例分析

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

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

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

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

由圖3可知：

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

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

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

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

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

圖4 適用性分析的試驗(yàn)流程Fig.4 Experimental process for applicability analysis

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

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

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

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

5 結(jié) 論

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

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

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

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