郭昕剛,霍金花,程 超,許連杰
(長春工業(yè)大學(xué) 計算機(jī)科學(xué)與工程學(xué)院, 吉林 長春 130012)
隨著現(xiàn)代工業(yè)的發(fā)展,數(shù)據(jù)驅(qū)動方法在過程監(jiān)控中發(fā)揮著重要作用[1-2]。在主流方法中,多元統(tǒng)計過程監(jiān)控(multivariate statistical process monitoring,MSMP)被用于監(jiān)控多元復(fù)雜工業(yè)過程中的故障[3]。傳統(tǒng)的多元統(tǒng)計監(jiān)控主要包括主成分分析法(principal component analysis,PCA)、偏最小二乘法(partial least square,PLS)和獨立成分分析[4](independent component analysis,ICA)。
工業(yè)過程的復(fù)雜性和未知的動態(tài)特性使得靜態(tài)過程的監(jiān)測效果較差。傳統(tǒng)的MSMP方法是靜態(tài)過程監(jiān)控方法,即當(dāng)前時間的樣本與過去的樣本沒有關(guān)聯(lián),樣本數(shù)據(jù)相互獨立,忽略了動態(tài)特性。Li等提出了一種部分動態(tài)主成分分析(dynamic principal component analysis,DPCA)來提高動態(tài)過程監(jiān)測能力[5]。針對過程動力學(xué)和數(shù)據(jù)非高斯統(tǒng)計的特點,動態(tài)獨立成分分析(dynamic independent component analysis,DICA)也相繼被提出[6-7]。上述方法以消除動力學(xué)為目標(biāo),往往假設(shè)所有變量具有相同的動態(tài)特性,然后對原始數(shù)據(jù)進(jìn)行擴(kuò)展,與傳統(tǒng)方法相比提高了性能,但產(chǎn)生了大量的冗余信息,處理動態(tài)特性較強(qiáng)的變量時性能較差。慢特征分析[8](slow feature analysis, SFA)可以從時間序列中提取緩慢變化的特征,表征變量變化的快慢程度,是一種有效的無監(jiān)督算法。Shang等提出了基于SFA的動態(tài)監(jiān)控,實現(xiàn)了運(yùn)行和控制性能的監(jiān)控[9]。上述方法雖然克服了動態(tài)缺陷,但只建立單一模型卻忽略了局部信息,導(dǎo)致大規(guī)模工業(yè)過程的診斷性能較差。
多模塊算法最早由Macgregor等提出,有效利用局部信息改進(jìn)過程分析將整個模型分為多個子模型[10]。Ge等根據(jù)主成分分析的不同方向?qū)υ紨?shù)據(jù)進(jìn)行劃分,將線性變量分配到同一個塊中,分塊過程中存在缺失和重疊的問題[11]。Tong等分別分析了變量與主元子空間和殘差空間的相關(guān)性,再將變量分配到相應(yīng)的子空間中實現(xiàn)分塊[12]。這些方法都以先驗知識作為前提,極大程度限制了該方法的應(yīng)用。KL(Kullback-Leibler)散度作為一種新的概率測度,可以測量兩個統(tǒng)計變量之間的差值,用來解決故障診斷問題[13-14]。Wang等利用KL散度將具有類似統(tǒng)計特征的變量劃分為一塊,在每個低維子空間中建立PCA模型,使用貝葉斯策略進(jìn)行融合[15]。周偉等在DPCA建模的基礎(chǔ)上,利用 KL 散度量化模型得到分向量概率分布之間的相似度,從而建立多塊模型實現(xiàn)對微小故障的診斷[16]。
在實際的工業(yè)生產(chǎn)過程中,大多數(shù)緩慢變化的過程數(shù)據(jù)也具有非線性時變的特性,其使正常數(shù)據(jù)出現(xiàn)偏移,導(dǎo)致出現(xiàn)誤警現(xiàn)象。柯亮等提出了基于滑動窗PCA的微小故障檢測方法,利用滑動窗口的策略對數(shù)據(jù)實時更新,使得故障數(shù)據(jù)與非故障數(shù)據(jù)分化明顯,實現(xiàn)了對微小故障的放大,提高了模型的自適應(yīng)能力[17]。
綜合以上問題,提出一種基于KL散度的多模塊滑動窗口慢特征分析故障診斷方法,利用KL散度的統(tǒng)計特性,建立多塊模型,在每個塊中應(yīng)用SFA來提取過程數(shù)據(jù)的不同動態(tài),利用滑動窗口得到最優(yōu)模型,克服依據(jù)先驗知識分塊的策略、單一模型不穩(wěn)定等問題。在田納西伊斯曼(Tennessee Eastman, TE)過程監(jiān)控中取得了較好的診斷效果。
(1)
由于KL散度的定義具有不對稱性,不能作用于距離的度量。 在實際的應(yīng)用中,人們將其改進(jìn)為對稱形式
(2)
假設(shè)存在m維時間序列輸入信號x(t)=[x1(t),x2(t),…,xm(t)]T, SFA的目標(biāo)是找到一組特征函數(shù)g(t)=[g1(t),g2(t),…,gm(t)]T,使得特征s(t)=g(x(t))緩慢變化[2]。s(t)=[s1(t),s2(t),…,sm(t)]T用來表示緩慢變化的特征,SFA的優(yōu)化問題表示如下
(3)
約束條件為
〈si〉t=0
(4)
(5)
?i≠j,〈sisj〉t=0
(6)
線性SFA中,每個緩慢變化的特征(slow features, SFs)是所有輸入數(shù)據(jù)的線性組合
s=Wx
(7)
式中,W=[w1,w2,…,wm]T是權(quán)重矩陣。
R=〈x(t)x(t)T〉t=UΛUT
(8)
式中,U是R的特征值,Λ是由特征值組成的對角矩陣。白化矩陣可表示為Q=Λ-1/2UT,白化過程為
z=Λ-1/2UTx=Qx
(9)
結(jié)合式(7)和式(9),SFs進(jìn)一步表示為
s=Wx=WQ-1z=Pz
(10)
式中,P=WQ-1,〈zzT〉t=Q〈xxT〉QT=I,〈z〉t=0。優(yōu)化問題(3)和約束條件進(jìn)一步表示為
(11)
(12)
(13)
(14)
為了滿足式(11)~(14),優(yōu)化問題再次轉(zhuǎn)化為求解正交矩陣P。
(15)
式中,P=[p1,p2,…,pm]是正交特征向量矩陣,對應(yīng)的特征值為Ω=diag(λ1,λ2,…,λm),且λ1≤λ2≤…≤λm。
權(quán)重W表示為
W=PQ=PΛ-1/2UT
(16)
最后,得到m個緩慢程度由大到小排列的SFs。
s=Wx=PQx=PΛ-1/2UTx
(17)
(18)
使用KL散度度量任意兩個變量間的相關(guān)性,并進(jìn)一步構(gòu)造KL散度分量,作為分塊的基礎(chǔ)。
對于訓(xùn)練數(shù)據(jù)X=[x1,x2, …,xm]T∈Rn×m,m和n表示變量數(shù)和樣本個數(shù)。取任意兩個變量xi和xj,且xi∈X,xj∈X,在正態(tài)分布的條件下,KL散度分量表示為
(19)
(20)
最后將訓(xùn)練數(shù)據(jù)分成兩個子模塊X=[x1,x2, …,xm]T= [X1,X2]T。
根據(jù)工業(yè)過程的動態(tài)特性,將子模塊X=[X1,X2]T擴(kuò)展d個時延,使當(dāng)前樣本與過去樣本相關(guān)聯(lián),得到增強(qiáng)的過程矩陣,在SFA建模過程中擴(kuò)展了動態(tài)特性,矩陣增強(qiáng)過程如式(21)所示。
(21)
式中,i是子模塊數(shù)量,n是過程變量的樣本量。
隨后進(jìn)行SFA離線建模,局部建模過程中使用滑動窗口算法訓(xùn)練最優(yōu)的離線模型,如算法1所示。
算法1 滑動窗口訓(xùn)練最優(yōu)模型
以上模型通過將樣本數(shù)據(jù)根據(jù)歐氏距離升序排列,使用滑動窗口算法對樣本進(jìn)行篩選,得到最優(yōu)的局部模型及相應(yīng)的權(quán)重矩陣W1,W2,提高了模型的自適應(yīng)能力。
為了實現(xiàn)對故障的檢測,分別構(gòu)造T2和S2統(tǒng)計量對故障進(jìn)行在線監(jiān)測,T2表示統(tǒng)計量在慢特征空間中的靜態(tài)變化,S2表示過程統(tǒng)計量的動態(tài)變化分布。T2統(tǒng)計量定義為
(22)
(23)
其中,子塊的SFs為si=Wixi,i=1,2。S2統(tǒng)計量定義為
(24)
(25)
算法2 SVDD建模及監(jiān)測
監(jiān)測結(jié)果融合后,得到閾值控制限D(zhuǎn)t=1和半徑中心DR。當(dāng)D≥Dt時,檢測出故障,反之D
基于KL散度的多模塊滑動窗口慢特征分析方法(multi-block moving window slow feature analysis method based on KL divergence,KL-MWSFA)的過程監(jiān)測詳細(xì)步驟如下。
離線訓(xùn)練:
步驟1:對訓(xùn)練數(shù)據(jù)X標(biāo)準(zhǔn)化,計算變量間的KL散度分量dKL。
步驟2:隨機(jī)取dKL中的兩個值作為中心,根據(jù)損失函數(shù)依次迭代更新中心數(shù)值,直到收斂,劃分為兩個子塊X=[X1,X2]T。
步驟3:將子塊進(jìn)行d個時延擴(kuò)展為動態(tài)矩陣,計算與測試數(shù)據(jù)的歐氏距離,將其升序排列。
步驟4:利用滑動窗口對每個子塊進(jìn)行SFA離線建模,分別計算T2,S2的閾值。
在線監(jiān)測:
步驟5:對故障數(shù)據(jù)標(biāo)準(zhǔn)化,擴(kuò)展d個時延,然后分塊。
步驟6:使用SFA進(jìn)行訓(xùn)練,得到T2,S2的閾值,作為SVDD的輸入Y。
步驟7:建立SVDD模型,計算測試向量Yt到超球半徑的距離,得到DR。
步驟8:控制限D(zhuǎn)t=1,當(dāng)D≥Dt時,發(fā)生故障,反之正常。
TE過程是美國化工公司在大量實際工程經(jīng)驗的基礎(chǔ)上,由Downs和Vogel[19]提出的一個化工過程模型,其產(chǎn)生的數(shù)據(jù)具有時變性、強(qiáng)耦合性以及非線性,廣泛應(yīng)用于過程控制、監(jiān)測和故障診斷研究,工藝流程如圖1所示。主要由反應(yīng)器、冷凝器、壓縮機(jī)、氣液分離器以及汽提塔5個操作控制部分構(gòu)成。其共有33個變量,正常情況下的基準(zhǔn)數(shù)據(jù)集包含500個樣本。故障數(shù)據(jù)共有960個樣本,前160個是正常數(shù)據(jù),第161到960個樣本為故障數(shù)據(jù)。本實驗選擇22個過程變量和11個操作變量作為輸入,用500個正常樣本建立模型。仿真共有21個故障類型,其中16種是已知故障,剩余5種是未知故障,實際的建模和監(jiān)測過程中,一般不包括由攪拌速度引起的故障21,故采用15種已知故障進(jìn)行測試,具體描述如表1所示。
圖1 TE過程流程圖Fig.1 TE process flow chart
表1 TE過程故障類型表
隨機(jī)選取兩個變量作為初始中心,計算KL散度分量構(gòu)成的對稱矩陣到初始中心的歐氏距離,離初始中心最近的分量重新確定為新的中心,依次迭代,最后收斂時將變量自動分成兩個子塊。具體分塊策略見第2.1節(jié),變量間的分布如圖2所示,分塊結(jié)果如表2所示。
將KL-MWSFA和核主成分分析(kernel
(a) 聚類前的變量分布(a) Distribution of variables before clustering
principal component analysis, KPCA),ICA,SVDD進(jìn)行對比分析,故障檢測率如表3所示,可以看出,對于故障1,2,4,5,6,7,8,9,10,11,12,13,KL-MWSFA方法都能有效檢測到這些故障。由于故障3,9,15的震級非常小,幾乎所有的方法都無法檢測到故障。但對于故障5,10,12,該方法具有較好的檢測性能。
故障5是冷凝器冷卻水入口溫度的階躍變化,導(dǎo)致冷凝器溫度發(fā)生變化[20]。故障5的監(jiān)測圖如圖3所示,圖3(a)在第300個采樣點以后,無法檢測出故障。圖3(b)I2統(tǒng)計量在第161個采樣點之前發(fā)生虛警現(xiàn)象,圖3(c)在故障正常檢測一段時間后,大部分統(tǒng)計量向控制限下方波動,檢測效果一般。圖3(d)在第161個采樣點檢測到故障后,一直持續(xù)到最后一個樣本,檢測效果較好。
表3 KPCA,DICA,SVDD和KL-MWSFA的故障檢測率Tab.3 Fault detection rate of KPCA, DICA, SVDD and KL-MWSFA
(a) KPCA監(jiān)測結(jié)果(a) Monitoring results of KPCA
(b) DICA監(jiān)測結(jié)果(b) Monitoring results of DICA
(c) SVDD監(jiān)測結(jié)果(c) Monitoring results of SVDD
(d) KL-MWSFA監(jiān)測結(jié)果(d) Monitoring results of KL-MWSFA圖3 故障5監(jiān)測結(jié)果Fig.3 Monitoring results for fault 5
故障10是C進(jìn)料溫度的隨機(jī)變化,主要影響汽提塔壓力[21]。圖4為KPCA,DICA,SVDD和KL-MWSFA的監(jiān)測圖,圖4(a)的T2統(tǒng)計樣本大部分在控制限下面,檢測效果不佳,圖4(c)的監(jiān)控數(shù)據(jù)在控制限附近波動,也無法清楚地檢測到故障的發(fā)生,由表3可見DICA和KL-MWSFA檢測效果較好,但在圖4(b)的子圖中,我們可以看到個別正常樣本出現(xiàn)虛警,只有圖4(d)在檢測故障的同時沒有發(fā)生虛警狀況。
故障12是冷凝器冷卻水入口溫度的隨機(jī)變化[22]。由圖5(a)可知,T2監(jiān)測性能不好且在閾值附近上下波動,圖5(b)雖然可以檢測到故障,但在正常工況時發(fā)生虛警,圖5(c)在發(fā)生虛警的同時監(jiān)測效果沒有圖5(b)的效果好,而圖5(d)的監(jiān)測性能在檢測中均得到改善,沒有虛警并檢測到所有故障。
(a) KPCA監(jiān)測結(jié)果(a) Monitoring results of KPCA
(b) DICA監(jiān)測結(jié)果(b) Monitoring results of DICA
(c) SVDD監(jiān)測結(jié)果(c) Monitoring results of SVDD
(d) KL-MWSFA監(jiān)測結(jié)果(d) Monitoring results of KL-MWSFA圖4 故障10監(jiān)測結(jié)果Fig.4 Monitoring results for fault 10
(a) KPCA監(jiān)測結(jié)果(a) Monitoring results of KPCA
(b) DICA監(jiān)測結(jié)果(b) Monitoring results of DICA
(c) SVDD監(jiān)測結(jié)果(c) Monitoring results of SVDD
(d) KL-MWSFA監(jiān)測結(jié)果(d) Monitoring results of KL-MWSFA圖5 故障12監(jiān)測結(jié)果Fig.5 Monitoring results for fault 12
本文提出了一種基于KL散度的多模塊滑動窗口慢特征分析方法,該方法使用KL散度算法度量每兩個變量間的相似度,利用最小誤差平方和準(zhǔn)則對變量間的距離依次迭代,實現(xiàn)了對樣本的無監(jiān)督分塊,并在每個子塊中建立SFA監(jiān)測模型,解決了單一模塊的不穩(wěn)定性問題,同時引入滑動窗口建立最優(yōu)子模型,利用SVDD將子塊監(jiān)測結(jié)果融合。并將所提方法應(yīng)用于TE過程,驗證了該方法的有效性。