基于貝葉斯推斷的威布爾型產(chǎn)品可靠性評估

0 引言

以運(yùn)載火箭、風(fēng)電機(jī)組、新型裝備等為代表的產(chǎn)品由于可靠性試驗(yàn)周期短，實(shí)驗(yàn)成本高，實(shí)驗(yàn)數(shù)據(jù)短缺，在對其進(jìn)行可靠性評估等量化工作時面臨小子樣問題[1]。最大似然估計法等傳統(tǒng)可靠性評估方法在面對小子樣可靠性評估問題時容易產(chǎn)生較大偏差[2]。而貝葉斯推斷能融合并利用各種先驗(yàn)信息（包括歷史數(shù)據(jù)、相似產(chǎn)品數(shù)據(jù)以及專家經(jīng)驗(yàn)）與試驗(yàn)信息，并用試驗(yàn)信息不斷更新先驗(yàn)信息，最終得到可靠性評估模型中待估參數(shù)和可靠性指標(biāo)的估計值與置信區(qū)間，能有效控制小樣本數(shù)據(jù)帶來的不確定性[3]。因此，貝葉斯推斷在小子樣可靠性評估中應(yīng)用廣泛。

在已有研究多集中于指數(shù)分布型產(chǎn)品或簡化威布爾型產(chǎn)品的情況下[4-6]，對更一般雙參數(shù)的威布爾型產(chǎn)品建立基于貝葉斯推斷的可靠性評估模型。威布爾分布函數(shù)無共軛先驗(yàn)分布，假設(shè)待估參數(shù)服從Gamma分布[7]，保證其正定性，利用先驗(yàn)信息與專家經(jīng)驗(yàn)確定Gamma分布中超參數(shù)，得到待估參數(shù)先驗(yàn)概率密度函數(shù)，結(jié)合似然函數(shù)生成后驗(yàn)概率密度函數(shù)，應(yīng)用馬爾可夫鏈蒙特卡洛法具體為NUTS抽樣模擬待估參數(shù)的后驗(yàn)概率密度函數(shù)，求得待估參數(shù)估計值與置信區(qū)間。以某型號風(fēng)力發(fā)電機(jī)齒輪箱的8次故障數(shù)據(jù)為算例，證明該方法的正確性與有效性，并完成對其可靠性評估。

1 貝葉斯推斷

貝葉斯推斷能將先驗(yàn)信息與觀察到的故障數(shù)據(jù)相結(jié)合，以推斷待估參數(shù)ξ（ξ1，ξ2，…，ξn）的數(shù)學(xué)特征。貝葉斯推斷的核心是貝葉斯定理，貝葉斯定理可以實(shí)現(xiàn)通過收集到的數(shù)據(jù)y（y1，y2，…，yn）對已有信息更新[8，9]。如下所示：

其中y（y1，y2，…，yn），表示觀測到隨機(jī)變量y的數(shù)據(jù)向量；p（ξ｜y）為后驗(yàn)密度函數(shù)；p（ξ）為先驗(yàn)密度函數(shù)；f（y｜ξ）為數(shù)據(jù)的抽樣密度函數(shù)即似然函數(shù)；m（y）為數(shù)據(jù)的邊緣密度函數(shù)。根據(jù)貝葉斯定理，后驗(yàn)分布可成比例表示為似然函數(shù)與先驗(yàn)分布的乘積，即

其中m（y）與參數(shù)ξ無關(guān)，被視為比例常數(shù)的一部分。獲得觀測數(shù)據(jù)前，對參數(shù)ξ（ξ1，ξ2，…，ξn）的信息集中于p（ξ），在獲得觀測數(shù)據(jù)y（y1，y2，…，yn）后生成似然函數(shù)，并更新生成后驗(yàn)分布。

1.1 Weibull分布

Weibull分布作為常用可靠性分析統(tǒng)計分布，對機(jī)械產(chǎn)品的故障時間數(shù)據(jù)擬合能力強(qiáng)[10，11]，故障時間t的概率密度函數(shù)為

其中，β是形狀參數(shù)；θ是尺度參數(shù)，也稱為特征壽命。待估參數(shù)為β與θ。

可靠度函數(shù)為

累積分布函數(shù)為

平均故障間隔時間為

1.2 似然函數(shù)

對于故障時間T（t1，t2，…，tN），T~Weibull（θ，β），似然函數(shù)為各故障時間概率密度函數(shù)乘積，即

1.3 待估參數(shù)先驗(yàn)分布與后驗(yàn)分布

當(dāng)Weibull分布的形狀參數(shù)和尺度參數(shù)均未知時，不存在適用的共軛先驗(yàn)分布。對于θ與β，假設(shè)其先驗(yàn)分布為Gamma分布，且兩個參數(shù)相互獨(dú)立，即

θ~Gamma（α1，λ1），β~Gamma（α2，λ2）

α1，α2＞0為形狀參數(shù)；λ1，λ2＞0為尺寸參數(shù)。對于超參數(shù)α1，α2＞0與λ1，λ2＞0可根據(jù)產(chǎn)品可靠性先驗(yàn)信息，以矩估計與專家經(jīng)驗(yàn)相結(jié)合的方法求解。

此時θ概率密度函數(shù)（條件先驗(yàn)）為f（θ｜α1，λ1）=，均值對θ概率密度函數(shù)（條件先驗(yàn)）化簡可得p（θ｜α1，λ1）∝θα1-1exp（-λ1θ）。θ條件后驗(yàn)密度函數(shù)為：

對于以上條件后驗(yàn)密度函數(shù)由于無法用解析方法求得，需要應(yīng)用數(shù)值模擬方法—馬爾可夫鏈蒙特卡洛方法[12，13]。

2 馬爾可夫鏈蒙特卡洛法

馬爾可夫鏈蒙特卡洛法（MCMC）也稱為馬爾可夫鏈蒙特卡洛抽樣器，其將馬爾可夫鏈與蒙特卡洛法相結(jié)合。其 中 馬爾可夫鏈?zhǔn)侵敢幌盗须S 機(jī)樣本ξ（ξ1，ξ2，…，ξn），P（ξn+1｜ξ1，ξ2，ξ3，…，ξn）=P（ξn+1｜ξn）。即ξn+1取值只與當(dāng)前取值ξn有關(guān)，而與之前取值無關(guān)。在馬爾可夫鏈中相鄰樣本（ξi，ξi+1）之間可以發(fā)生轉(zhuǎn)換，其轉(zhuǎn)換概率P（ξi→ξi+1）由轉(zhuǎn)換概率分布函數(shù)確定。在齊次馬爾可夫鏈中，生成的樣本分布將逐漸收斂，并趨于平穩(wěn)[14，15]。在貝葉斯推斷背景下，該平穩(wěn)分布即為待估參數(shù)后驗(yàn)分布[16，17]。但馬爾可夫鏈的初始樣本并不來源于平穩(wěn)分布，因此也不能代表該平穩(wěn)分布，該過程稱為老煉階段（burn-in），應(yīng)從馬爾可夫鏈中剔除。

設(shè)老練期生成樣本數(shù)Nburn-in，馬爾可夫鏈生成總樣本數(shù)N。當(dāng)算法進(jìn)行了Nburn-in次后，產(chǎn)生的仿真序列可以認(rèn)為是來自后驗(yàn)分布的樣本，而Nburn-in次之前的樣本屬于算法的老煉階段，不能認(rèn)為來自后驗(yàn)分布。因此，MCMC創(chuàng)造了一個馬爾可夫鏈，該馬爾可夫鏈的不變分布即為所求的后驗(yàn)概率分布。常用MCMC算法有Metropolis-Hastings算 法、Gibbs抽 樣 與Hamiltonian蒙特卡洛方法[18]。其中Hamiltonian蒙特卡羅（Hamiltonian Monte Carlo，HMC）方法也稱為混合蒙特卡羅（Hybrid Monte Carlo）方法，其采用系統(tǒng)動力學(xué)而不是概率分布來建議馬爾可夫鏈中的未來狀態(tài)，使得馬爾可夫鏈能更有效地探索目標(biāo)分布從而實(shí)現(xiàn)更快的收斂。無掉頭抽樣器（No-U-Turn Sampler，NUTS）作為HMC的拓展，使用成對平均算法[19]，可以在完全不進(jìn)行任何手動調(diào)整的情況下運(yùn)行，生成的樣本至少與精細(xì)的手動調(diào)整HMC一樣。

對威布爾型產(chǎn)品，參數(shù)向量中包含2個元素（θ，β），由1.3節(jié)知每個元素的條件后驗(yàn)分布為p（θ｜T，β，α1，λ1），p（β｜T，θ，α2，λ2）。

3 算例分析

以某型號風(fēng)力發(fā)電機(jī)齒輪箱的8次故障數(shù)據(jù)為算例[20]，驗(yàn)證上述可靠性評估模型。風(fēng)電機(jī)組齒輪箱故障狀態(tài) 數(shù) 據(jù)（4707，6167，9887，9991，13712，16543，18007，20793），均為停機(jī)故障，單位：時。使用最小二乘法計算得到，θ=19000，β=2.3。

對于先驗(yàn)分布，應(yīng)用矩估計法并借助專家經(jīng)驗(yàn)，可確定以下方程式：

解得，λ1=3.6；α1=64800，λ2=1，α2=1。將結(jié)果代入θ條件后驗(yàn)密度函數(shù)與β條件后驗(yàn)密度函數(shù)。

3.1 NUTS抽樣

采用基于Python平臺的PyMC3進(jìn)行NUTS抽樣，仿真環(huán)境為PyCharm Edu，抽樣次數(shù)設(shè)置為5000次，Nburn-in=1000，抽樣結(jié)果如表1所示。

表1 NUTS抽樣結(jié)果

即θ后驗(yàn)均值17998.02，標(biāo)準(zhǔn)差72.72；β后驗(yàn)均值2.11，標(biāo)準(zhǔn)差0.63。

生成后驗(yàn)密度直方圖和最高密度區(qū)域圖如圖1。

圖1 θ、β后驗(yàn)密度直方圖與最高密度區(qū)域

生成密度估計圖和痕跡圖如圖2、圖3。

圖2 θ密度估計圖和痕跡圖

圖3 β密度估計圖和痕跡圖

3.2 結(jié)果分析

文獻(xiàn)[20]分別使用最小二乘法根據(jù)8次故障停機(jī)數(shù)據(jù)及文獻(xiàn)改進(jìn)方法根據(jù)8次停機(jī)故障數(shù)據(jù)加12次非停機(jī)故障數(shù)據(jù)對雙參數(shù)進(jìn)行估計，結(jié)果如表2所示。

表2 三種方法參數(shù)估計對比

通過對比可知，提出基于貝葉斯推斷的威布爾型產(chǎn)品可靠性評估方法在利用8次停機(jī)故障數(shù)據(jù)情況下，估計所得θ均值介于最小二乘法與文獻(xiàn)[20]改進(jìn)方法。在可利用故障數(shù)據(jù)一致情況下，本文方法得到的尺度參數(shù)θ估計值比最小二乘法更精確，相較文獻(xiàn)[20]方法更加合理。對于形狀參數(shù)β，估計值大于1，說明失效率隨時間增加，根據(jù)本文方法得出的故障規(guī)律與最小二乘法、文獻(xiàn)[20]改進(jìn)方法一致。

通過生成抽樣能量圖，可以看出能量（Marginal Energy Distribution）及能量轉(zhuǎn)換（Energy Transition Distribution）兩個分布范圍接近，寬窄相近，說明NUTS抽樣中失去的信息不重要，后驗(yàn)估計準(zhǔn)確。

3.3 可靠性評估

圖4 抽樣能量圖

將NUTS抽樣產(chǎn)生的形狀參數(shù)與尺度參數(shù)帶入可靠性評估模型（4）（5）（6）（7），可得故障時間t的概率密度函數(shù)為：

可靠度函數(shù)為：

R（t）=exp[-（t/17998）2.11]

累積分布函數(shù)為：

F（t）=1-exp[-（t/17998）2.11]

平均故障間隔時間為：

由上述可靠性模型計算得出某型號風(fēng)力發(fā)電機(jī)齒輪箱的平均故障間隔時間為15943h，完成了可靠性評估。

4 結(jié)論

本文針對更一般的雙參數(shù)威布爾型產(chǎn)品建立基于貝葉斯推斷的可靠性評估模型，在無共軛先驗(yàn)分布情況下，均假設(shè)其待估參數(shù)服Gamma從先驗(yàn)分布，推導(dǎo)出待估參數(shù)后驗(yàn)概率密度，而非簡化為指數(shù)分布型可靠性模型。充分利用已有數(shù)據(jù)與專家經(jīng)驗(yàn)，應(yīng)用MCMC法中NUTS抽樣求解待估參數(shù)點(diǎn)估計值與置信區(qū)間，實(shí)現(xiàn)復(fù)雜密度函數(shù)的仿真模擬。經(jīng)過算例分析可知該模型的輸出結(jié)果正確可靠，既擴(kuò)大貝葉斯推斷的應(yīng)用范圍又能保證求解精度。模型通過調(diào)整先驗(yàn)分布可以求解不同威布爾型產(chǎn)品的待估參數(shù)與可靠度函數(shù)，具有可重用性。對于不同威布爾型產(chǎn)品，借助先驗(yàn)信息與專家經(jīng)驗(yàn)確定待估參數(shù)先驗(yàn)分布后，通過修改PyMC3中參數(shù)信息，可用于求解待估參數(shù)，有助于對小子樣產(chǎn)品的可靠性評估與可靠性管理。

圖5 故障時間概率密度函數(shù)

圖6 可靠度函數(shù)

在本文基礎(chǔ)上，可以進(jìn)一步研究融合多源數(shù)據(jù)的基于貝葉斯推斷可靠性評估，量化專家經(jīng)驗(yàn)控制主觀性，充分利用各種先驗(yàn)信息，進(jìn)行更準(zhǔn)確的小子樣產(chǎn)品可靠性評估與優(yōu)化。