(浙江理工大學 機械與自動控制學院，杭州 310018)

0 引言

現代工業(yè)過程正朝著智能化的方向發(fā)展，為了保證產品的質量、避免生產安全事故的發(fā)生，實時高效的故障檢測系統在工業(yè)過程中的作用變得越來越重要。主元分析(Principal Component Analysis, PCA) 作為一種簡單實用的多元統計過程監(jiān)控方法，能夠有效地剔除冗余信息，只需要正常工況下的歷史數據即可建立模型，在工業(yè)過程的故障檢測和診斷中得了廣泛的應用[1-3]。考慮到過程變量本身的時序自相關性，Ku等[4]提出了動態(tài)主元分析方法(Dynamic Principal Component Analysis, DPCA)，通過用帶有時間滯后特性的變量構造出動態(tài)數據矩陣后，利用PCA方法提取相關主元，提高了故障檢測的精度。

但是通過PCA方法得到的負荷向量中的元素通常是非零的，不僅不利于特征的解釋和提取，計算量也非常大。為了獲取稀疏負荷向量，稀疏主元分析方法(Sparse Principal Component Analysis, SPCA) 近年來得到了極大的發(fā)展，并在生物學、醫(yī)學等各個領域得到了廣泛的應用[5-6]。Jolliffe[7]提出了SCoTLASS算法來獲得稀疏主元。Zou等[8]提出利用LASSO懲罰項來獲得稀疏主元。Leng等[9]提出了簡單自適應主元分析，用自適應的LASSO懲罰項取代傳統的Elastic Net。Kang等[10]進一步的提出了自適應稀疏主元分析(ASPCA)。彭必燦等[11]，劉洋等[12]，Gajjar等[13]采用SPCA實現了工業(yè)過程的故障檢測及診斷。Gajjar等[13]并提出一種前向選擇方法來確定稀疏主元的非零負荷數目。

本文結合DPCA和SPCA兩者的特點，提出了一種基于稀疏動態(tài)主元分析法(Sparse Dynamic Principal Component Analysis, SDPCA)的故障檢測方法。該方法先通過疊加時間滯后變量的方式建立測量數據的動態(tài)增廣矩陣，然后再通過Lasso約束函數獲取增廣矩陣的稀疏主元。此外，本文還提出了一種新的稀疏主元非零負荷數目的確定方法，該方法考慮到了稀疏主元之間的關聯性，對文獻[13]中提出的前向選擇法進行了改進，在保證累積貢獻率的情況下，進一步提高了稀疏度。最后，通過數值仿真和田納西伊斯曼過程來驗證所提方法的性能，并與PCA、DPCA和SPCA方法進行對比。

1 基于主元分析的故障檢測

1.1 主元分析

給定一個具有n個觀測值和m個過程變量的數據集X∈Rn×m，對X進行特征值分解，可得：

(1)

式中，Λ∈Rm×m為特征值的對角矩陣，且對角數值沿對角線遞減。V∈Rm×m為酉矩陣，其列向量為主元的負荷向量。

選取V的前k列得到一個新的負荷矩陣P∈Rm×k，可得到得分矩陣：

T=XP

(2)

(3)

(4)

主元數目的求取可以通過設定貢獻度百分比(Cumulative Percent Variance, CPV)來完成，

(5)

式中,CL為設定值，滿足上式成立的k的最小值即為滿足當前貢獻率閾值的最優(yōu)解。

1.2 主元分析的故障檢測

HotellingT2和Q統計量常被用于基于多變量統計法的故障檢測。T2統計量由計算主元的得分向量在空間中的馬氏距離獲得，而Q統計量則是殘差向量在殘差空間中的歐式距離，

(6)

(7)

式中,xi為數據集在i時刻的m維觀測向量，ei為該觀測向量的殘差向量。

通常，T2統計量的閾值用自由度為k和n-k的F分布計算，Q統計量的閾值則可用χ2分布計算出來：

(8)

2 基于稀疏動態(tài)主元分析的故障檢測

2.1 動態(tài)主元分析

工業(yè)過程一般都具有較強的動態(tài)特性，因此過程數據具有強烈的序列相關性。為了解決工業(yè)過程的動態(tài)特性，DPCA通過擴展觀測矩陣，可以消除數據的自相關性，從而提高故障檢測的精度[14-15]。

通過在每個觀測向量后疊加l個滯后的觀測向量來獲得增廣矩陣，構建方法如下：

(10)

式中是數據集在t時刻的m維觀測值，n是總的樣本數目。然后將該增廣矩陣代替原觀測矩陣來進行主元分析。一般情況下，在工業(yè)過程控制中序列的滯后參數選擇1或2[16]。

DPCA通過增廣矩陣可以產生更多的信息關聯，比靜態(tài)PCA更適用于動態(tài)工業(yè)過程的故障檢測。但通過DPCA得到的主元數目過多，使得變量之間的關系更不容易解釋，并導致計算量過大。

2.2 基于稀疏動態(tài)主元分析的故障檢測

本文在DPCA的基礎上，提出稀疏動態(tài)主元分析法(SDPCA)，通過結合了DPCA和SPCA方法的優(yōu)點，以提高故障檢測的精度，并減少實時計算量，其主要步驟包括：

1)獲得動態(tài)增廣矩陣；

2)確定主元中的非零負荷數目；

3)對增廣矩陣進行稀疏主元求解。

首先，通過公式(10)獲得設計矩陣，這里參數可采用交叉驗證進行選取。然后再確定非零負荷的數目，本文提出了一種新的非零負荷數目的確定方法，并在下一小節(jié)中給出具體的過程。最后，在確定的非零負荷數目下，將降維問題轉化為回歸最優(yōu)化問題，對其進行稀疏求解。這里，我們采用類似于Zou等[8]提出的SPCA中的優(yōu)化算法來獲取增廣矩陣的稀疏主元，即在DPCA模型上增加LASSO懲罰項：

SubjecttoATA=Ik=k.

(11)

(12)

(13)

由于稀疏后的負荷向量不一定是正交的，T2統計量及其控制限的計算如下：

(14)

(15)

其中，γ是稀疏得分向量的協方差矩陣。而Q統計量的計算方式則與PCA模型一致，可采用公式(7)和公式(9)。

2.3 非零負荷數目的確定

如何選擇每個主元的非零負荷數目一直是稀疏主元分析的難點。過多的非零負荷會造成計算上的繁瑣，所以需要盡可能地減少數目，但同時又要保證剩下的非零負荷擁有足夠多的相關信息。

Jolliffe指出當主元的解釋方差相近時，約束標準的選擇對模型的最終效果影響不大[17]。在此基礎上，Gajjar等[13]針對SPCA模型和對應的PCA模型的解釋方差之間的聯系，提出了一種非零負荷數目的前向選擇算法，具體步驟如下：

1)通過公式(1)求解得特征值矩陣。

2)計算可獲得稀疏主元的公式(11)，其中主元個數為k，最后一個主元的非零負荷數目為q(k和q的初始值為1) 。

3)如果η1≥90%，則確定q為最后一個負荷向量的非零負荷數目；否則q=q+1，返回步驟2)。

4)當k值與傳統PCA(DPCA)的主元個數相同時，或是稀疏后的主元貢獻度超過一定閾值時，停止計算。否則k=k+1，q=1，返回步驟2)。

其中，η1為當前SPCA主元的解釋方差與所對應特征值的比值，計算公式為：

(16)

但需要注意的是，在 Gajjar等的方法中，每個稀疏主元的非零負荷數目的選擇只和它所對應的一個特征值有密切聯系，而忽略了與它之前的特征值的關系，可能導致解釋方差過低。為此，本文將此方法進行了改進，提出了一種新的基于前向選擇的非零負荷數目算法如下：

1)通過公式(1)求解得特征值矩陣。

2)計算可獲得稀疏主元的公式(11)，其中主元個數為k，最后一個主元的非零負荷數目為q(k和q的初始值為1) 。

3)如果η1≥80%且η2≥90%，則確定q為最后一個負荷向量的非零負荷數目；否則q=q+1，返回步驟2)。

4)當稀疏后的主元貢獻度超過85%時，停止計算。否則k=k+1，q=1，返回步驟2)。

其中，η2為當前k個SPCA(或SDPCA)的方差的和與所對應的k特征值的和的比值，計算公式為：

(17)

在本文提出的算法中，主要依靠來確定非零負荷的數目，因此不僅包含了之前的主元的解釋方差信息，也將參數考慮進來，以防止當前的解釋方差過低，加強了每個稀疏主元的非零向量數目與特征值之間的聯系。此外，還可直接地通過調整貢獻度的大小來確定主元的數目。

3 仿真

3.1 數值仿真

為了表明方法的有效性，采用文獻[4]的數值例子如下：

y(k)=z(k)+v(k).

(18)

(19)

其中：輸入變量w是均值為零，方差為1的白噪聲，v是均值為零，方差為0.1的白噪聲。輸入變量u和輸出變量y，均為可測變量。采集正常情況下的100個樣本，用于建模。為了模擬故障，在非正常情況下，在第10個樣本上引入故障。其中，故障1是，w1變量添加一個單位階躍。故障2是，將階躍的幅值增加至3。

為了減少隨機變量的影響，本次實驗一共生成了1000組數據集。在設置PCA和SPCA的主元個數為3，DPCA和SDPCA的主元個數為4時，DPCA和SDPCA的動態(tài)遲滯設置為1。此時PCA和DPCA模型的貢獻率可達到98%以上，SPCA的貢獻率在92%至96%之間，SDPCA的貢獻率在89%至94%之間。4種方法的置信度均設置為99%。表1為分別用PCA、SPCA、DPCA和SDPCA4種方法所得到的故障檢測率(Fault Detection Rate, FDR)平均值。

表1 各種故障檢測方法的檢測率 %

從表1中可以看出，4種故障檢測方法在T2統計量下的檢測結果差別并不大，但SDPCA方法用Q統計量測得的檢測率明顯高于其它3種方法。并且從SPCA與PCA、SDPCA與DPCA的對比可以看出，經過稀疏處理之后的檢測效果會有小的提升。

3.2 TE過程仿真

田納西-伊斯曼過程是一種模仿真實化學工業(yè)現場的標準控制過程，被廣泛的用作控制、優(yōu)化及故障診斷的仿真對象。該過程主要包括了5個操作單元，分別是反應器、冷凝器、氣液分離器、循環(huán)壓縮機、汽提塔。過程數據的工況被人為設定21種故障情況，每種故障情況以3分鐘為采樣間隔分別采集了960組采樣數據，每組數據含有53個特征變量，每組數據的在第161個采樣時刻引入不同的故障。本文選取了第1至22個測量變量及第42至第52個控制變量，共33個特征變量用于故障檢測。

表2 TE過程在T2統計量下的檢測率 %

這里，將SDPCA方法的動態(tài)遲滯設置為1。經過本文提出的非零負荷數目的選擇方法，數據的SDPCA模型一共產生了30個稀疏主元。每個主元的非零負荷的個數分別為24, 23, 9, 4, 5, 4, 2, 2, 2, 2, 4, 2, 2, 3, 2, 2, 2, 1, 1, 2, 1, 1, 1, 2, 1, 4, 2, 2, 2 和 2，貢獻率為85.48%。DPCA方法同樣選取85%的貢獻度，其模型含有24個主元。SPCA與PCA方法的檢測數據主元數目均為14。在置信度設為99%情況下，T2統計量檢測的結果如表2所示，Q統計量的檢測結果如表3所示。表2和表3中，每個故障的最優(yōu)檢測率用加粗字體表示。

通過表2可以看出，4種方法在T2統計檢測結果中，SDPCA在其中的16個故障數據下的檢測效果都是最優(yōu)的。特別是故障3、故障9和故障15這幾個傳統方法難以檢測的故障，在其它3種方法下的檢測率均低于2%，而通過SDPCA方法在T2統計量下可以得到高于5%的檢測率。通過表3可以看出，在Q統計檢測結果中，SDPCA方法在11個故障下的檢測率都是最高的。尤其在故障10和故障16的檢測中，SDPCA在兩種統計量下的檢測效果均明顯優(yōu)于其它3種方法。圖1～4為故障10在4種方法下的Q統計量檢測結果。

表3 TE過程在Q統計量下的檢測率 %

此外，稀疏算法可以有效減少非零負荷的數目， TE過程的數據經過稀疏動態(tài)處理之后，SDPCA模型的所有主元中所包含的非零負荷數目為116，而DPCA模型下共包含了792個非零負荷。非零負荷的大量減少，使得SDPCA相比DPCA具有更高的實時計算效率。表4列出了同一臺服務器使用MATLAB軟件分別用4種方法采用T2統計量處理TE過程故障數據所用的時間。從表中可以看出經過數據動態(tài)處理后的SDPCA和DPCA兩種方法所花費的時間更長一些，但經過稀疏處理后的SDPCA方法比DPCA方法的檢測時間更短一些。

表4 TE過程的檢測時間

圖1 TE過程故障10在PCA下的檢測結果

圖2 TE過程故障10在DPCA下的檢測結果

圖3 TE過程故障10在SPCA下的檢測結果

圖4 TE過程故障10在SDPCA下的檢測結果

4 結束語

結合動態(tài)主元分析和稀疏主元分析的特點，本文提出了一種基于稀疏動態(tài)主元分析的故障檢測方法，并進一步改進了非零負荷數目的前向選擇算法。該方法不僅考慮到了實際工業(yè)數據的動態(tài)特性，同時也大量地減少了非零負荷的數量，從而提高了計算效率，并增加了特征的可解釋性。通過數值和TE過程的仿真結果表明，所提方法能夠獲得良好的動態(tài)過程的故障檢測性能。未來的研究將進一步研究非零負荷數目的優(yōu)化選擇，以及將動態(tài)稀疏方法推廣到其它的多變量統計方法中。