李彥強　江杰

摘要：該文借助WINBUGS軟件運用MCMC（MarkovChainMonteCarlo，馬爾科夫鏈蒙特卡洛）估計方法，采用Gibbs抽樣原理，在貝葉斯理論的框架下，結(jié)合某船舶柴油機氣缸蓋歷史的和更新的磨損故障時間數(shù)據(jù)，有效提高了磨損故障概率分布的精度。首先介紹了柴油機對船舶航行安全的重要性及其常見故障類型，確定了采用貝葉斯框架下的MCMC抽樣方法對磨損故障分布進行仿真計算；其次闡述了通過最大似然估計法獲得先驗分布參數(shù)的過程，并說明了MCMC方法的原理及其常用的抽樣方法；最后，引

WINBUGS軟件并采用此軟件進行建模分析，進行基于Gibbs抽樣的MCMC方法的柴油機磨損故障分布計算，以實際算例證明此方法提高了柴油機磨損故障概率分布的計算精度。

關(guān)鍵詞：柴油機；磨損故障；貝葉斯網(wǎng)絡(luò)；MCMC；WinBUGS

中圖分類號：TP311

文獻標識碼：A

文章編號：1009-3044（2017）10-0190-05

1.概述

柴油機是船舶動力裝置的關(guān)鍵設(shè)備，一旦其發(fā)生故障，將會嚴重影響船舶的正常運行，并有可能直接或間接地造成無法估量的經(jīng)濟損失，甚至可能損壞其他關(guān)鍵設(shè)備，嚴重時還會危及到船員和旅客的人身安全。船舶柴油機的工作環(huán)境和條件往往多變而惡劣，由于不穩(wěn)定的環(huán)境與人為操作，船舶柴油機常在負荷與轉(zhuǎn)速多變的工況下運轉(zhuǎn)。柴油機的結(jié)構(gòu)復(fù)雜性決定了柴油機故障類型的復(fù)雜性，負荷與轉(zhuǎn)速多變又會直接導(dǎo)致柴油機的許多零部件常處于高溫、高壓、高速的環(huán)境下工作，加之人為因素導(dǎo)致的潤滑不足、操作不規(guī)范等，柴油機零部件故障多發(fā)、復(fù)雜，具有不確定性和多樣性。

船舶柴油機的故障按照子系統(tǒng)來分，主要分為噴油系統(tǒng)故障、燃料系統(tǒng)故障、渦輪增壓系統(tǒng)故障、冷卻系統(tǒng)故障以及潤滑系統(tǒng)故障等。如果按照導(dǎo)致柴油機故障的原因來分類，那么主要有磨損、沖擊振動、腐蝕、老化等。其中，磨損故障時柴油機的主要故障類型，約占船舶柴油機故障比率的37.5%t31。船舶柴油機磨損故障一旦發(fā)生，若沒有及時、準確地排故和維修，就無法保障船舶動力裝置的正常運行，甚至?xí)?dǎo)致嚴重的安全事故和經(jīng)濟損失。但是，過于頻繁地對船舶柴油機進行檢修，不僅會帶來不必要的人為差錯，還會導(dǎo)致維護費用急劇上升，不符合船舶管理的經(jīng)濟性要求。因此，掌握柴油機磨損故障的發(fā)生規(guī)律、合理地安排檢修工作至關(guān)重要。

國內(nèi)外在這一領(lǐng)域有不少相關(guān)研究，但大多集中在柴油機故障診斷技術(shù)可靠性分析方面，對磨損故障的具體分布研究不夠細化，本文在這方面進行了一定的研究，結(jié)合MC-MC方法有效提高了船舶柴油機磨損故障分布的計算精度。

柴油機磨損故障的失效前平均時間（Mean Time To Fail-tire，MTI'F）的概率分布可以通過歷史的統(tǒng)計數(shù)據(jù)獲得，但由于船舶柴油機工作環(huán)境復(fù)雜多變，根據(jù)歷史數(shù)據(jù)計算得出的概率分布僅能代表綜合環(huán)境下的情況，對于近期具體的工作環(huán)境下的故障分布來說不夠準確。如果能結(jié)合近期的MTTF數(shù)據(jù)，對由歷史數(shù)據(jù)獲得的柴油機磨損故障MTTF分布進行更新，就可以得到更加符合現(xiàn)有實際情況的概率分布，有效提高磨損故障發(fā)生的預(yù)測精度。

在信息更新方面，貝葉斯網(wǎng)絡(luò)（Bayesian Network，BN）具有獨特的優(yōu)勢。貝葉斯網(wǎng)絡(luò)是表示隨機變量間依賴和獨立關(guān)系的網(wǎng)絡(luò)模型，也稱為貝葉斯網(wǎng)、因果概率網(wǎng)絡(luò)或信念網(wǎng)絡(luò)，它是一個有向非循環(huán)的圖形結(jié)構(gòu)，其中節(jié)點代表問題域的隨機變量，節(jié)點間的有向邊代表變量之間的直接依賴關(guān)系，每個節(jié)點都附有一個概率分布。貝葉斯定理可以結(jié)合先驗信息和樣本信息，獲得更新的后驗分布，提高預(yù)測精度，如圖1所示。但是，貝葉斯統(tǒng)計的困難在于，后驗分布是復(fù)雜的、高維的分布，難以進行直接計算，所以隨著計算機技術(shù)與MCMC方法的發(fā)展，后驗分布的計算被大大簡化。

常見的船舶柴油機零部件故障分布規(guī)律有指數(shù)分布、正態(tài)分布、指數(shù)正態(tài)分布和威布爾分布。超載下工作或過熱導(dǎo)致的柴油機故障、人為失誤以及偶然操作不當引起的突發(fā)故障等均適用于指數(shù)分布的故障分布規(guī)律；由幾種相對獨立、作用均勻的微小差異量因素造成的零部件故障，如燃燒室組件、氣缸和活塞等的磨損、噴油系統(tǒng)沉淀故障等，其故障概率分布函數(shù)大多為正態(tài)分布；油水系統(tǒng)和齒輪傳動系統(tǒng)故障、滾珠軸承的磨損故障等服從威布爾分布規(guī)律。由于不同故障分布問題的在本文的研究方法下類似，所以，選取正態(tài)分布進行后續(xù)計算說明。

3.后驗分布參數(shù)獲取

在貝葉斯估計中，后驗均值即為參數(shù)的估計。但在實際計算時，參數(shù)后驗分布往往難以計算，解析法求解幾乎不可能，計算量龐大而過程復(fù)雜。近十幾年來，計算機技術(shù)高速發(fā)展，加之貝葉斯方法的不斷改進完善，特別是MCMC方法以及Win-BUGS（Bayesian inference Using Gibbs Sampling）軟件的發(fā)展和應(yīng)用，原本復(fù)雜的高維計算問題可以通過從后驗密度抽樣取平均的方法方便快速地得到解決。原本傳統(tǒng)蒙特卡洛的基本思想是，當所要求解的問題是某種事件出現(xiàn)的概率，或者是某個隨機變量的期望值時，它們可以通過某種“試驗”的方法，得到這種事件出現(xiàn)的頻率，或者這個隨機變數(shù)的平均值，并用它們作為問題的解。在處理高維、復(fù)雜的問題時，MCMC方法相對于經(jīng)典的蒙特卡洛方法有明顯的優(yōu)越性，可以在非標準形式、各變量之間互相不獨立的情況下進行分布的模擬，在近十幾年內(nèi)被廣泛應(yīng)用于貝葉斯統(tǒng)計、顯著性檢驗等方面。

3.1 MCMC方法

MCMC方法的核心思想是反復(fù)對目標分布進行抽樣，得到目標分布的一組樣本數(shù)據(jù)，然后再利用抽樣得到的樣本對目標分布的性質(zhì)進行分析和討論。也就是說，通過建立一個平穩(wěn)分布為π（x）的馬爾科夫鏈來得到π（x）的樣本，基于這些樣本對目標分布進行各種統(tǒng)計推斷。

在MCMC方法中，可以根據(jù)轉(zhuǎn)移核的不同構(gòu)造方法，將主要的MCMC方法分為三種：

（1）Metropolis抽樣

Metropolis抽樣是最先被使用的MCMC方法，基本思想是從一個對稱分布抽取樣本，計算接收樣本的概率，根據(jù)求得的概率值決定是否接收樣本，再抽取新的樣本，重復(fù)這個過程，直至收斂。

（2）Metropolis-Hastings抽樣

Hastings把最初的Metropolis抽樣進行了進一步推廣，不要求抽取樣本的分布為對稱分布。

（3）Gibbs抽樣

Gibbs抽樣是MCMC方法中應(yīng)用最廣泛的一種，它是Me-tropolis-Hastings抽樣的一種特殊情況。

3.2基于WinBUGS軟件的柴油機零部件磨損故障分布計算

WinBUGS是英國劍橋公共衛(wèi)生研究所的MRC BiostatisticsUnit推出的用MCMC方法進行貝葉斯推理的專用軟件包。WinBUGS可以方便快捷地對各種或常用或復(fù)雜的模型和分布進行Gibbs抽樣，貝葉斯網(wǎng)絡(luò)的各變量節(jié)點之間的關(guān)系也可以通過簡單直觀的有向圖模型進行描述。除此之外，WinBUGS還可以給出指定參數(shù)的Gibbs動態(tài)抽樣圖，用Smoothing方法得到后驗分布的核密度估計圖、抽樣值的自相關(guān)圖等，這使得我們可以對馬爾科夫鏈的收斂性進行直觀的判斷。參數(shù)后驗分布的均值、標準差、蒙特卡洛平均標準誤差（MC error）等信息也可以直接得到。

在解決柴油機磨損故障MTFF的分布計算問題時，可以按照以下幾個步驟進行后驗分布參數(shù)的估計：

1）建立Doodle模型并檢驗

在WinBUGS中，Doodle模型通過節(jié)點、箭頭和平板等元素來進行構(gòu)建。模型中的每個節(jié)點都具備相應(yīng)的詳細屬性，在構(gòu)建Doodle模型時就要指定每個節(jié)點的詳細屬性。

檢驗Doodle模型的目的是檢驗WinBUGS對所建立的模型是否識別，即檢驗建立的模型是否符合軟件要求。當出現(xiàn)報錯信息時，需要修正模型，直至顯示"model is syntactically cot-rect”。

2）載人數(shù)據(jù)、模型編譯和設(shè)定初始值

數(shù)據(jù)輸入可以按照以下方式進行定義：

list（t=c（98，110，77，56，102，89，126，97，74，52），

x=c（6，16，54，3 3，25，4 68 890），N=10）

這種數(shù)據(jù)的表示格式被稱為S-PLUS格式。當輸入的數(shù)據(jù)為多維數(shù)組形式，具體的數(shù)據(jù)輸入格式可以參照軟件中的Help按鈕，其中有進行詳細的介紹。

當數(shù)據(jù)載入后，就可以通過Compile按鈕繼續(xù)進行模型的編譯，對數(shù)據(jù)的完整性和一致性進行檢驗。最終，寫入初始值，通過load inits按鈕進行設(shè)定。

3）變量監(jiān)控和模型迭代

在Inference選項中的Sample按鈕下，設(shè)置所需監(jiān)控的變量，逐個輸入并set。接著在Update中輸入迭代次數(shù)，進行模型迭代。

4）結(jié)果輸出

在之前打開的樣本監(jiān)控工具窗下輸入，便可以得到所有設(shè)置節(jié)點的運行結(jié)果。至此，一個模型的構(gòu)建到運行結(jié)束。

這四個步驟可以用圖1來表示：

下面，本文將針對船舶柴油機氣缸蓋磨損故障的MTTF進行詳細的算例說明，計算其更符合實際、更準確的概率分布。

4.實例應(yīng)用

某船舶柴油機氣缸蓋磨損故障的MTYF歷史數(shù)據(jù)與更新數(shù)據(jù)如表1、表2所示：

根據(jù)前面介紹的極大似然估計，利用現(xiàn)在最常用的質(zhì)量管理統(tǒng)計Minitab，結(jié)合歷史數(shù)據(jù)可以得出氣缸蓋磨損故障MTTF服從Normal（103.4，12.78²），如圖2所示：

Auto C（自相關(guān)函數(shù)）：自相關(guān)函數(shù)用于表征參數(shù)自相關(guān)程度的大小，波動范圍越大，表示自相關(guān)程度越大?？墒褂米韵嚓P(guān)圖來檢測收斂性，如果自相關(guān)程度低，那么在迭代次數(shù)較少時就已收斂。此處自相關(guān)程度較低，基本在迭代次數(shù)20以內(nèi)就已收斂，故后驗分布計算值可信。

stats中的mean代表所求節(jié)點的后驗平均值，通俗來說就是我們要得到的結(jié)果。樣本數(shù)據(jù)越多，后驗估計越準確?？梢杂嬎忝總€參數(shù)的Monte Carlo error（蒙特卡羅標準平均誤差）來評定后驗估計的精度。一般來說，仿真應(yīng)一直進行到參數(shù)的MC error小于樣本SD的5%，即可滿足精度要求。此處結(jié)果均滿足精度要求。

綜合以上，鏈的收斂性可以通過以下幾步來進行判斷：

1）初步可由history和trace所畫的圖來初步觀察，若值在一個平行的區(qū)域內(nèi)而沒有強烈的季節(jié)性，那么鏈可能收斂，此處明顯可能收斂；

2）通過auto C自相關(guān)函數(shù)進一步判斷收斂性，此處自相關(guān)程度低，在迭代次數(shù)很少時就已收斂；

3）接著，可根據(jù)running quantiles運行分位數(shù)來判斷收斂性，圖中分位點穩(wěn)定，暗示已收斂；

4）最后，根據(jù)stats中的MC error比對應(yīng)的后驗標準差要小很多，那么就可以認為鏈收斂。鏈收斂就說明后驗密度被精確估計。

某船舶柴油機氣缸蓋磨損故障MTYF：

先驗分布：Normal（103.4，12.78²）

后驗分布：Normal（98.1，12.49²）

由圖3可知，得到的后驗分布更接近實際總體數(shù)據(jù)的分布，在只利用歷史先驗分布和更新數(shù)據(jù)的條件下提高了分布精度。近期實際運行中得到的MTYF更新數(shù)據(jù)對由歷史數(shù)據(jù)得到的MTTF分布進行了修正，得到的后驗分布更貼近目前的使用實際情況。