一種增量式貝葉斯算法及篦冷機故障診斷

劉兆倫 張春蘭 武 尤 王海羽 劉 彬

1.燕山大學河北省特種光纖與光纖傳感實驗室，秦皇島，066004 2.燕山大學信息科學與工程學院，秦皇島，066004

0 引言

貝葉斯網絡(Bayesian network，BN)作為不確定性知識表達和推理中的重要理論模型，已在數據挖掘和故障診斷等領域得到廣泛應用[1]。目前故障診斷研究大部分采用批量式算法，該算法對已有數據進行一次性學習得到網絡結構，新數據到來時拋棄已有結構，將新舊數據整合到一起，在一個更大的數據集下重新學習，需要較大的存儲空間和計算時間，造成了資源浪費，考慮到實際應用中很多數據不斷產生，批量式算法已不能滿足要求。而增量式算法依次處理數據樣本，能夠將變化環境中的新采樣納入模式中，且更新在先前基礎上完成，是一種以新數據的出現為順序來更新學習結果的學習方式[2]。

許多學者對批量式算法在故障診斷方面的應用進行了研究。劉彬等[3]提出了一種基于BN的改進結構算法并建立了具有較高診斷準確率的回轉窯故障診斷模型，但該模型沒有考慮新數據產生時模型結構和參數的變化情況，不具有靈活性；ZHU等[4]提出了一種基于BN的電力系統故障診斷方法，該方法建立了基于簡化BN的三元模型，用于輸電系統故障區段的估計，新數據到來時需對所有數據進行重新學習，不能適應外界環境變化，說明批量式算法已經不能滿足實際應用需求。近年來，部分學者對增量式算法進行了研究。田鳳占等[5]將EM(expectation maximization)算法和遺傳算法引入增量學習過程，定義新變異算子且擴展傳統交叉算子，提出了一種包含隱變量的增量學習算法，提高了學習精度，但該算法在存儲空間方面不占優勢；ALCOBE[6]提出了一種將批量爬山(hill-climbing，HC)算法轉換為增量形式的啟發式方法，在一定程度上節省了時間和存儲空間；王飛等[7]提出了一種基于遺傳算法的增量學習方法，需在時效性和準確率之間折中選擇；ZHU等[8]建立了一種基于BN和增量學習的農作物疾病動態診斷的數學模型，但該模型實時性較差；YASIN等[9]對MMHC算法進行改進，得到了一種增量貝葉斯結構學習算法，但該算法不能滿足應用需求。考慮到實際應用，隨著新數據的增加會出現兩種情況：①新舊數據來自同一分布，只需更新參數；②新舊數據來自不同分布，參數偏離穩定，數據與結構表現出不適應，當這種不適應達到一定程度時，則需更新結構。由此增量學習方法可分為兩類：一類是結構不變，只進行參數調整；另一類是先更新結構，后更新參數。

針對上述現狀，本文將專家知識與數據學習相結合得到原結構和原參數，構造WTUN (whether to update network)函數用于判斷新數據集到達時原結構是否改變，若是，則定義Affect函數用于衡量數據對結構中各節點的平均影響程度，利用爬山算子對影響較大節點的馬爾可夫毯結構進行修正得到候選結構，利用改進的評分函數選擇評分最大的結構作為最優結構，無論結構是否改變均利用EM算法更新參數，得到一種貝葉斯增量學習算法(incremental learning algorithm, ILA)。將該算法應用到篦冷機工藝參數建模中，利用模型分析參數間的影響，對原模型進行增量維護，最后對二次風溫進行故障診斷，找出導致故障發生的原因并采取措施，可提升熱回收效率從而達到節約能源的目的。

1 ILA算法

1.1 ILA算法的構建

數據學習要求獲得大量數據樣本，且所得結構不準確，專家知識學習主觀性較強，易產生較大誤差，將兩種方法相結合會相對準確且能提高效率[10]，本文結合數據學習與專家知識獲得原貝葉斯結構G?，利用最大似然估計(maximum likelihood estimation,MLE)算法學習得到原參數θ?ijk?。

在實際應用中，新數據連續產生，每隔Δt?時間采集一次新數據，將進行批量式學習的數據作為初始舊數據集D?o，可表示為D?o=f?(D?p,t?0)；Δt?1時間內產生的數據作為新數據集D?n，而Δt?2時間內新數據到達時將其作為新數據集D?n，Δt?1時間內的數據則作為舊數據，依此類推，學習過程中只保留當前的D?o和D?n，過程如下：

(1)

其中，t?0為初始時間，k?=1,2，…,m?，m?為新數據集的個數。BIC函數如下：

(2)

其中，P?(D?/G?)表示似然概率，即給定結構情況下關于數據的后驗概率；q?i?是變量X?i?父節點取值組合的數目；r?i?是變量X?i?取值的數目；N?ijk?是X?i?取k?值、X?i?的父節點取j?值時的樣本數目，也稱為充分統計因子。由式(2)可以判斷數據與結構的擬合程度，其實質是對每個樣本在結構中有相互影響關系的變量條件概率的求和[11]，依據其思想，構建WTUN函數：

WTUN?(D?|G?)=

(3)

依據馬爾可夫毯原則[12]對式(3)進行化簡。在BN中節點x?i?的馬爾可夫毯結構包括它的父節點、子節點及子節點的父節點，當節點x?i?的馬爾可夫毯給定時，在該網絡中x?i?與其他變量條件獨立，因此，網絡中的條件概率可化簡為

P?(x?i?|x?1,x?2,…，x?i?-1,x?i?+1,…,x?n?)=

P?(x?i?|MB?(x?i?))

(4)

依據乘法公式，可將式(4)變換為

P?(x?i?|MB?(x?i?))=P?(x?i?,MB?(x?i?))/P?(MB?(x?i?))

(5)

由于式(5)中P?(MB?(x?i?))不包含x?i?，即無論x?i?取何值，P?(MB?(x?i?))的取值相同，視為常量，根據隨機變量形式鏈式規則，將式(5)化簡為

(6)

P?(x?i?|MB?(x?i?))=bP?(x?i?|π?(x?i?))·

(7)

將式(7)代入式(6)，再代入式(4)可得

(8)

式中，π?(x?i?)為節點x?i?的父節點；x?j?為貝葉斯結構中節點x?i?的子節點；π?(x?j?)為節點x?i?的子節點的父節點。

由式(8)可看出，WTUN函數是對結構中有相互影響關系變量的條件概率的求和，與BIC函數計算結果相比，只是某些概率項的系數有所增大，故可用于判斷數據與原結構的擬合程度。對于給定的D?o、D?n，以及正數閾值β?∈[0,1]，若

(9)

則需更新結構，否則只更新參數。

若式(9)不成立，即原結構不需要調整，只需對BN進行參數學習，根據EM算法思想[13]，將原參數變為先驗參數，將以下結果作為新的參數：

(10)

若式(9)成立，則需要進行結構調整，在新的結構上更新參數。構造Affect函數，計算當前數據對各節點的平均影響程度，給定閾值δ?，存儲并記錄Affect值大于閾值δ?的節點集S?，通過爬山算子S?中各節點的馬爾可夫范圍進行修正，得到候選結構集G?q，利用改進評分函數對G?q進行評分，選取評分最高的結構作為最優結構G?good，并利用式(10)更新參數，完成BN的增量維護。

已知P?(D?n|G?)可以體現新數據與原結構的擬合程度，新數據對G?的所有影響由P?(D?n|G?)的變化相對于參數θ?ijk?的變化組成[14]，則當前數據對原結構中各節點的平均影響表示為

(11)

式中，A?為新數據對原結構每個節點的平均影響；B?為原數據對原結構每個節點的平均影響。

對A?進行化簡，有

(12)

其中，d?n為新數據集樣本個數；P?為結構對應的條件概率集，又已知

則式(12)可化簡為

(13)

同理可得B?，則Affect函數表達式為

(14)

θ?ijk?=P?(X?i?=k?|π?(X?i?=j?))

利用式(14)計算每個節點的Affect值，給定閾值δ?，存儲并記錄Affect值大于閾值δ?的節點集S?，利用增加弧、反轉弧、刪除弧三種操作集合op?={op?(1),op?(2),…,op?(k?)}修正S?中每個節點的馬爾可夫毯結構，得到候選結構集：

G?q?=op?(j?)(MB?(S?i?))

(15)

式中，op?(j?)為3種操作之一；MB?(S?i?)為節點集S?中某個節點的馬爾可夫毯結構。

對G?q中的結構進行評分：

(16)

式中，N?n、N?o分別為新舊數據的數量；λ?為調整新舊數據學習速度及趨勢的因素，如果新數據與原結構不匹配，則λ?增大，否則λ?減小；g?為結構集G?q中的個體；pen?(g)為結構g?的懲罰函數，防止數據與結構的過度擬合。

選取G?q中評分最大的結構作為最優增量結構G?good，即

(17)

得到G?good后，利用式(10)更新參數，完成BN的增量學習過程。

1.2 ILA算法實現

ILA算法的流程圖見圖1。對ILA算法分步進行簡要描述，算法步驟如下：

圖1 ILA算法流程圖Fig.1 The flow chart of ILA algorithm

(1)輸入初始數據集、批量式算法與專家知識結合獲得原貝葉斯結構G?，利用MLE算法學習得到原參數；

(2)檢測到新數據，利用WTUN函數判斷G?是否改變，若是則轉步驟(3)，否則利用EM算法更新參數，轉步驟(6)；

(3)利用Affect函數找到數據對節點影響大的節點集S?，并得到每個節點的馬爾可夫毯結構的集合G?p；

(4)在原結構中利用op?(j?)對步驟(3)得到的G?p中的結構進行修改，得到候選結構集G?q；

(5)使用Score函數對步驟(4)中的每個候選結構進行評分，以評分最大的結構作為增量結構G?good，利用EM算法更新參數；

(6)輸出增量學習后的結構和參數。

2 ILA算法的仿真實驗

2.1 WTUN函數有效性驗證

有新數據的WTUN函數值越大，說明貝葉斯網絡對新數據的擬合程度越好，隨著新數據的增加，與新舊結構的WTUN值在不斷變化，若與新結構的相適性越來越好，則可將WTUN函數用于判斷數據與結構擬合程度的可行性。

利用經典Asia、Car及Alarm網絡進行仿真驗證，以Asia網絡為例，其標準網絡如圖2所示。利用Netica軟件的smulate cases功能生成8000組數據作為原數據集D?o，修改概率分布表，在不產生環路的條件下，刪除和增加一定數量的邊后得到的網絡如圖3所示，用同樣的方法再產生8000組數據作為新數據集D?n，部分數據結果見表1。

圖2 原Asia網絡Fig.2 The original Asia network

圖3 新Asia網絡Fig.3 The new Asia network

表1 Asia網絡部分數據舉例

對圖2所示的網絡，分別加入0～100%的新數據，隨著新數據的增加，分別計算含不同比例新數據下的WTUN?(D?n|G?)及WTUN?(D?n|G?n)，對Car、Alarm與Asia進行相同操作，WTUN值的變化情況分別如圖4和圖5所示。

圖4 WTUN值隨原結構變化情況Fig.4 The WTUN value varies with the original structure

圖5 WTUN值隨新結構變化情況Fig.5 The WTUN value varies with the new structure

由圖4可知，隨新數據的增加，與原結構的WTUN值越來越小，即數據與原結構的適應性越來越差。由圖5可知，與新結構的WTUN值越來越大，即數據與新結構的擬合程度越來越好，而且變化比較平穩，從而驗證了WTUN函數的有效性。

2.2 Affect函數正確性驗證

利用Asia網絡驗證Affect函數的正確性。同樣利用圖2所示網絡作為原網絡，生成2 000組原數據，利用圖3所示網絡作為新網絡，生成8 000組新數據，即考慮10 000組數據中含20%原數據和80%新數據的情況，說明α?的值為0.8。利用式(14)計算數據對原網絡各個節點的Affect值，并給定閾值δ?為1，以節點Bronchitis為例說明式(14)的具體計算過程，該節點的原參數與新參數的學習結果分別見表2和表3。

當Bronchitis=present，Smoking=smoker時，節點Bronchitis的Affect值為

表2 節點Bronchitis的原條件概率

表3 節點Bronchitis的新條件概率

其中，497和156分別為D?o與D?n中滿足Bronchitis=present、Smoking=smoker的樣本數量。同理可計算另外3種情況下的數值，這4個值求和的結果即節點Bronchitis的Affect值。因此，可得到每個節點的Affect值，結果見表4。

表4 Asia網絡中數據對各節點的影響情況

由表4結果可知，Affect函數判斷數據對結構影響度的準確率為75%。為了排除數據的偶然性，在上述條件下進行20次試驗取平均值，最終得到準確率為72.5%，由此可以驗證Affect函數的正確性。

2.3 增量學習算法有效性仿真證明

在Affect函數正確性驗證實驗中，由表4可知，當加入8 000組新數據時，根據Affect值大小判斷需要修正的節點分別為D?、S?、B?、v?、L?、T?，對這些節點馬爾可夫毯結構用爬山算子進行修改，利用Score函數進行評分，得到評分最大的最優網絡結構，如圖6所示。

圖6 ILA算法得到的最優結構Fig.6 Optimal structure obtained by incremental algorithm

考慮評分比值、Log-less比值、時間比值3個指標，將ILA算法分別與批量爬山、增量爬山(incremental hill-climbing，IHC)算法[6]進行學習結果質量以及時間消耗兩方面的仿真對比分析。通過Nursery網絡，Alarm網絡生成的數據集對算法進行仿真，對每組數據集獨立運行算法 20 次，運行結果見表5和表6。評分比值、Log-less比值、時間比值定義分別為

(18)

(19)

(20)

表5 ILA與HC算法學習質量與時間消耗對比

表6 ILA與IHC算法學習質量與時間消耗對比

由表5和表6可看出ILA算法與批量HC算法及增量IHC算法的評分比值和Log-less比值很接近，說明本文所提算法與HC算法、IHC算法學習到的結構質量相當。時間比值與1相差越大，則ILA算法的時間越占優勢，表5表明本文所提算法相比批量式算法明顯減少了時間消耗，且數據量越大，本文算法的時間消耗越有優勢。表6表明ILA算法與IHC算法的時間消耗相當，個別情況稍遜于IHC算法。

根據Alarm網絡生成的訓練數據，分別讀入k?為500、1 000個事例進行一次增量學習，分別比較ILA算法與HC算法、增量遺傳算法(incremental genetic algorithm，IGA)[7]的存儲量，如圖7所示。

圖7 三種算法的存儲量對比Fig.7 The comparison of the storage capacity ofthree algorithms

由圖7可以看出，批量HC算法的存儲量隨事例數的增大呈線性增長，這是因為批量HC算法每次要將新舊數據合在一起重新學習，IGA 需要的存儲量先上升達到一個頂峰后下降的過程，雖然后期存儲量趨于常數，但前期存儲量較大，而本文所提算法通過Affect函數判斷出新數據對原結構哪些節點影響大，只對這些節點進行了局部修改，極其有效地提高了效率，減少了存儲空間，存儲開銷優于HC與IGA算法。綜上所述，ILA算法是有效的。

3 篦冷機工藝參數建模及故障診斷

篦冷機是新型干法水泥生產線的關鍵設備，擔負著出窯熟料的冷卻、輸送和熱回收任務[15]，其工作狀態直接影響到水泥廠的效益及水泥用戶的利益，但其工藝參數較多且影響錯綜復雜，且故障具有很大不確定性，人工排查方式效率很低。由于BN強大的推理能力及方便的決策機制，故將ILA算法應用到篦冷機熟料換熱系統，對診斷模型進行維護。首先根據篦冷機原理確定BN的變量，其次結合專家知識與數據學習建立篦冷機原始模型，然后利用ILA算法的增量機制更新初始模型，最后利用變量消元算法進行故障診斷。

3.1 變量選取與工藝參數建模

篦冷機結構如圖8所示，根據其工藝原理選取變量，其范圍與狀態分類的對應關系見表7。利用MATLAB對這8個參數的數據進行離散化處理，然后通過專家知識與SHC算法[16]結合的方法建立原始模型，利用EM算法對所建立的結構和原數據進行參數學習，得到原始模型，如圖9所示。

圖8 篦冷機結構示意圖Fig.8 Schematic diagram of grate cooler

表7 變量范圍與變量狀態分類的對應關系

圖9 篦冷機原始模型Fig.9 Original model of grate cooler

3.2 基于ILA算法的診斷模型更新

在原始模型的基礎上，檢測到新數據時，利用ILA算法，加入1 000組新數據和加入10 000組新數據時的學習結果分別如圖10、圖11所示，實現了診斷模型的更新。

圖10 加入1 000組新數據時增量學習結果Fig.10 Incremental learning results when adding1 000 new sets of data

圖11 加入10 000組新數據時增量學習結果Fig.11 Incremental learning results when adding10 000 new sets of data

由圖10和圖11可以看出，1 000組新數據不足以引起結構的改變，只是參數進行了更新，而加入10 000組新數據則引起了結構的變化，同時參數也發生改變。根據初始模型，利用ILA、IHC算法進行增量學習，利用SHC算法進行批量學習得到的診斷模型如圖12～圖14所示。由圖12～圖14可知：利用SHC得到的診斷模型偏差較大，如忽略了喂料量對二次風溫的影響，錯誤添加了二室風溫對二室篦下壓力的影響。

圖12 利用ILA算法得到的診斷模型Fig.12 Diagnosis model based on ILA algorithm

圖13 利用IHC算法得到的診斷模型Fig.13 Diagnosis model based on IHC algorithm

圖14 利用SHC算法得到的診斷模型Fig.14 Diagnosis model based on SHC algorithm

圖15 三種算法得到篦冷機結構的時間消耗Fig.15 Comparison of the three algorithms for the timeconsumption of the structure of the grate cooler

圖16 三種算法得到篦冷機結構的存儲量Fig.16 Comparison of the three algorithms for thestoragecapacity of the structure of the grate cooler

3種算法得到的模型的運行時間以及存儲開銷對比如圖15、圖16所示。由圖15、圖16可知，在篦冷機熟料換熱系統中，隨著數據量的增大，批量式SHC算法的運行時間越來越長，存儲量呈線性增長，而增量ILA算法時間消耗增加不明顯，稍遜于增量IHC算法，雖然開始存儲量稍遜于IHC算法，但當觀測數據量達到8 000左右時，與IHC算法相當，基本可以滿足要求。

3.3 熟料二次風溫故障診斷

二次風溫是篦冷機的重要參數，由篦下高壓冷卻風機與高溫熟料進行熱交換得到，二次風在窯頭負壓作用下，進入回轉窯供窯內燃燒，二次風溫偏低會使熱回收效率變差，造成燃料煤的大量浪費，同時產生大量廢氣污染環境，提高二次風溫可、節省燃煤，達到節約能源的目的。如監測到二次風溫忽然偏低，由圖12可知導致二次風溫忽然偏低的原因為V?、R?、M?及G?p，以10 000組數據為例，利用變量消元法計算后驗概率：

α?(0.0648 0.1174 0.8178)

α?(0.3524 0.2191 0.4285)

α?0.2416 0.6545 0.1039

α?0.3248 0.2159 0.4393

通過比較后驗概率α?可知，二次風溫忽然降低是由篦速偏高導致的，應及時減小篦速，最大程度上減少燃煤的浪費與氮氧化物等污染廢氣的排放。從數據集中篩選出滿足二次風溫偏小征兆的數據1 000、3 000、5 000、8 000組進行測試性診斷實驗，記錄每一組數據對應的診斷結果，部分驗證結果見表8。

表8 ILA算法的部分驗證結果

由診斷結果可得到各算法在每一組測試集的正確診斷故障數：

H?=D?ture/D?test

(21)

其中，H?為正確診斷率；D?ture為診斷正確的數據數；D?test為測試集大小。最終得到3種算法對比結果，如圖17所示。由圖17可知：在每組測試數據中，利用ILA對二次風溫進行故障診斷的準確率要高于SHC算法以及IHC算法，能夠基本滿足實際生產中二次風溫故障診斷的需求。

圖17 三種算法診斷準確率對比Fig.17 Comparison of diagnostic accuracy of three algorithms

4 結論

仿真實驗結果表明，ILA算法能有效完成網絡更新，一定程度上節省了時間和存儲空間，在真實數據環境下建立的篦冷機診斷模型具有較高的診斷準確率，避免了由于二次風溫故障診斷不及時導致的浪費資源、污染環境等問題，為篦冷機換熱系統故障診斷提供了新的思路。