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

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

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

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

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

(1)

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

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

(2)

其中k表示時(shí)延。

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

(3)

(4)

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

(5)

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

(6)

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

(7)

s.t.wTΣXw=1

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

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

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

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

WTW=In∈Rn×n

(8)

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

(9)

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

(10)

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

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

(11)

(12)

(13)

式(12)表示數(shù)據(jù)x在可預(yù)測(cè)主元子空間的投影，式(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ì)量定義分別為：

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

L2=xTWdΛ-1WdTx

(14)

(15)

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

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

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

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

(16)

(17)

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

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

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

圖1 TE過(guò)程流程

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

[1] 李晗,蕭德云.基于數(shù)據(jù)驅(qū)動(dòng)的故障診斷方法綜述[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] 陳玉東,施頌椒.動(dòng)態(tài)系統(tǒng)的故障診斷方法綜述[J].化工自動(dòng)化及儀表, 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.