陳紅霞，張俊峰*，馬愛博，李宏悅，李晨光

(1. 內(nèi)蒙古工業(yè)大學(xué)機(jī)械工程學(xué)院 呼和浩特 010051；2. 內(nèi)蒙古自治區(qū)先進(jìn)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室 呼和浩特 010051)

重型數(shù)控機(jī)床對現(xiàn)代軍工、航天、船舶、軌道交通、發(fā)電設(shè)備等重工業(yè)發(fā)展有重要意義[1]。目前，研究人員大多借鑒中小型數(shù)控機(jī)床的可靠性模型對重型數(shù)控機(jī)床進(jìn)行可靠性分析[2]。重型數(shù)控機(jī)床具有加工工況復(fù)雜、故障溯源困難、維護(hù)成本較高、故障樣本數(shù)據(jù)匱乏等特點(diǎn)，這對重型數(shù)控機(jī)床的可靠性評估及后續(xù)研究都造成極大困難[3]。小樣本數(shù)據(jù)下的重型機(jī)床可靠性評估，是目前可靠性研究的一個(gè)難題[4]。

數(shù)控機(jī)床的可靠性研究最早開始于20 世紀(jì)70 年代[5]，前蘇聯(lián)專家對建立機(jī)床參數(shù)模型和機(jī)床工藝可靠性方面進(jìn)行了相關(guān)研究[6]。近些年來，一些專家將貝葉斯理論引入小樣本條件下的重型數(shù)控機(jī)床可靠性建模[7]，根據(jù)先驗(yàn)信息和主觀經(jīng)驗(yàn)進(jìn)行模型參數(shù)估計(jì)[8]。文獻(xiàn)[9]采用威布爾雙參數(shù)和三參數(shù)分布對重型數(shù)控機(jī)床零部件或系統(tǒng)建模[9]。文獻(xiàn)[10]采用極大似然估計(jì)法和Edgeworth 級(jí)數(shù)法對重型數(shù)控機(jī)床的故障模型進(jìn)行參數(shù)估計(jì)，并解釋參數(shù)估計(jì)結(jié)果與數(shù)據(jù)擬合的關(guān)系。文獻(xiàn)[11]利用第二類極大似然方法得到折合因子的貝葉斯估計(jì)[11]。文獻(xiàn)[12]在馬爾科夫鏈蒙特卡洛方法(Markov Chain Monte Carlo, MCMC)方法的基礎(chǔ)上針對數(shù)據(jù)的多源性，建立多源異種數(shù)據(jù)融合模型，為小樣本的可靠性建模提供新思路。上述研究對小樣本數(shù)據(jù)下的可靠性建模問題提出了不同的方法，給重型數(shù)控機(jī)床可靠性建模及評估奠定了理論基礎(chǔ)。

本文基于改進(jìn)的貝葉斯方法對重型數(shù)控機(jī)床進(jìn)行可靠性建模與評估。首先，選取雙參數(shù)的威布爾分布對機(jī)床故障數(shù)據(jù)進(jìn)行擬合，采用極大似然參數(shù)估計(jì)法和貝葉斯參數(shù)估計(jì)法對可靠性模型求解。其次，對傳統(tǒng)的貝葉斯估計(jì)法進(jìn)行改進(jìn)，得到分層貝葉斯模型，通過對比分析及仿真實(shí)驗(yàn)，證明改進(jìn)后的貝葉斯模型在小樣本數(shù)據(jù)條件下更具有優(yōu)越性。最后，利用分層貝葉斯模型對該機(jī)床進(jìn)行可靠性評估，得到其理論故障間隔時(shí)間。

1 數(shù)控機(jī)床的可靠性模型參數(shù)估計(jì)

1.1 傳統(tǒng)可靠性模型

一般地，對系統(tǒng)建立可靠性模型均采用雙參數(shù)威布爾分布。以重型數(shù)控機(jī)床故障間隔時(shí)間為參數(shù)，經(jīng)過變換得到服從威布爾分布的可靠性模型。此時(shí)，機(jī)床可靠性服從威布爾分布的概率密度函數(shù)和概率分布函數(shù)分別為：

式中，m和η 分別為威布爾分布的形狀參數(shù)和尺度參數(shù)；t為機(jī)床的故障時(shí)間。

1.2 極大似然估計(jì)法

一般地，采用極大似然估計(jì)對分布函數(shù)進(jìn)行參數(shù)估計(jì)。當(dāng)已知分布模型中有未知參數(shù)θ1,θ2,···,θk時(shí)，可以建立似然函數(shù)表達(dá)式：

式中，ti為實(shí)際采集的數(shù)據(jù)子樣(i=1, 2, ···,n)。

在威布爾分布模型中，未知參數(shù)為m和η，將式(1)的威布爾分布概率密度函數(shù)代入式(3)中，建立的極大似然函數(shù)為：

1.3 基于貝葉斯理論的參數(shù)估計(jì)

由貝葉斯參數(shù)估計(jì)理論可知，在對分布函數(shù)的參數(shù)估計(jì)中，所有未知參數(shù)都被認(rèn)為是隨機(jī)變量。求解待估參數(shù)值時(shí)，需要先根據(jù)觀測數(shù)據(jù)和主觀經(jīng)驗(yàn)設(shè)定其模型的先驗(yàn)分布，進(jìn)而求出參數(shù)估計(jì)結(jié)果。其中，先驗(yàn)分布的選擇優(yōu)劣會(huì)直接影響模型參數(shù)估計(jì)的精度。

基于貝葉斯理論對重型數(shù)控機(jī)床建立可靠性模型，具體思路為在研究對象樣本數(shù)據(jù)和先驗(yàn)信息的基礎(chǔ)上，建立其貝葉斯模型，計(jì)算模型參數(shù)的后驗(yàn)分布并得出相應(yīng)的概率分布。

結(jié)合貝葉斯定理分析可知，分布模型參數(shù)的后驗(yàn)分布可以表示為：

式中，θ 為待估計(jì)的參數(shù)；y為基于觀測數(shù)據(jù)的信息；P(θ)為 參 數(shù) θ 的先驗(yàn) 概 率；P(y|θ)是設(shè)定 參 數(shù)θ 條件下的數(shù)據(jù)所服從的后驗(yàn)分布；P(y|θ)/P(y)為事件y對參數(shù)設(shè)定的支持程度，即似然函數(shù)。

因?yàn)檎龖B(tài)分布是威布爾分布的共軛先驗(yàn)分布模型，所以設(shè)定威布爾分布的參數(shù)m和η 服從正態(tài)分布，這里用N來表示。

經(jīng)過推導(dǎo)，得到重型機(jī)床基于貝葉斯理論的可靠性模型，如式(10)所示：

1.4 層次貝葉斯模型

為了優(yōu)化貝葉斯參數(shù)估計(jì)的方法，將待估參數(shù)的先驗(yàn)分布進(jìn)一步分化，得到層次貝葉斯模型。層次貝葉斯模型理論上可以分為N層，但隨著分化層數(shù)的增加，模型的計(jì)算效率逐漸降低。為了保證計(jì)算效率，采用兩層貝葉斯模型進(jìn)行分析。

由一般貝葉斯模型可知，m～N(α1,β1)為參數(shù)m的先驗(yàn)分布， η～N(α2,β2) 為 參數(shù) η的先驗(yàn)分布。將α2～unif(α1′,β1′),β2～unif(α2′,β2′)分別作為參數(shù)m, η 的先驗(yàn)分布，這樣就得到了更進(jìn)一層的貝葉斯模型，這里稱 α2,β2為超參數(shù)，unif 表示均勻分布，表達(dá)如下：

經(jīng)過上式推導(dǎo)，即可得到基于層次貝葉斯理論的可靠性模型。

1.5 貝葉斯模型的求解方法

求解貝葉斯模型的參數(shù)需要計(jì)算高維積分的邊緣后驗(yàn)分布，而且隨著模型組成的復(fù)雜會(huì)使參數(shù)求解變得更加困難。通過引入 MCMC 方法對模型進(jìn)行數(shù)據(jù)模擬、迭代以提升參數(shù)求解的效率。該方法中包含多種抽樣方法，包括Metropolis-Hastings 抽樣方法和 Gibbs 抽樣方法。其中，Metropolis-Hastings抽樣方法更為高效，具體操作流程如圖1 所示。

圖1 蒙特卡洛仿真流程

通過MCMC 法處理得到模型待估參數(shù)的后驗(yàn)分布，取其期望值作為參數(shù)估計(jì)的結(jié)果，基于排序取分位數(shù)法，可以求出待估參數(shù)的置信區(qū)間。

1.6 參數(shù)估計(jì)結(jié)果分析

為了分析參數(shù)估計(jì)結(jié)果的優(yōu)劣，本文采用統(tǒng)計(jì)學(xué)中的均方根誤差值(normalized root mean square error, NRMSE)，NRMSE 值越小則結(jié)果精度更優(yōu)。計(jì)算方法為：

式中：F(ti)是機(jī)床的累計(jì)故障概率，F(xiàn)(ties)是將參數(shù)估計(jì)值代入累積故障概率函數(shù)后所得到的累積故障概率值。

2 建立重型數(shù)控機(jī)床的可靠性模型

重型機(jī)床故障數(shù)據(jù)樣本量較小、數(shù)據(jù)跨度較大、故障發(fā)生具有不確定性，若采用一般統(tǒng)計(jì)學(xué)方法無法對此類問題進(jìn)行有效分析，本文基于貝葉斯理論進(jìn)行分析，可以一定程度優(yōu)化分析結(jié)果。

2.1 故障數(shù)據(jù)預(yù)處理

本文采集某企業(yè)近3 年內(nèi)，某型號(hào)重型數(shù)控機(jī)床的停機(jī)維修記錄共58 條。通過數(shù)據(jù)預(yù)處理，得到該機(jī)床的故障間隔時(shí)間，其中部分結(jié)果為22 h，23 h，24 h，···，912 h，936 h，1 248 h。

2.2 參數(shù)估計(jì)結(jié)果

將故障間隔時(shí)間代入機(jī)床的雙參數(shù)威布爾分布可靠性模型中，利用最大似然參數(shù)估計(jì)法和貝葉斯參數(shù)估計(jì)法分別對分布模型的參數(shù)進(jìn)行求解。采用最大似然參數(shù)估計(jì)法得到的結(jié)果如表1 所示，用下標(biāo)MLE 來表示。

表1 最大似然參數(shù)估計(jì)結(jié)果

將機(jī)床的故障間隔時(shí)間從大到小排列，代入式(13)所示的中位秩公式：

計(jì)算該重型機(jī)床實(shí)際的累計(jì)失效概率密度，得到其點(diǎn)位分布圖，如圖2 所示。

圖2 樣本的數(shù)據(jù)點(diǎn)分布圖

將參數(shù)估計(jì)結(jié)果代入機(jī)床的可靠性模型，得到基于威布爾分布的機(jī)床擬合累計(jì)失效概率曲線，如圖3 所示，圖中虛線為95%置信區(qū)間的上下限。通過比較兩條曲線的重合度來分析模型的擬合精度。

圖3 最大似然估計(jì)下的威布爾模型分布函數(shù)

由圖3 可知，通過極大似然估計(jì)得到的機(jī)床可靠性模型分布曲線，部分與參考數(shù)據(jù)點(diǎn)位重合，沒有重合的部分也基本分布在置信區(qū)間內(nèi)。

基于貝葉斯理論得到的參數(shù)估計(jì)結(jié)果如表2 所示，其中，用下標(biāo)GBR 表示一般貝葉斯參數(shù)估計(jì)的結(jié)果，用下標(biāo)HBR 表示分層貝葉斯參數(shù)估計(jì)的結(jié)果。

表2 一般貝葉斯和層次貝葉斯參數(shù)估計(jì)結(jié)果

通過對參數(shù)進(jìn)行多層迭代，輸入不同的初始值，構(gòu)造馬爾科夫鏈然后觀察其抽樣迭代過程是否收斂，運(yùn)算迭代過程如圖4、圖5 所示。將α1=1,α2=200,β1=β2=0.5代 入 式(10)，α1=1,β1=0.5,α1′=α2′=0,β1′=100,β2′=1代入式(11)，得到待估參數(shù)的后驗(yàn)分布，如圖6、圖7 所示。

圖4 基于一般貝葉斯模型的仿真迭代軌跡圖

圖5 基于分層貝葉斯模型的仿真迭代軌跡圖

圖6 基于一般貝葉斯模型的參數(shù)后驗(yàn)分布概率密度函數(shù)

圖7 基于分層貝葉斯模型的參數(shù)后驗(yàn)分布概率密度函數(shù)

從待估參數(shù)的后驗(yàn)分布圖中可以看出，分層貝葉斯參數(shù)估計(jì)的結(jié)果更加集中。將參數(shù)估計(jì)結(jié)果分別代入式(1)、式(13)中，得到圖8、圖9 所示的分布曲線。由此可知，貝葉斯模型與傳統(tǒng)的模型相比，其擬合曲線與真實(shí)數(shù)據(jù)的重合度更高。

圖8 一般貝葉斯估計(jì)下的曲線分布

圖9 分層貝葉斯估計(jì)下的曲線分布

2.3 結(jié)果分析

將機(jī)床的3 種可靠性模型代入式(12)，計(jì)算得到NRMSE 值，結(jié)果分別為MLE=0.097 1，GBR=0.121 0，HBR=0.090 0。對比結(jié)果可知：分層貝葉斯模型相對于其他兩個(gè)模型來說誤差值更小，數(shù)據(jù)擬合效果更優(yōu)；基于極大似然估計(jì)法模型的NRMSE值低于一般貝葉斯模型的NRMSE 值。雖然此處傳統(tǒng)可靠性模型精度比一般貝葉斯模型更優(yōu)，可是隨著經(jīng)驗(yàn)信息的豐富，優(yōu)化參數(shù)先驗(yàn)信息的設(shè)定后，可以逐漸改善參數(shù)估計(jì)的精確性。

2.4 不同樣本量下的3 種參數(shù)估計(jì)方法對比

改變數(shù)據(jù)樣本的容量，分別使用上述3 種參數(shù)估計(jì)方法，經(jīng)過計(jì)算求出不同樣本容量下的標(biāo)準(zhǔn)均方根誤差值，結(jié)果如圖10 所示。

圖10 不同樣本量下的建模NRMSE 值

從圖中可以看出，隨著樣本容量的增加，經(jīng)過3 種參數(shù)估計(jì)方法得到的可靠性模型NRMSE 值均在下降，說明建模誤差在降低，精度越來越好，當(dāng)數(shù)據(jù)量趨于大數(shù)據(jù)樣本時(shí)，貝葉斯方法依然有良好的可靠性。但當(dāng)數(shù)據(jù)量較小時(shí)，其NRMSE 值更小，分層貝葉斯參數(shù)估計(jì)法有顯著優(yōu)勢。特別是在樣本量趨于0 時(shí)，基于極大似然估計(jì)建模得到的NRMSE 值近似于1，而基于貝葉斯理論建模得到NRMSE 值依然小于1，說明在樣本量很小時(shí)，傳統(tǒng)的方法幾乎失去精確性，而基于貝葉斯理論得到的結(jié)果依然有效。綜上所述，在數(shù)據(jù)量不足時(shí)，分層貝葉斯模型的擬合效果更優(yōu)。

2.5 仿真驗(yàn)證

基于matlab 軟件，將服從形狀參數(shù)Beta=2 及尺度參數(shù)Eta=1 的雙參數(shù)威布爾分布函數(shù)隨機(jī)篩選10 組數(shù)據(jù)，把產(chǎn)生的小樣本量代入到3 種參數(shù)估計(jì)方法中，得到表3 的結(jié)果。

表3 仿真估計(jì)結(jié)果對比

經(jīng)過仿真分析也可以驗(yàn)證分層貝葉斯在進(jìn)行小樣本參數(shù)估計(jì)方面的優(yōu)勢，其估計(jì)的參數(shù)相對來說最接近真實(shí)值。

3 可靠性評估

平均故障間隔時(shí)間(mean time between failures,MTBF)是常用的可靠性評價(jià)指標(biāo)之一，是指產(chǎn)品從發(fā)生故障到下一次發(fā)生故障的平均工作時(shí)間，其表達(dá)式為：

式中，f(t)為概率密度函數(shù)。

將層次貝葉斯參數(shù)估計(jì)結(jié)果代入式(14)后，求得機(jī)床平均故障間隔時(shí)間為M TBF=210.215 h。分析可知，計(jì)算得到的MTBF 值小于所采集數(shù)據(jù)中的最大值1 248.000 h，210.215 h 大于所采集數(shù)據(jù)中的60%，經(jīng)過排序分析處于中上位置。由此可見，平均故障間隔時(shí)間對于機(jī)床預(yù)防故障和維修保養(yǎng)有參考意義，說明了貝葉斯理論對于重型機(jī)床可靠性建模的有效性。

4 結(jié) 束 語

本文在重型數(shù)控機(jī)床故障數(shù)據(jù)小樣本條件下，基于貝葉斯理論建立服從威布爾分布的可靠性模型，并對傳統(tǒng)貝葉斯模型進(jìn)行改進(jìn)，得到分層貝葉斯模型。引入MCMC 算法求解貝葉斯模型的參數(shù)估計(jì)結(jié)果，對比極大似然估計(jì)法與貝葉斯參數(shù)估計(jì)的優(yōu)劣，結(jié)果顯示這兩種參數(shù)估計(jì)結(jié)果近似，但貝葉斯參數(shù)估計(jì)法得到的參數(shù)置信區(qū)間更窄。采用NRMSE 分析參數(shù)估計(jì)結(jié)果精度，隨著樣本數(shù)據(jù)容量增大，傳統(tǒng)參數(shù)估計(jì)法與貝葉斯參數(shù)估計(jì)法趨于一致，當(dāng)數(shù)據(jù)樣本較小時(shí)，分層貝葉斯模型得到的建模精度更高。對貝葉斯理論模型進(jìn)行仿真實(shí)驗(yàn)，結(jié)果表明貝葉斯參數(shù)估計(jì)是有效的。利用貝葉斯可靠性模型對該重型數(shù)控機(jī)床進(jìn)行可靠性評估，在機(jī)床故障數(shù)據(jù)小樣本條件下，其評估結(jié)果更精確有效，對重型數(shù)控機(jī)床的可靠性研究有重要意義。