胡 宇,唐小峰,文永康,鄒 建

(成都天奧測控技術(shù)有限公司，成都 611731)

0 引言

貝葉斯網(wǎng)絡(luò)方法是由貝葉斯公式發(fā)展而來，貝葉斯公式是用于解決“逆概”的問題一種方法，即通過“結(jié)果”計算“原因”的概率[1]。貝葉斯方法作為概率統(tǒng)計學(xué)的重要分支，廣泛應(yīng)用于機器學(xué)習(xí)、數(shù)據(jù)挖掘、數(shù)據(jù)融合、醫(yī)療診斷、郵件過濾和產(chǎn)品客戶需求分析等場景。故障診斷的過程通常是利用測試“結(jié)果”推理出故障“原因”，這和處理“逆概”問題的概念基本相同，因此，貝葉斯公式和貝葉斯網(wǎng)絡(luò)也經(jīng)常應(yīng)用于故障診斷領(lǐng)域。劉欽文利用貝葉斯網(wǎng)絡(luò)處理傳統(tǒng)多信號模型中故障傳遞和故障檢測中的不確定性問題[2]，文獻(xiàn)[3]使用樸素貝葉斯算法分類提高了滾動軸承故障診斷的準(zhǔn)確率，文獻(xiàn)[4]提出基于重要度分級的概念優(yōu)化傳統(tǒng)貝葉斯網(wǎng)絡(luò)模型有利于故障特性的提取，文獻(xiàn)[5]提出基于貝葉斯網(wǎng)絡(luò)的民機起落架系統(tǒng)故障診斷方法提升了診斷效率和精度。

基于知識的故障診斷方法是目前最為流行的故障診斷方法，這類方法主要是通過故障征兆和系統(tǒng)模型等方式來描述故障現(xiàn)象與故障原因之間的聯(lián)系，系統(tǒng)模型可分為定量模型和定性模型，復(fù)雜系統(tǒng)的精確定性模型往往難以獲得[6]。常用的故障樹和相關(guān)性矩陣都是基于定性模型的故障診斷方法，故障分辨率低[7]，處理復(fù)雜系統(tǒng)多重并發(fā)故障的能力較差。貝葉斯網(wǎng)絡(luò)具有強大的不確定性推理能力，可以利用歷史數(shù)據(jù)建立較為精確的故障概率模型，因此在復(fù)雜系統(tǒng)故障診斷中具有一定的優(yōu)勢[8]。風(fēng)能作為一種無污染、可再生的綠色能源，近年來得到了快速發(fā)展[9]，我國建立了大量的風(fēng)力發(fā)電場并投入使用。經(jīng)過長期運轉(zhuǎn)，風(fēng)力發(fā)電機維護(hù)成本越來越高，對狀態(tài)監(jiān)測和故障診斷的要求也越來越高，傳統(tǒng)手段捉襟見肘。風(fēng)力發(fā)電系統(tǒng)是集機電、機械、電氣和自動控制技術(shù)為一體的復(fù)雜系統(tǒng)，故障原因不確定度較高，可采用貝葉斯網(wǎng)絡(luò)方法進(jìn)行故障診斷。

1 貝葉斯網(wǎng)絡(luò)故障診斷模型

1.1 網(wǎng)絡(luò)結(jié)構(gòu)

貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)用于定性的描述各網(wǎng)絡(luò)節(jié)點之間的關(guān)聯(lián)，在貝葉斯網(wǎng)絡(luò)模型中采用“有向無環(huán)圖-DAG”來實現(xiàn)。貝葉斯網(wǎng)絡(luò)DAG由一系列代表節(jié)點的圓圈和代表關(guān)聯(lián)的箭頭組成，“有向”是指節(jié)點之間的連接是有方向的，“無環(huán)”是指從任一節(jié)點出發(fā)的連接不會再回到該節(jié)點，如圖1所示。

圖1 有向無環(huán)圖

貝葉斯網(wǎng)絡(luò)中的節(jié)點通常代表某個事件或隨機變量的某種狀態(tài)，在故障診斷模型中可以是故障原因、故障單元、故障模塊、故障狀態(tài)、故障現(xiàn)象、和測試結(jié)果等等。節(jié)點間的連線代表了因果關(guān)系，連線起端的節(jié)點為“因”，稱為“父節(jié)點”，連線終端箭頭指向的節(jié)點為“果”，稱為“子節(jié)點”。圖1中節(jié)點“A”、“B”連接并指向“C”，表示事件“A”、“B”有概率導(dǎo)致事件“C”發(fā)生，或者表示變量“A”、“B”的不同狀態(tài)組合會使“C”進(jìn)入某種狀態(tài)。沒有“父節(jié)點”的節(jié)點為根節(jié)點，也稱為一級節(jié)點，一級節(jié)點的子節(jié)點為二級節(jié)點，以此類推n級節(jié)點的子節(jié)點為n+1級節(jié)點，當(dāng)某個節(jié)點擁有多個父節(jié)點時，n取其中最大的那個值。例如：圖1中節(jié)點“A”、“B”為一級節(jié)點，節(jié)點“C”為二級節(jié)點，節(jié)點“D”為三級節(jié)點，節(jié)點“E”既是二級節(jié)點“C”的子節(jié)點，也是三級節(jié)點“D”的父節(jié)點，取最大值“3”再加“1”，因此節(jié)點“E”為四級節(jié)點。

未知的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)難以通過樣本數(shù)據(jù)學(xué)習(xí)的方法獲得，一般采用工程經(jīng)驗進(jìn)行構(gòu)建[6]。以某照明電路為例，電路包括電池、3只保險管和指示燈，保險管“FUSE2”和“FUSE3”并聯(lián)后與保險管“FUSE1”串聯(lián)，當(dāng)電路中所有器件正常時指示燈會點亮，如圖2所示。

圖2 照明電路原理圖

假設(shè)三只保險的連接關(guān)系未知，但經(jīng)驗表明其故障狀態(tài)會影響指示燈的亮/滅狀態(tài)。以變量F1、F2、F3分別代表保險管“FUSE1”、“FUSE2”、“FUSE3”的故障狀態(tài)，變量T代表指示燈“LED”的燃亮狀態(tài)。F1、F2、F3的故障會導(dǎo)致T熄滅，即存在因果關(guān)系，“故障”為因，“熄滅”為果，電路的網(wǎng)絡(luò)結(jié)構(gòu)如圖3所示。

圖3 照明電路貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)圖

1.2 條件概率分布CPD

貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)定性的描述了節(jié)點之間的關(guān)聯(lián)，為了更加準(zhǔn)確的描述節(jié)點的關(guān)聯(lián)程度，還需要用到條件概率參數(shù)。條件概率是指事件B已經(jīng)發(fā)生的條件下，事件A發(fā)生的概率，記為P(A|B)[1]。貝葉斯網(wǎng)絡(luò)中的條件概率是指父節(jié)點發(fā)生或處于某種狀態(tài)下，子節(jié)點發(fā)生或進(jìn)入某種狀態(tài)的概率。當(dāng)子節(jié)點擁有多個父節(jié)點時，需要關(guān)注父節(jié)點全部組合狀態(tài)下的子節(jié)點進(jìn)入各個狀態(tài)的概率，這些概率值的集合即為條件概率分布。圖3所示網(wǎng)絡(luò)結(jié)構(gòu)中， F1、F2、F3是代表“FUSE1”、“FUSE2”、“FUSE3”的故障狀態(tài)的變量，每個變量有兩個取值，“1”為“故障”，“2”為“正?！?。T是代表指示燈“滅”、“亮”的變量，共兩個取值，“1”為“滅”，“2”為“亮”。當(dāng)“FUSE1”、“FUSE2”、“FUSE3”均“故障”時，指示燈“滅”，即F1=1、F2=1、F3=1時T=1的概率為1，根據(jù)歷史數(shù)據(jù)，照明電路的條件概率分布如表1所示。

表1 照明電路CPD

貝葉斯網(wǎng)絡(luò)中的根節(jié)點狀態(tài)存在不確定度，但它沒有父節(jié)點，因此不能利用條件概率來描述根節(jié)點所處狀態(tài)的概率，在貝葉斯統(tǒng)計推斷中通過先驗概率來表達(dá)根節(jié)點的狀態(tài)置信度。照明電路貝葉斯網(wǎng)絡(luò)的根節(jié)點是反映保險管故障狀態(tài)的變量，保險管故障的概率取決于其可靠性，假設(shè)在使用一段時間過后，保險管“FUSE1”故障的概率為10%，無故障的概率為90%，則先驗概率記為：P(F1=1)=0.1，P(F1=2)=0.9。設(shè)三只保險管型號和生產(chǎn)批次一致，可靠性指標(biāo)完全相同，則其先驗概率如表2所示。

表2 保險管先驗概率

1.3 推理算法

貝葉斯推理的基本公式包括：

P(A|B)=P(AB)/P(B)[1]

(1)

P(A)=∑iP(A|Bi)/P(Bi)[1]

(2)

P(Bi|A)=P(A|Bi)/P(Bi)/∑jP(A|Bj)/P(Bj)[1]

(3)

公式(1)為條件概率公式，表示事件B條件下發(fā)生事件A的概率為事件A和事件B同時發(fā)生的概率除以事件B發(fā)生的概率。公式(2)為全概率公式，表示發(fā)生事件A的概率為事件B各種狀態(tài)下發(fā)生事件A的概率乘以事件B該狀態(tài)下概率的乘積之和。公式(3)為貝葉斯公式，由公式(1)、(2)推導(dǎo)得出，用于已知條件概率P(A|B)和先驗概率P(B)時，計算后驗概率P(B|A)。

貝葉斯網(wǎng)絡(luò)模型中最重要的公式為聯(lián)合概率分布公式，用于表述網(wǎng)絡(luò)結(jié)構(gòu)中所有節(jié)點事件同時發(fā)生的概率，可由公式(1)推導(dǎo)得出。圖1所示的網(wǎng)絡(luò)結(jié)構(gòu)的聯(lián)合概率表達(dá)為：

PJ=P(A)*P(B)*P(C|A,B)*P(D|C)*P(E|C,D)

(4)

此求解方法為消元法，是一種精確求解的方法，即通過求和某個節(jié)點各種狀態(tài)的先驗概率與條件概率的乘積，消除聯(lián)合概率表達(dá)式中非待求解的節(jié)點。

例如：圖3所示照明電路的聯(lián)合概率分布為：PJ=P(F1,F2,F3,T) =P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)；

指示燈“LED”不亮?xí)r保險管“FUSE1”故障的概率表述為：P(F1=1|T=1)，求解過程為：

P(F1=1|T=1) =P(F1=1,T=1)/P(T=1) =

∑F1=1,T=1,F2,F3PJ/∑T=1,F1,F2,F3PJ

(5)

∑F1=1,T=1,F2,F3PJ=

∑F1=1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)

(6)

∑F1,T=1,F2,F3PJ=

∑F1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)

(7)

將{F1,T,F2,F3}依次取值{1,1,1,1}、{1,1,1,2}、{1,1,2,1}、{1,1,2,2}代入式(6),得出：

∑F1=1,T=1,F2,F3PJ=[P(F1=1)*P(F2=1)*P(F3=1)*P(T=1|F1=1,F2=1,F3=1)]+[P(F1=1)*P(F2=1)*P(F3=2)*P(T=1|F1=1,F2=1,F3=2)] +[P(F1=1)*P(F2=2)*P(F3=1)*P(T=1|F1=1,F2=2,F3=1)] + [P(F1=1)*P(F2=2)*P(F3=2)*P(T=1|F1=1,F2=2,F3=2)]=0.1*0.1*0.1*1 + 0.1*0.1*0.9*1 + 0.1*0.9*0.1*1 + 0.1*0.9*0.9*1 =0.001+0.009+0.009+0.081 = 0.1。

將{F1,T,F2,F3}依次取值{1,1,1,1}、{1,1,1,2}、{1,1,2,1}、{1,1,2,2}、{2,1,1,1}、{2,1,1,2}、{2,1,2,1}、{2,1,2,2}、代入式(7),得出：

∑F1,T=1,F2,F3PJ=∑F1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)= [P(F1=1)*P(F2=1)*P(F3=1)*P(T=1|F1=1,F2=1,F3=1)] +[P(F1=1)*P(F2=1)*P(F3=2)*P(T=1|F1=1,F2=1,F3=2)] +[P(F1=1)*P(F2=2)*P(F3=1)*P(T=1|F1=1,F2=2,F3=1)] + [P(F1=1)*P(F2=2)*P(F3=2)*P(T=1|F1=1,F2=2,F3=2)] + [P(F1=2)*P(F2=1)*P(F3=1)*P(T=1|F1=2,F2=1,F3=1)] + [P(F1=2)*P(F2=1)*P(F3=2)*P(T=1|F1=2,F2=1,F3=2)] + [P(F1=2)*P(F2=2)*P(F3=1)*P(T=1|F1=2,F2=2,F3=1)] +[P(F1=2)*P(F2=2)*P(F3=2)*P(T=1|F1=2,F2=2,F3=2)] = 0.1*0.1*0.1*1 + 0.1*0.1*0.9*1 + 0.1*0.9*0.1*1 + 0.1*0.9*0.9*1 + 0.9*0.1*0.1*1 + 0.9*0.1*0.9*0 + 0.9*0.9*0.1*0 + 0.9*0.9*0.9*0 = 0.001+0.009+0.009+0.081+0.009+0+0+0 = 0.109。

將上述結(jié)果代入式(5)，得出：

P(F1=1|T=1)=0.1/0.109=91.74%，同理可求出P(F2=1|T=1)=17.43%、P(F3=1|T=1)= 17.43%，即當(dāng)指標(biāo)燈不亮的時候保險管“FUSE1”故障的概率為91.74%，“FUSE2”或“FUSE3”故障的概率為17.43%。

2 風(fēng)力發(fā)電機貝葉斯網(wǎng)絡(luò)MATLAB建模

2.1 網(wǎng)絡(luò)結(jié)構(gòu)建模

風(fēng)力發(fā)電機故障高發(fā)組件、部件包括：冷卻風(fēng)扇、集電環(huán)組件、編碼器組件、接地碳刷組件、發(fā)電機本體、發(fā)電機軸承、測溫電阻等，冷卻風(fēng)扇由電機、電容、電纜、機械結(jié)構(gòu)等零部件組成。集電環(huán)組件由主碳刷、集電環(huán)和其他配件等組成；編碼器組件主要包括編碼器、電纜、機械結(jié)構(gòu)和附件等零部件；接地碳刷組件由接地碳刷和機械結(jié)構(gòu)組成；發(fā)電機本體包括接線端子、繞組、冷卻系統(tǒng)和線纜等零部件；發(fā)電機軸承包括軸承、線纜等零部件；測溫電阻包括繞組測溫電阻、軸承測溫電阻和冷卻系統(tǒng)測溫電阻等。發(fā)電機運行過程中可實時監(jiān)視電機轉(zhuǎn)速、繞組溫度、軸承溫度和冷卻系統(tǒng)溫度等參數(shù)，用于評估系統(tǒng)是否正常運行。風(fēng)力發(fā)電機系統(tǒng)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示。

圖4 風(fēng)力發(fā)電機系統(tǒng)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)圖

風(fēng)力發(fā)電機系網(wǎng)絡(luò)節(jié)點較多，人工處理較繁瑣，因此采用Matlab FullBNT-1.0.7貝葉斯網(wǎng)絡(luò)工具箱進(jìn)行建模。根據(jù)FullBNT的要求將網(wǎng)絡(luò)中的節(jié)點從小到大進(jìn)行編號，且父節(jié)點的編號必須小于子節(jié)點編號，從1開始共50個節(jié)點編號如圖4所示。風(fēng)力發(fā)電機系統(tǒng)網(wǎng)絡(luò)結(jié)構(gòu)Matlab建模部分代碼如下：

N = 50;

dag = zeros(N,N);

N1=1;N2=2;N3=3;N4=4;N5=5;N6=6;N7=7;N8=8;N9=9;N10=10;

N11=11;N12=12;N13=13;N14=14;N15=15;N16=16;N17=17;N18=18;N19=19;N20=20;

N21=21;N22=22;N23=23;N24=24;N25=25;N26=26;N27=27;N28=28;N29=29;N30=30;

N31=31;N32=32;N33=33;N34=34;N35=35;N36=36;N37=37;N38=38;N39=39;N40=40;

N41=41;N42=42;N43=43;N44=44;N45=45;N46=46;N47=47;N48=48;N49=49;N50=50;

dag(N1,N36)=1

dag(N2,N36)=1

dag(N3,N36)=1

dag(N4,N36)=1

dag(N5,N36)=1

………………………………

g=bnt.dag;

draw_graph(g);

其中：代碼dag(Na,Nb)表示節(jié)點a指向節(jié)點b，代碼運行結(jié)果如圖5所示。

圖5 風(fēng)力發(fā)電機系統(tǒng)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)圖

2.2 概率分布建模

在進(jìn)行概率分布建模前，需要對網(wǎng)絡(luò)中各節(jié)點的狀態(tài)取值類型和數(shù)量進(jìn)行定義。風(fēng)力發(fā)電機系統(tǒng)貝葉斯網(wǎng)絡(luò)各節(jié)點均為離散變量，節(jié)點44：“繞組溫度”、節(jié)點45：“軸承溫度”、節(jié)點46：“系統(tǒng)冷卻溫度”均有4個狀態(tài)，分別取值“1”代表“超量程”、“2”代表“正?！薄ⅰ?”代表“過熱”、“4”代表“其他故障”。節(jié)點44、45、46之外的節(jié)點均只有2個狀態(tài)，取值“1”代表“故障”、“2”代表“正?！?。風(fēng)力發(fā)電機系統(tǒng)貝葉斯網(wǎng)絡(luò)中共有35個一級節(jié)點，8個二級節(jié)點，6個三級節(jié)點，1個四級節(jié)點。根據(jù)歷史故障數(shù)據(jù)統(tǒng)計，部分1級節(jié)點先驗概率如表3所示。

表3 1級節(jié)點先驗概率

二級及以上節(jié)點的條件概率的數(shù)量m，取決于其父節(jié)點的數(shù)量n以及節(jié)點狀態(tài)的數(shù)量k，設(shè)kn+1為該節(jié)點自身狀態(tài)數(shù)量：

(8)

例如節(jié)點N41有節(jié)點N22、節(jié)點N23兩個父節(jié)點，節(jié)點N41、N22、N23均只有兩個狀態(tài)，因此N41共有8個條件狀態(tài)參數(shù)；節(jié)點50有8個父節(jié)點，5個節(jié)點有兩種狀態(tài)，3個節(jié)點為4種狀態(tài)，其條件概率參數(shù)達(dá)4 096個。為了將條件概率分布參數(shù)導(dǎo)入FullBNT工具箱，需按一定的規(guī)律對各節(jié)點的條件概率參數(shù)進(jìn)行排序，建立每個節(jié)點的條件概率表。條件概率表用于描述當(dāng)前節(jié)點及相關(guān)節(jié)點處于不同狀態(tài)組合條件下的概率，表的列從左至右以節(jié)點編號從小到大的順序排列，最后一列為條件概率值；表的行從上之下以節(jié)點變量狀態(tài)值從小到大的順序排列。節(jié)點N41“接地碳刷組件”的條件概率表如表4所示，其中第一行代表節(jié)點22、節(jié)點23、節(jié)點41依次取值{1,1,1}時的概率為0.999 99。

表4為完成其余節(jié)點的條件概率表，再將這些概率參數(shù)填入名為FLFDJCPT.xls的EXCEL表格中，便于Matlab導(dǎo)入，部分?jǐn)?shù)據(jù)如表5所示。

表4 節(jié)點N41條件概率表

風(fēng)力發(fā)電機系統(tǒng)條件概率分布建模Matlab部分代碼如下：

discrete_nodes = 1:N;

node_sizes=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 2 2 2 2];

bnt = mk_bnet(dag,node_sizes,'names',{'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10','N11','N12','N13','N14','N15','N16','N17','N18','N19','N20','N21','N22','N23','N24','N25','N26','N27','N28','N29','N30','N31','N32','N33','N34','N35','N36','N37','N38','N39','N40','N41','N42','N43','N44','N45','N46','N47','N48','N49','N50'},'discrete',discrete_nodes);

bnt.CPD{N1} = tabular_CPD(bnt,N1,[0.08258 0.91742]);

bnt.CPD{N2} = tabular_CPD(bnt,N2,[0.00965 0.99035]);

………………………………

bnt.CPD{N35} = tabular_CPD(bnt,N35,[0.00021 0.99979]);

表5 風(fēng)力發(fā)電機條件概率分布

bnt.CPD{N36} = tabular_CPD(bnt,N36,[xlsread('FLFDJCPT.xls','sheet1','A2:A65')]);

bnt.CPD{N37} = tabular_CPD(bnt,N37,[xlsread('FLFDJCPT.xls','sheet1','B2:B33')]);

bnt.CPD{N38} = tabular_CPD(bnt,N38,[xlsread('FLFDJCPT.xls','sheet1','C2:C33')]);

bnt.CPD{N39} = tabular_CPD(bnt,N39,[xlsread('FLFDJCPT.xls','sheet1','D2:D17')]);

bnt.CPD{N40} = tabular_CPD(bnt,N40,[xlsread('FLFDJCPT.xls','sheet1','E2:E65')]);

bnt.CPD{N41} = tabular_CPD(bnt,N41,[xlsread('FLFDJCPT.xls','sheet1','F2:F9')]);

bnt.CPD{N42} = tabular_CPD(bnt,N42,[xlsread('FLFDJCPT.xls','sheet1','G2:G65')]);

bnt.CPD{N43} = tabular_CPD(bnt,N43,[xlsread('FLFDJCPT.xls','sheet1','H2:H9')]);

bnt.CPD{N44} = tabular_CPD(bnt,N44,[xlsread('FLFDJCPT.xls','sheet1','I2:I65')]);

bnt.CPD{N45} = tabular_CPD(bnt,N45,[xlsread('FLFDJCPT.xls','sheet1','J2:J65')]);

bnt.CPD{N46} = tabular_CPD(bnt,N46,[xlsread('FLFDJCPT.xls','sheet1','K2:K65')]);

bnt.CPD{N47} = tabular_CPD(bnt,N47,[xlsread('FLFDJCPT.xls','sheet1','L2:L17')]);

bnt.CPD{N48} = tabular_CPD(bnt,N48,[xlsread('FLFDJCPT.xls','sheet1','M2:M5')]);

bnt.CPD{N49} = tabular_CPD(bnt,N49,[xlsread('FLFDJCPT.xls','sheet1','N2:N17')]);

bnt.CPD{N50} = tabular_CPD(bnt,N50,[xlsread('FLFDJCPT.xls','sheet1','O2:O4097')]);

代碼“discrete_nodes = 1:N;”表示所有節(jié)點均為離散變量；代碼“node_sizes =[x x x x x]”描述各網(wǎng)絡(luò)節(jié)點的狀態(tài)取值數(shù)量，其中節(jié)點44、45、46狀態(tài)數(shù)量為4，其余節(jié)點均為2；代碼“bnt.CPD{Nx} = xxx”用于導(dǎo)入各節(jié)點的概率參數(shù)，其中先驗概率直接鍵入，條件概率參數(shù)較多，通過讀取“FLFDJCPT.xls”文件導(dǎo)入。

3 實驗結(jié)果與分析

3.1 故障診斷結(jié)果

在完成風(fēng)力發(fā)電機貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)和條件概率分布建模后，即可使用FullBNT工具箱進(jìn)行推理。例如需要得出系統(tǒng)故障時各零部件故障的概率，則可將節(jié)點50置為1，其他節(jié)點設(shè)為隱含狀態(tài)，調(diào)用“聯(lián)合樹推理引擎”分別求解各節(jié)點取值為“1”時的概率，Matlab代碼如下：

engine = jtree_inf_engine(bnt);

evidence =cell(1,N);

evidence{N50} = 1;

[engine, loglik] = enter_evidence(engine,evidence);

fori=1:50

marg = marginal_nodes(engine,i);

NBP(i)=marg.T(1)

end

bar(NBP)

代碼執(zhí)行結(jié)果如圖6所示。

圖6 風(fēng)力發(fā)電機系統(tǒng)零部件故障概率分布

從圖6可以看出，當(dāng)系統(tǒng)故障時，除去節(jié)點50本身，節(jié)點36“冷卻風(fēng)扇”故障概率最高，節(jié)點47“集電環(huán)組件”和節(jié)點1“冷卻風(fēng)扇電機”故障概率次之。節(jié)點36“冷卻風(fēng)扇”是故障概率最高的組件，節(jié)點1“冷卻風(fēng)扇電機”是故障概率最高的零部件。

假設(shè)節(jié)點44“繞組溫度”超量程，節(jié)點46“冷卻系統(tǒng)溫度”超量程，求其他節(jié)點故障概率，則將代碼“evidence{N50} = 1;”改為“evidence{N44} = 1; evidence{N46} = 1;”，執(zhí)行結(jié)果如圖7所示。

圖7 繞組及冷卻系統(tǒng)溫度超量程故障概率分布

從圖7可以看出，當(dāng)節(jié)點44 “繞組溫度”和節(jié)點46“冷卻系統(tǒng)溫度”同時出現(xiàn)超量程故障時，除去節(jié)點44和節(jié)點46本身，故障概率最高的節(jié)點依次為節(jié)點26、節(jié)點50和節(jié)點42。節(jié)點26“線纜”是節(jié)點42“發(fā)電機本體”組成零部件，節(jié)點50“風(fēng)力發(fā)電機系統(tǒng)”為系統(tǒng)本身，節(jié)點26故障會導(dǎo)致節(jié)點42組件故障，進(jìn)而導(dǎo)致節(jié)點50系統(tǒng)故障。節(jié)點26“發(fā)電機本體-線纜”是故障概率最高的零部件，故障概率為99.55%，遠(yuǎn)高于其他零部件的故障概率。

3.2 與D矩陣診斷結(jié)果的對比

3.2.1D矩陣診斷示例1

D矩陣故障診斷方法是典型的定性模型方法，其核心是建立故障原因與故障現(xiàn)象的依賴關(guān)系[10]，再通過故障現(xiàn)象推導(dǎo)故障原因。D矩陣依賴模型與貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)類似，同樣利用帶有箭頭的連線來表達(dá)故障原因與故障現(xiàn)象之間的關(guān)聯(lián)。若某故障原因會導(dǎo)致某故障現(xiàn)象出現(xiàn)，則在它們之間通過箭頭進(jìn)行連接。D矩陣故障診斷通常將依賴模型轉(zhuǎn)化到一個布爾矩陣中，利用“1”、“0”來代表“依賴”或“不依賴”關(guān)系。風(fēng)力發(fā)電機系統(tǒng)節(jié)點50與其父節(jié)點的依賴矩陣如表6所示。

表6 節(jié)點50依賴矩陣

表6中第一列為故障原因的節(jié)點編號，第一行為故障現(xiàn)象的節(jié)點編號，單元格中的“1”表示對應(yīng)行的首列的節(jié)點故障會導(dǎo)致對應(yīng)列的首行的節(jié)點出現(xiàn)故障。如節(jié)點41所在行中對應(yīng)列節(jié)點41和節(jié)點50的單元格為“1”，則表示節(jié)點41故障時，節(jié)點41自身和節(jié)點50會出現(xiàn)故障。D矩陣診斷推理基本規(guī)則為：當(dāng)依賴矩陣中某一故障現(xiàn)象存在，則具有依賴關(guān)系的故障原因可能存在；反之，當(dāng)依賴矩陣某一故障現(xiàn)象不存在時，則具有依賴關(guān)系的故障原因肯定不存在。如表6所示第一行的最后一列節(jié)點50出現(xiàn)故障時，則該列單元格為“1”的對應(yīng)行的節(jié)點41、節(jié)點42、節(jié)點44、節(jié)點45、節(jié)點46、節(jié)點47、節(jié)點48、節(jié)點49、節(jié)點50故障現(xiàn)象均可能存在。同理，通過節(jié)點41、節(jié)點42、節(jié)點44、節(jié)點45、節(jié)點46、節(jié)點47、節(jié)點48、節(jié)點49的依賴矩陣，可以推導(dǎo)出他們的父節(jié)點也可能故障。因此，節(jié)點50“風(fēng)力發(fā)電機系統(tǒng)”故障時，采用D矩陣診斷方法推理得到的結(jié)果是：節(jié)點1到節(jié)點49均可能存在故障，即系統(tǒng)中所有零部件都可能故障。

3.2.2D矩陣診斷示例2

對于節(jié)點44 “繞組溫度”和節(jié)點46“冷卻系統(tǒng)溫度”超量程故障，D矩陣診斷過程如下。首先建立節(jié)點44、節(jié)點46及其父節(jié)點的依賴矩陣，如表7所示。

根據(jù)表7所示依賴關(guān)系，當(dāng)節(jié)點44、節(jié)點46故障時，可直接推導(dǎo)出節(jié)點26、27、33、35、36可能故障，再通過節(jié)點36故障可以推導(dǎo)出節(jié)點1、2、3、4、5可能故障。因此，節(jié)點44 “繞組溫度”和節(jié)點46“冷卻系統(tǒng)溫度”超量程故障，D矩陣診斷結(jié)果為：節(jié)點1、2、3、4、5、26、27、33、35、36均有可能故障。

表7 節(jié)點44、46依賴矩陣

表8 診斷結(jié)果對比表

3.2.3 診斷結(jié)果對比

通過上述兩個示例可以看出，D矩陣診斷通常只能得出可能的故障原因的集合，而貝葉斯網(wǎng)絡(luò)診斷不僅可以給出可能的故障原因，還能給出其置信度。貝葉斯網(wǎng)絡(luò)故障診斷結(jié)果與D矩陣故障診斷結(jié)果對比如表8所示。

5 結(jié)束語

貝葉斯網(wǎng)絡(luò)基于概率統(tǒng)計學(xué)理論具有嚴(yán)謹(jǐn)?shù)睦碚摶A(chǔ)，是處理不確定關(guān)系問題的首選模型，適用于復(fù)雜系統(tǒng)的故障診斷及維修決策。通過對風(fēng)力發(fā)電機系統(tǒng)的建模和模擬故障診斷結(jié)果可以看出，相比于傳統(tǒng)相關(guān)性建模故障診斷方法具有更高的故障分辨率，其輸出的故障原因置信度可作為復(fù)雜系統(tǒng)維護(hù)保障的有力參考。貝葉斯網(wǎng)絡(luò)建模較D矩陣建模復(fù)雜，需要大量的統(tǒng)計數(shù)據(jù)支撐，數(shù)據(jù)獲取和處理難度較大，實際應(yīng)用時還要考慮模型修正、概率參數(shù)隨時間的動態(tài)變化等復(fù)雜因素。目前貝葉斯網(wǎng)絡(luò)在裝備故障診斷領(lǐng)域的實際應(yīng)用案例比較少，相關(guān)方法、模型和算法有待優(yōu)化、完善，還需投入大量的研究。