林圣才 楊煜普 屈衛(wèi)東

(上海交通大學(xué)電子信息與電氣工程學(xué)院自動化系系統(tǒng)控制與信息處理教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240)

近二十年來，隨著現(xiàn)代化工及冶金等工業(yè)過程的日益大規(guī)?；蛷?fù)雜化，工業(yè)過程的安全問題越來越受到人們的關(guān)注。復(fù)雜的工業(yè)過程往往難以用精確的物理模型去描述[1]，因此基于多元統(tǒng)計(jì)分析的故障診斷方法應(yīng)運(yùn)而生，并在工業(yè)過程中獲得了成功的應(yīng)用[2～5]。Wise B M等將主元分析方法(PCA)引入了過程監(jiān)控[6]，Lee J M等在PCA的基礎(chǔ)上提出了核PCA方法并將其用于故障診斷[7]。PCA方法是從觀測數(shù)據(jù)中提取與統(tǒng)計(jì)無關(guān)的主元，通過構(gòu)造統(tǒng)計(jì)量對過程狀況進(jìn)行監(jiān)控統(tǒng)計(jì)，判斷過程是否出現(xiàn)故障，它要求數(shù)據(jù)服從高斯分布。但是實(shí)際工業(yè)過程往往并不滿足這個條件，同時PCA方法無法反映過程的動態(tài)時序特性，這在一定程度上影響了它的故障檢測準(zhǔn)確率?？深A(yù)測元分析(Forecastable Component Analysis,ForeCA)作為一種新的統(tǒng)計(jì)信號處理方法[8]，克服了這個不足。它是一種全新的用于多變量時序相關(guān)信號的降維與特征提取方法，它能從已有的數(shù)據(jù)中捕捉到系統(tǒng)的動態(tài)特性，并以此來預(yù)測系統(tǒng)運(yùn)行變化的趨勢，因此所提取的特征更能從本質(zhì)的上描述工業(yè)過程。

筆者將可預(yù)測元分析方法引入到故障檢測中，通過所挖掘的可預(yù)測元提取出觀測信號中的可預(yù)測分量，構(gòu)造兩種統(tǒng)計(jì)量對其進(jìn)行統(tǒng)計(jì)監(jiān)控。該方法克服了主元分析方法需要數(shù)據(jù)服從高斯分布且無法反映過程時序特性的不足，能夠預(yù)測系統(tǒng)運(yùn)行變化的趨勢，反映出系統(tǒng)的動態(tài)特性，提升故障檢測的效果。在TE過程上的仿真結(jié)果表明了該方法的可行性和有效性。

設(shè)矩陣X∈Rn×m，可預(yù)測元分析的基本思想是尋找到一個線性變換WT∈Rk×n，使得：

(1)

W為負(fù)荷矩陣，它的列向量表示負(fù)荷向量，彼此相互正交。

γy(k)=E(yt-μy)(yt-k-μy)T,k∈R

(2)

其中k表示時延。

定義單變量平穩(wěn)過程的譜密度為對其自協(xié)方差函數(shù)的傅里葉變換：

(3)

(4)

熵越大則平穩(wěn)過程的后續(xù)變化越難被預(yù)測，且白噪聲無法被預(yù)測，因此可得：

(5)

根據(jù)式(5)定義平穩(wěn)過程的可預(yù)測度為：

(6)

對于多變量二階平穩(wěn)過程Xt，考慮線性變換yt=wTXt，其中w(w∈Rn)是W的列向量，即可預(yù)測元，此時yt就可以看成是一個單變量的二階平穩(wěn)過程。Goerg G給出了ForeCA的最優(yōu)化問題[8]：

(7)

s.t.wTΣXw=1

在求解式(7)問題時，首先使用加權(quán)交疊平均(WOSA)譜估計(jì)法對隨機(jī)過程進(jìn)行譜密度估計(jì)[9]，然后使用EM-Like算法求取可預(yù)測元[10]。通過文獻(xiàn)[8]給出的算法可以計(jì)算出一組按照可預(yù)測度由高到低順序排列的可預(yù)測元(可預(yù)測元個數(shù)可以指定，一般不大于平穩(wěn)過程的變量個數(shù))，進(jìn)而得到線性變換矩陣WT。

2 基于可預(yù)測元分析的故障檢測

首先選取一段正常工況生產(chǎn)下的觀測數(shù)據(jù)Yn×m，其中n為變量個數(shù)，m為采樣點(diǎn)數(shù)(時間序列)，由于變量使用的量綱不同，因此需要對觀測數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，處理后的數(shù)據(jù)記為Xn×m。對Xn×m運(yùn)用ForeCA算法，選取可預(yù)測元的個數(shù)等于觀測數(shù)據(jù)中變量的個數(shù)n，使用算法后得負(fù)荷向量wi∈Rn,i=1,2,…,n，每個負(fù)荷向量wi對應(yīng)的可預(yù)測度為Ωi，進(jìn)而得到線性變換矩陣：

WT=[w1,w2,…,wn]T∈Rn×n

WTW=In∈Rn×n

(8)

定義累積可預(yù)測度貢獻(xiàn)率為：

(9)

根據(jù)式(9)的定義，一般取Ψ(Ω)≥85%，由此可求得k值。定義前k個可預(yù)測元對應(yīng)的負(fù)荷向量wi(wi∈Rn,i=1,2,…,k)為可預(yù)測主元負(fù)荷向量，并令：

(10)

定義WdT為可預(yù)測主元負(fù)荷矩陣。通過從可預(yù)測元矩陣中選取可預(yù)測主元一方面可降低矩陣的維數(shù)和計(jì)算量，另一方面可構(gòu)造新的用于故障檢測的統(tǒng)計(jì)量。式(9)給出了一種選取可預(yù)測主元的方法，除此之外，還可以使用交叉驗(yàn)證的方式。

當(dāng)用于在線數(shù)據(jù)時，設(shè)某次采樣得到的數(shù)據(jù)為x,x∈Rn，可得：

(11)

(12)

(13)

式(12)表示數(shù)據(jù)x在可預(yù)測主元子空間的投影，式(13)表示數(shù)據(jù)x在殘差子空間的投影。

根據(jù)式(11)、(13)定義兩種統(tǒng)計(jì)指標(biāo)：L2統(tǒng)計(jì)量與SPE統(tǒng)計(jì)量。L2統(tǒng)計(jì)量與SPE統(tǒng)計(jì)量定義分別為：

通常，我們都會不自覺地將一張攝影作品歸入某個特定的時代框架里。然而荒誕的是，這張照片似乎想要從一切短暫的年代歸屬之中退出去。一方面，觀眾感覺這張照片反映的是一個詭異而險(xiǎn)惡的“舊世界”。它以某種方式將這種痕跡留在照片上：斑點(diǎn)、劃痕，諸如此類來自玻璃負(fù)片時代的特征。威特金常常采用一種高度直覺化的方式來完成照片制作的物理過程，比如刮擦負(fù)片、漂白或是調(diào)節(jié)畫面顏色。

L2=xTWdΛ-1WdTx

(14)

(15)

其中Λ表示由前k個可預(yù)測主元對應(yīng)的可預(yù)測度組成的對角陣。L2統(tǒng)計(jì)量是通過可預(yù)測模型內(nèi)部的可預(yù)測元模的波動來反映系統(tǒng)的變化情況，SPE統(tǒng)計(jì)量則表示一個觀測數(shù)據(jù)到可預(yù)測模型空間的距離，反映了測量值對模型的偏離程度。

由于所構(gòu)造的統(tǒng)計(jì)量并不一定嚴(yán)格服從正態(tài)分布，因此可以采用核密度估計(jì)法[11,12]對統(tǒng)計(jì)量進(jìn)行密度估計(jì)，選取合適的置信水平進(jìn)而確定統(tǒng)計(jì)量的控制限。

當(dāng)使用ForeCA算法提取出可預(yù)測系統(tǒng)運(yùn)行變化趨勢的特征，并構(gòu)造出上述兩種統(tǒng)計(jì)量后，將此可預(yù)測模型運(yùn)用于在線數(shù)據(jù)，對其進(jìn)行檢驗(yàn)：如果檢驗(yàn)結(jié)果在相應(yīng)統(tǒng)計(jì)量的控制限以下，則說明目前系統(tǒng)工作在可預(yù)測模型所預(yù)測的變化范圍之內(nèi)，即系統(tǒng)工作正常；反之，則說明目前系統(tǒng)的工作狀態(tài)已經(jīng)偏離可預(yù)測模型所預(yù)測的變化范圍，因此有理由判斷系統(tǒng)已經(jīng)出現(xiàn)了故障。

當(dāng)使用L2統(tǒng)計(jì)量與SPE統(tǒng)計(jì)量檢測到過程系統(tǒng)出現(xiàn)故障后，需要定位出系統(tǒng)發(fā)生異常的位置，這里采用貢獻(xiàn)圖法[13]來處理這個問題。對于L2統(tǒng)計(jì)量與SPE統(tǒng)計(jì)量，對每個變量定義如下的貢獻(xiàn)值：

(16)

(17)

利用貢獻(xiàn)圖法可以知道過程變量對當(dāng)前狀態(tài)的貢獻(xiàn)，貢獻(xiàn)值最大的變量很可能是這次故障發(fā)生的位置所在，因此可以及時給予報(bào)警。

3 TE平臺仿真實(shí)驗(yàn)

TE實(shí)驗(yàn)平臺是Downs J J和Vogel E F于1993年提出的[14]。TE實(shí)驗(yàn)平臺在Eastman化學(xué)公司的世界工藝流程上做了少許改動，可以很好地模擬現(xiàn)實(shí)中的復(fù)雜工況，其流程如圖1所示。

圖1 TE過程流程

選取正常的樣本數(shù)據(jù)500個，每個樣本點(diǎn)包含33個變量，分別為22個連續(xù)變量XMEAS(1)～XMEAS(22)和前11個控制變量XMV(1)～XMV(11)，變量的具體含義參見文獻(xiàn)[14]。首先對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，使用ForeCA算法對數(shù)據(jù)進(jìn)行處理，之后對測試數(shù)據(jù)進(jìn)行檢測，測試數(shù)據(jù)集包含960個樣本點(diǎn)，每個樣本點(diǎn)包含33個變量。前160個樣本點(diǎn)為正常數(shù)據(jù)，后800個樣本點(diǎn)為故障數(shù)據(jù)。

將ForeCA法和傳統(tǒng)PCA法進(jìn)行對比，在PCA法中，選取T2統(tǒng)計(jì)量與SPE統(tǒng)計(jì)量，表1列出了8個故障的檢測準(zhǔn)確率。

表1 ForeCA和PCA方法故障檢測準(zhǔn)確率比較 %

從表1中可以看出，L2統(tǒng)計(jì)量的檢測準(zhǔn)確率要高于T2統(tǒng)計(jì)量，F(xiàn)oreCA算法的SPE統(tǒng)計(jì)量和PCA算法的SPE統(tǒng)計(jì)量相比各有優(yōu)勢。圖2給出了IDV(10)和IDV(20)的故障檢測準(zhǔn)確率對比圖。仿真實(shí)驗(yàn)表明了ForeCA方法在故障檢測中的可行性與有效性。

a. IDV(10) b. IDV(20)

選取IDV(10)作為典型故障進(jìn)行詳細(xì)分析。IDV(10)的發(fā)生是由于TE過程中供料C的溫度產(chǎn)生了隨機(jī)變化。當(dāng)系統(tǒng)某時刻出現(xiàn)此變化后，控制回路會補(bǔ)償這個變化，進(jìn)而導(dǎo)致過程的總體特征產(chǎn)生變化。圖3顯示的是一個可預(yù)測元提取的IDV(10)特征，虛線表示前160個正常數(shù)據(jù)的均值，點(diǎn)劃線為后800個故障數(shù)據(jù)的均值，兩者的均值幾乎相同，但是故障數(shù)據(jù)的方差產(chǎn)生了非常大的變化，這說明提取的特征很好地抓住了過程變化的總體方差特征，同時顯示出一定的周期性，因此可對過程下一次的狀態(tài)做出合理預(yù)測，這兩點(diǎn)都為ForeCA用于故障檢測提供了保證。圖3提取的IDV(10)特征也正好符合L2與SPE統(tǒng)計(jì)量所顯示的IDV(10)檢測圖的兩個波峰的特點(diǎn)。

圖3 可預(yù)測元提取的IDV(10)特征

仍然以IDV(10)為例。假設(shè)在檢測出系統(tǒng)發(fā)生故障后，使用貢獻(xiàn)圖法確定故障變量(圖4)，圖中變量18的貢獻(xiàn)值最大，由此可以推斷故障很可能是變量18異常導(dǎo)致的，即解吸塔溫度異常。綜合考慮TE過程的所有故障，只有IDV(10)能直接導(dǎo)致解吸塔溫度發(fā)生變化，因此在很大程度上可以認(rèn)為系統(tǒng)發(fā)生了IDV(10)。

圖4 IDV(10)發(fā)生時貢獻(xiàn)圖法的診斷結(jié)果

4 結(jié)束語

針對傳統(tǒng)PCA算法具有的需要數(shù)據(jù)服從高斯分布且丟失過程動態(tài)特性的缺點(diǎn)，將ForeCA應(yīng)用于過程監(jiān)控領(lǐng)域，選取可預(yù)測主元，構(gòu)造新的統(tǒng)計(jì)量，建立了完整的基于ForeCA的故障診斷方法。ForeCA方法可以從觀測數(shù)據(jù)中提取過程的動態(tài)特性，預(yù)測過程以后的運(yùn)行變化趨勢，這在一定程度上提高了建模精度，因而具有較好的故障檢測能力。最后在TE過程上的仿真結(jié)果表明了ForeCA方法的可行性和有效性。在后續(xù)工作中，可以通過改進(jìn)監(jiān)控統(tǒng)計(jì)量和ForeCA算法來進(jìn)一步提高故障檢測能力。

[1] 李晗,蕭德云.基于數(shù)據(jù)驅(qū)動的故障診斷方法綜述[J].控制與決策,2011,26(1):1～9.

[2] Kimura D,Nii M,Yamaguchi T,et al.Fuzzy Nonlinear Regression Analysis Using Fuzzified Neural Networks for Fault Diagnosis of Chemical Plants[J]. JACIII, 2011, 15(3): 336～344.

[3] Kano M, Nakagawa Y. Data-based Process Monitoring, Process Control, and Quality Improvement: Recent Developments and Applications in Steel Industry[J]. Computers & Chemical Engineering, 2008, 32(1): 12～24.

[4] 陳玉東,施頌椒.動態(tài)系統(tǒng)的故障診斷方法綜述[J].化工自動化及儀表, 2001, 28(3): 1～14.

[5] Zhang Y,Zhang Y.Fault Detection of Non-Gaussian Processes Based on Modified Independent Component Analysis[J]. Chemical Engineering Science, 2010, 65(16): 4630～4639.

[6] Wise B M, Ricker N L, Veltkamp D F, et al. A Theoretical Basis for the Use of Principal Component Models for Monitoring Multivariate Processes[J]. Process Control and Quality,1990,1(1): 41～51.

[7] Lee J M, Yoo C K, Choi S W, et al. Nonlinear Process Monitoring Using Kernel Principal Component Analysis[J].Chemical Engineering Science,2004, 59(1): 223～234.

[8] Goerg G.Forecastable Component Analysis[C].Proceedings of the 30th International Conference on Machine Learning. Atlanta,GA,USA:ICML,2013: 64～72.

[9] Nuttall A H,Carter G C.Spectral Estimation Using Combined Time and Lag Weighting[J].Proceedings of the IEEE,1982, 70(9): 1115～1125.

[10] Dempster A P, Laird N M, Rubin D B. Maximum Likelihood from Incomplete Data via the EM Algorithm[J].Journal of the Royal Statistical Society, 1977, 39(1): 1～38.

[11] Botev Z I, Grotowski J F, Kroese D P. Kernel Density Estimation via Diffusion[J]. The Annals of Statistics, 2010, 38(5): 2916～2957.

[12] Wand M P, Jones M C. Kernel Smoothing[M]. Boca Raton: Crc Press,1994.

[13] Yoon S, MacGregor J F. Fault Diagnosis with Multivariate Statistical Models Part I: Using Steady State Fault Signatures[J]. Journal of Process Control, 2001, 11(4): 387～400.

[14] Downs J J, Vogel E F. A Plant-wide Industrial Process Control Problem[J]. Computers & Chemical Engineering, 1993, 17(3): 245～255.