李 榕, 李 元
(沈陽化工大學(xué) 信息工程學(xué)院, 遼寧 沈陽 110142)
隨著工業(yè)生產(chǎn)過程越來越復(fù)雜多變,生產(chǎn)過程中對于設(shè)備的安全性和可靠性要求也愈加嚴(yán)格[1-2]。為了更多地捕捉過程信息,往往安裝多樣的傳感器來獲取大量反映過程信息的數(shù)據(jù)。通過分析正常數(shù)據(jù)和不同故障運行模式下的數(shù)據(jù),實現(xiàn)基于數(shù)據(jù)驅(qū)動的故障診斷[3-5]。
工業(yè)過程數(shù)據(jù)往往具有維數(shù)高、復(fù)雜性強等特點,若直接對獲取的原始數(shù)據(jù)進行故障診斷,不僅消耗計算機資源多,而且診斷效果較差,因此,在工業(yè)過程數(shù)據(jù)上,對原始數(shù)據(jù)首先進行降維處理,去除數(shù)據(jù)冗余信息、壓縮數(shù)據(jù)維度從而減少存儲空間,消除數(shù)據(jù)噪聲進而提高分類器的分類準(zhǔn)確度。主成分分析(principal component analysis, PCA)[6]作為工業(yè)過程數(shù)據(jù)中最經(jīng)典的降維方法,能夠處理高維、嘈雜和高度相關(guān)的數(shù)據(jù),因其優(yōu)異的表現(xiàn)而被廣泛應(yīng)用于故障診斷領(lǐng)域。陳輝等[7]為了解決在設(shè)備故障分類識別中易受輸入樣本相關(guān)性影響的問題,提出PCA-BP神經(jīng)網(wǎng)絡(luò)方法對設(shè)備進行故障預(yù)測。張弛等[8]提出主成分分析和優(yōu)化參數(shù)支持向量機的智能變電站故障診斷模型,用以解決目前智能變電站故障診斷結(jié)構(gòu)復(fù)雜、樣本數(shù)據(jù)量小的問題。但是,當(dāng)數(shù)據(jù)分布具有非高斯和非線性特征時,由于使用二階統(tǒng)計量和線性假設(shè),PCA此時表現(xiàn)不佳。核主成分分析(kernel principal component analysis, KPCA)[9-11]在一定程度拓寬了PCA的應(yīng)用范疇,利用非線性映射將輸入數(shù)據(jù)映射到高維特征空間,在特征空間執(zhí)行PCA。李夢瑤等[12]用KPCA降低數(shù)據(jù)維度、剔除冗余信息,從而實現(xiàn)齒輪箱的故障診斷。武文棟等[13]提出KPCA結(jié)合麻雀搜索算法優(yōu)化核極限學(xué)習(xí)機的方法,用以提高光伏陣列故障診斷的精度。然而,在實際工業(yè)過程中,數(shù)據(jù)通常同時存在著高斯分布和非高斯分布特性,而KPCA在選取主成分的過程中只使用過程變量的方差信息,因此缺少處理非高斯數(shù)據(jù)的能力,對一些復(fù)雜的系統(tǒng)故障處理效果也就不理想。
2010年,Jenssen[14]提出一種新的降維方法——核熵成分分析方法(kernel entropy component analysis, KECA)。對比KPCA只考慮方差信息只能提取數(shù)據(jù)的高斯特征。KECA[15-17]以信息熵為信息衡量的指標(biāo),使用高階統(tǒng)計量最大程度保留原始數(shù)據(jù)的熵值,同時,其能夠捕捉到數(shù)據(jù)內(nèi)部的結(jié)構(gòu)信息。如今,KECA在干旱監(jiān)測、柴油機故障識別、風(fēng)電功率預(yù)測等領(lǐng)域被廣泛應(yīng)用[18-20]。但目前未查到該方法在化工工業(yè)過程多類型故障診斷領(lǐng)域的應(yīng)用。
針對傳統(tǒng)算法在多類型故障診斷率較低的問題,本文提出一種KECA-CSK的故障診斷方法。KECA進行多種數(shù)據(jù)分析時,不同類別的數(shù)據(jù)樣本往往互相垂直,而同一類別的數(shù)據(jù)樣本往往共線,基于此特點設(shè)計一種基于余弦相似度的K均值聚類器,建立診斷模型,從而對未知類型的故障數(shù)據(jù)進行判別。本文對田納西-伊斯曼化工過程(Tennessee-Eastman,TE)進行模擬實驗,驗證KECA-CSK均值方法在故障診斷領(lǐng)域的有效性。
定義數(shù)據(jù)集D={x1,x1,…,xn}的概率密度函數(shù)p(x),其Renyi熵H(p)表示為
(1)
由式(1)可以看出,Renyi熵估計取決于數(shù)據(jù)集的概率密度函數(shù),不需要對數(shù)據(jù)的分布進行限定,算法的適用性更強。因為對數(shù)函數(shù)本身具有單調(diào)性,因此只需考慮
(2)
(3)
(4)
式中:K為n×n的核矩陣;I為n×1的單位向量。核矩陣包含了二次Renyi熵估計值所需的所有元素。對核矩陣進行特征分解,
K=EΛET。
(5)
式中:Λ=diag(λ1,λ2,…,λn)表示特征值矩陣;E=(e1,e2,…,en)表示特征向量矩陣。將式(5)代入式(4),得到
(6)
由式(6)可得,Renyi熵可以表示為n個分量的累積,各個分量都同時包含特征值與特征向量。
將輸入數(shù)據(jù)集D={x1,x2,…,xn},通過非線性φ映射到核特征空間定義為xi→φ(xi)(i=1,2,…,n),該特征空間的數(shù)據(jù)可表示為Φ=[φ(x1),φ(x2),…,φ(xn)]。由于KECA可以被看作最大限度地保留Renyi熵估計的主軸在包含原始數(shù)據(jù)最多信息的情況下構(gòu)造成的子空間,即將數(shù)據(jù)映射到由k個KPCA主軸張成的子空間Uk上,對ζi值從大到小排序,選取前k個特征值和特征向量,轉(zhuǎn)換后的數(shù)據(jù)表示為
(7)
對于新的樣本xnew在Uk上的投影可表示為
(8)
K均值算法由于其原理簡單、易于理解、容易實現(xiàn)以及理論可靠等優(yōu)點而被廣泛應(yīng)用于不同領(lǐng)域。由于TE數(shù)據(jù)經(jīng)KECA特征提取后不同類的數(shù)據(jù)分布在不同的角度方向上,因此本文使用樣本到各個質(zhì)心的夾角余弦值來衡量數(shù)據(jù)樣本之間的相似性,即余弦值越大,表明數(shù)據(jù)與該類別的數(shù)據(jù)樣本相似度越高。
假設(shè)兩個包含m維特征的向量,表示為α=(a1,a2,…,am),b=(b1,b2,…,bm),則夾角余弦值可表示為
(9)
當(dāng)α,b方向一致時,數(shù)據(jù)之間夾角越小,夾角余弦值就越接近于1,數(shù)據(jù)相似度也就越高;同理,當(dāng)兩者之間夾角越大,夾角余弦值就越小,數(shù)據(jù)相似度也就越低。
基于CSK的方法步驟如下:
1) 初始質(zhì)心的確定。選取余弦相似度最低的數(shù)據(jù)樣本作為前2個初始質(zhì)心,然后再選擇與前2個初始質(zhì)心余弦相似度最低的數(shù)據(jù)樣本作為第3個初始質(zhì)心,以此類推,直至確定C個初始質(zhì)心;
2) 將數(shù)據(jù)樣本劃分到與C個初始質(zhì)心余弦相似度最高的類別中;
3) 更新每個數(shù)據(jù)類別中的類質(zhì)心;
4) 重復(fù)步驟(2)~(3), 當(dāng)質(zhì)心位置變化小于給定的閾值,循環(huán)結(jié)束,保留各類別中心。
CSK均值方法運用于故障診斷領(lǐng)域,首先是將混合的不同故障類別數(shù)據(jù)進行聚類,聚類結(jié)束以后,保留不同故障類別的聚類中心,構(gòu)建診斷模型。對新獲得的未知類別的測試數(shù)據(jù)進行故障診斷時,確定此故障數(shù)據(jù)樣本與哪一類別的聚類中心余弦值最大,那么該樣本的故障類別屬性就與那一類故障樣本的類別屬性一致。
本文KECA中核函數(shù)選取高斯核函數(shù),表達式為
k(x,xt)=exp(-0.5||x-xt||2/δ2)。
(10)
式中,δ為待優(yōu)化的核參數(shù)。
在聚類分析中,盡可能使得同類數(shù)據(jù)之間相互靠近,不同類數(shù)據(jù)之間相互遠離。因此,核參數(shù)選取規(guī)則定義如下:
1) 類內(nèi)離散度。
(11)
2) 類間離散度。
(12)
式中,C為數(shù)據(jù)樣本類別個數(shù)。
使用基于角結(jié)構(gòu)的類內(nèi)與類間離散度的差作為選取的準(zhǔn)則函數(shù),即:
(13)
根據(jù)Shi等[21]提出的參數(shù)值選取方法,將原始空間樣本歐氏距離中值的10%~20%作為參數(shù)選取范圍,然后利用所定義的準(zhǔn)則函數(shù)選取合適的核參數(shù)。
基于KECA-CSK均值故障診斷方法分為離線建模和在線診斷2部分,KECA-CSK均值方法的診斷流程圖如圖1所示。
1) 獲取訓(xùn)練數(shù)據(jù)共m個變量n個樣本X=[x1,x2,…,xn]∈Rm×n,并進行標(biāo)準(zhǔn)化處理;
2) 構(gòu)建KECA模型,根據(jù)式(7)計算得分矩陣Feca;
3) 計算Feca向量之間的余弦值,確定C個初始聚類中心;
4) 根據(jù)余弦相似度大小,將數(shù)據(jù)劃分到相對應(yīng)的類群中;
5) 更新C個類群的聚類中心;
6) 重復(fù)步驟4)~5),各類群中心變化小于閾值,則停止,獲得聚類中心。
1) 測試數(shù)據(jù)xnew并進行標(biāo)準(zhǔn)化處理;
2) 測試數(shù)據(jù)向KECA模型進行投影,通過式(8)計算得分矩陣Fknew;
3) 計算Fknew與建模時各類群聚類中心的余弦值,將其劃為余弦值最大的類群。
TE過程是由美國Eastman化學(xué)公司的Downs和Vogel共同開發(fā)的TE Benchmark仿真模擬平臺,由該平臺產(chǎn)生TE數(shù)據(jù),近年來被廣泛應(yīng)用于工業(yè)生產(chǎn)故障檢測與診斷[22-26]。TE工業(yè)流程主要包括反應(yīng)器、冷凝器、循環(huán)壓縮機、氣液分離器和汽提塔等5個操作單元,整個過程中共涉及8種物料成分,分別為主要參加反應(yīng)的氣體進料U、C、D、E;惰性不可溶進料B;反應(yīng)副產(chǎn)品F以及反應(yīng)液態(tài)主產(chǎn)物G和H。該平臺總共預(yù)設(shè)了21種故障,其中,IDV1~IDV7為階躍故障,IDV8~IDV12為隨機變化故障,IDV13為慢偏移故障,IDV14~IDV15為堵塞干擾故障,IDV16~IDV21為未知型故障。
本文選取TE過程故障1、故障5和故障18作為研究對象。故障具體描述如表1所示。在訓(xùn)練集中,正常數(shù)據(jù)樣本是在25 h 運行仿真下獲得,樣本總數(shù)為500;故障數(shù)據(jù)樣本是在24 h運行仿真下獲得,樣本總數(shù)為480。本文選取訓(xùn)練集3×400×52,其中3是3類故障類型,400是數(shù)據(jù)樣本個數(shù),52是變量數(shù)。在測試集中,正常數(shù)據(jù)樣本是在48 h運行仿真下獲得,采樣間隔時間為3 min,樣本總數(shù)為960;故障數(shù)據(jù)樣本是在48 h運行仿真下獲得,故障在8 h的時候引入,樣本總數(shù)為960,其中,前160個數(shù)據(jù)樣本為正常數(shù)據(jù)樣本。本文選取測試集為3×400×52。
表1 故障數(shù)據(jù)描述Table 1 Fault data description
3.2.1 核參數(shù)確定
以TE過程為背景進行仿真實驗,訓(xùn)練數(shù)據(jù)選取故障類型1、5和18共1 200個數(shù)據(jù)樣本,測試數(shù)據(jù)為同樣的未知標(biāo)簽的3類數(shù)據(jù)共1 200個數(shù)據(jù)樣本。經(jīng)計算,原始空間樣本歐氏距離中值為50.015,核參數(shù)的選擇區(qū)間為δ=50.015×u,u=0.1∶0.01∶0.2,利用準(zhǔn)則函數(shù)選取合適的核參數(shù),核參數(shù)的取值與準(zhǔn)則函數(shù)值的關(guān)系如圖2所示。為了便于可視化以及診斷分析,同時確定k、l為聚類數(shù)3。
圖2 不同核參數(shù)下的準(zhǔn)則函數(shù)值Fig.2 Criterion function values under different kernel parameters
觀察圖2可得,當(dāng)核參數(shù)取歐氏距離中值的17%時,準(zhǔn)則函數(shù)最大,因此本實驗中核參數(shù)取值為8.502 55。
3.2.2 實驗分析
為了驗證KPCA(KPCA中核參數(shù)的取值方式與KECA取值方式一致)與KECA之間的差異性,將所建立的核矩陣分解后,得到特征值和Renyi熵值,并將其歸一化,如圖3(a)所示,KPCA依據(jù)特征值大小選取第1、2、3主成分,KECA中基于Renyi熵值貢獻選取第1、2、6主成分,即特征值大的不一定取得Renyi熵值貢獻也大。根據(jù)KPCA與KECA兩種方法選擇各自不同的主成分后將數(shù)據(jù)進行投影,如圖3(b)和3(c),KPCA投影將3類數(shù)據(jù)分開后,數(shù)據(jù)分布較為分散,也沒有出現(xiàn)明顯的角結(jié)構(gòu)特性,而KECA分析中,數(shù)據(jù)之間具有明顯的角結(jié)構(gòu)特征,數(shù)據(jù)分布較為緊密,有利于CSK方法對數(shù)據(jù)進行聚類。
圖3 3種故障聚類實驗分析Fig.3 Analysis of three fault clustering experiments
對核矩陣的第1、2、3和6特征向量進行分析,如圖4所示。
圖4 訓(xùn)練數(shù)據(jù)核主成分Fig.4 Kernel principal components of training data
觀察圖4可得,KECA根據(jù)Renyi熵貢獻大小,選擇第e1、e2、e6特征向量。e1、e2、e6特征向量中均有一段對應(yīng)于其中一個數(shù)據(jù)類別的正值或負(fù)值,而對于其他類別值均接近于零,即KECA所選取的特征向量攜帶數(shù)據(jù)集的簇結(jié)構(gòu)信息,因此,轉(zhuǎn)換數(shù)據(jù)集的角結(jié)構(gòu)也將會出現(xiàn),如圖4(c)所示。對于e3特征向量,其在零附近上下波動,對于Renyi熵貢獻也基本為0,KECA按照Renyi熵值大小也就不會選取該特征向量作為投影方向。KPCA對數(shù)據(jù)特征提取時,根據(jù)特征值降序排列選取e1、e2、e3特征向量。然而,e3特征向量分布沒有刻畫出故障18的數(shù)據(jù)的特點,不能獲取數(shù)據(jù)的結(jié)構(gòu)信息。因此,對比KPCA方法,KECA能夠提取不同類數(shù)據(jù)的差異性,獲取數(shù)據(jù)內(nèi)部存在的結(jié)構(gòu)信息,使得轉(zhuǎn)換后的數(shù)據(jù)保持較好的可區(qū)分性。
3.2.3 訓(xùn)練數(shù)據(jù)聚類
基于上述分析后,使用本文所提出的CSK聚類器,對訓(xùn)練數(shù)據(jù)進行聚類,獲取各類群的聚類中心,聚類結(jié)果如圖5所示。
圖5 聚類結(jié)果Fig.5 Clustering result plot
圖5(a)是訓(xùn)練數(shù)據(jù)KECA-CSK聚類結(jié)果,3種類型樣本的聚類正確率分別為95.5%,96.0%和35.0%。圖5(b)是KPCA的聚類結(jié)果,3種故障類型數(shù)據(jù)樣本的聚類正確率分別為89.25%,54.25%和28.00%。KPCA聚類效果相比KECA聚類效果較差。訓(xùn)練數(shù)據(jù)聚類結(jié)果的正確率直接影響對未知類別故障數(shù)據(jù)的診斷正確率。
3.2.4 未知類別的故障數(shù)據(jù)診斷
對于包含與建模數(shù)據(jù)相同的測試樣本,將其投影到建模模型獲得得分向量,計算得分向量與建模結(jié)果中各類群中心的余弦值,將其劃為余弦值最大的那一類別中。若測試數(shù)據(jù)的故障類型不同于建模數(shù)據(jù)的故障類型時,將其代入其他數(shù)據(jù)建立的診斷模型中,圖6是KECA與KPCA的診斷結(jié)果。
圖6 故障診斷結(jié)果Fig.6 Troubleshooting results
圖6(a)是KECA的診斷結(jié)果, 3種未知故障類型的測試數(shù)據(jù)樣本基本均能得到有效判別, 總數(shù)據(jù)樣本的診斷正確率為80.08%。 KECA特征提取時候, 捕捉到的數(shù)據(jù)集群結(jié)構(gòu)信息使得CSK方法在訓(xùn)練數(shù)據(jù)上能夠取得較好的聚類結(jié)果, 進而影響對測試數(shù)據(jù)的診斷結(jié)果。 圖6(b)是KPCA診斷結(jié)果, 診斷結(jié)果較差, 總數(shù)據(jù)樣本的診斷正確率為59.08%, 驗證了KECA在工業(yè)過程故障診斷領(lǐng)域的優(yōu)越性。
針對實際工業(yè)過程中,傳統(tǒng)核主成分多故障類型診斷率低問題,本文提出一種基于KECA-CSK均值的故障診斷方法,并在TE化工過程數(shù)據(jù)取得較好的診斷效果。KECA對數(shù)據(jù)進行特征提取,轉(zhuǎn)換后的數(shù)據(jù)不同類別數(shù)據(jù)之間保持角結(jié)構(gòu)特性,在此基礎(chǔ)上設(shè)計一種基于余弦相似度的K均值聚類方法,從而獲得較高準(zhǔn)確率的建模結(jié)果,提高對未知類型數(shù)據(jù)判別的正確率。 后續(xù)工作集中于探究該方法在其他工業(yè)過程的應(yīng)用是否也能取得較好的診斷結(jié)果。