(浙江理工大學(xué) 機(jī)械與自動(dòng)控制學(xué)院，杭州 310018)

0 引言

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

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

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

1 基于主元分析的故障檢測(cè)

1.1 主元分析

給定一個(gè)具有n個(gè)觀測(cè)值和m個(gè)過(guò)程變量的數(shù)據(jù)集X∈Rn×m，對(duì)X進(jìn)行特征值分解，可得：

(1)

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

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

T=XP

(2)

(3)

(4)

主元數(shù)目的求取可以通過(guò)設(shè)定貢獻(xiàn)度百分比(Cumulative Percent Variance, CPV)來(lái)完成，

(5)

式中,CL為設(shè)定值，滿(mǎn)足上式成立的k的最小值即為滿(mǎn)足當(dāng)前貢獻(xiàn)率閾值的最優(yōu)解。

1.2 主元分析的故障檢測(cè)

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

(6)

(7)

式中,xi為數(shù)據(jù)集在i時(shí)刻的m維觀測(cè)向量，ei為該觀測(cè)向量的殘差向量。

通常，T2統(tǒng)計(jì)量的閾值用自由度為k和n-k的F分布計(jì)算，Q統(tǒng)計(jì)量的閾值則可用χ2分布計(jì)算出來(lái)：

(8)

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

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

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

通過(guò)在每個(gè)觀測(cè)向量后疊加l個(gè)滯后的觀測(cè)向量來(lái)獲得增廣矩陣，構(gòu)建方法如下：

(10)

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

DPCA通過(guò)增廣矩陣可以產(chǎn)生更多的信息關(guān)聯(lián)，比靜態(tài)PCA更適用于動(dòng)態(tài)工業(yè)過(guò)程的故障檢測(cè)。但通過(guò)DPCA得到的主元數(shù)目過(guò)多，使得變量之間的關(guān)系更不容易解釋?zhuān)?dǎo)致計(jì)算量過(guò)大。

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

本文在DPCA的基礎(chǔ)上，提出稀疏動(dòng)態(tài)主元分析法(SDPCA)，通過(guò)結(jié)合了DPCA和SPCA方法的優(yōu)點(diǎn)，以提高故障檢測(cè)的精度，并減少實(shí)時(shí)計(jì)算量，其主要步驟包括：

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

2)確定主元中的非零負(fù)荷數(shù)目；

3)對(duì)增廣矩陣進(jìn)行稀疏主元求解。

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

SubjecttoATA=Ik=k.

(11)

(12)

(13)

由于稀疏后的負(fù)荷向量不一定是正交的，T2統(tǒng)計(jì)量及其控制限的計(jì)算如下：

(14)

(15)

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

2.3 非零負(fù)荷數(shù)目的確定

如何選擇每個(gè)主元的非零負(fù)荷數(shù)目一直是稀疏主元分析的難點(diǎn)。過(guò)多的非零負(fù)荷會(huì)造成計(jì)算上的繁瑣，所以需要盡可能地減少數(shù)目，但同時(shí)又要保證剩下的非零負(fù)荷擁有足夠多的相關(guān)信息。

Jolliffe指出當(dāng)主元的解釋方差相近時(shí)，約束標(biāo)準(zhǔn)的選擇對(duì)模型的最終效果影響不大[17]。在此基礎(chǔ)上，Gajjar等[13]針對(duì)SPCA模型和對(duì)應(yīng)的PCA模型的解釋方差之間的聯(lián)系，提出了一種非零負(fù)荷數(shù)目的前向選擇算法，具體步驟如下：

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

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

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

4)當(dāng)k值與傳統(tǒng)PCA(DPCA)的主元個(gè)數(shù)相同時(shí)，或是稀疏后的主元貢獻(xiàn)度超過(guò)一定閾值時(shí)，停止計(jì)算。否則k=k+1，q=1，返回步驟2)。

其中，η1為當(dāng)前SPCA主元的解釋方差與所對(duì)應(yīng)特征值的比值，計(jì)算公式為：

(16)

但需要注意的是，在 Gajjar等的方法中，每個(gè)稀疏主元的非零負(fù)荷數(shù)目的選擇只和它所對(duì)應(yīng)的一個(gè)特征值有密切聯(lián)系，而忽略了與它之前的特征值的關(guān)系，可能導(dǎo)致解釋方差過(guò)低。為此，本文將此方法進(jìn)行了改進(jìn)，提出了一種新的基于前向選擇的非零負(fù)荷數(shù)目算法如下：

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

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

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

4)當(dāng)稀疏后的主元貢獻(xiàn)度超過(guò)85%時(shí)，停止計(jì)算。否則k=k+1，q=1，返回步驟2)。

其中，η2為當(dāng)前k個(gè)SPCA(或SDPCA)的方差的和與所對(duì)應(yīng)的k特征值的和的比值，計(jì)算公式為：

(17)

在本文提出的算法中，主要依靠來(lái)確定非零負(fù)荷的數(shù)目，因此不僅包含了之前的主元的解釋方差信息，也將參數(shù)考慮進(jìn)來(lái)，以防止當(dāng)前的解釋方差過(guò)低，加強(qiáng)了每個(gè)稀疏主元的非零向量數(shù)目與特征值之間的聯(lián)系。此外，還可直接地通過(guò)調(diào)整貢獻(xiàn)度的大小來(lái)確定主元的數(shù)目。

3 仿真

3.1 數(shù)值仿真

為了表明方法的有效性，采用文獻(xiàn)[4]的數(shù)值例子如下：

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

(18)

(19)

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

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

表1 各種故障檢測(cè)方法的檢測(cè)率 %

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

3.2 TE過(guò)程仿真

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

表2 TE過(guò)程在T2統(tǒng)計(jì)量下的檢測(cè)率 %

這里，將SDPCA方法的動(dòng)態(tài)遲滯設(shè)置為1。經(jīng)過(guò)本文提出的非零負(fù)荷數(shù)目的選擇方法，數(shù)據(jù)的SDPCA模型一共產(chǎn)生了30個(gè)稀疏主元。每個(gè)主元的非零負(fù)荷的個(gè)數(shù)分別為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，貢獻(xiàn)率為85.48%。DPCA方法同樣選取85%的貢獻(xiàn)度，其模型含有24個(gè)主元。SPCA與PCA方法的檢測(cè)數(shù)據(jù)主元數(shù)目均為14。在置信度設(shè)為99%情況下，T2統(tǒng)計(jì)量檢測(cè)的結(jié)果如表2所示，Q統(tǒng)計(jì)量的檢測(cè)結(jié)果如表3所示。表2和表3中，每個(gè)故障的最優(yōu)檢測(cè)率用加粗字體表示。

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

表3 TE過(guò)程在Q統(tǒng)計(jì)量下的檢測(cè)率 %

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

表4 TE過(guò)程的檢測(cè)時(shí)間

圖1 TE過(guò)程故障10在PCA下的檢測(cè)結(jié)果

圖2 TE過(guò)程故障10在DPCA下的檢測(cè)結(jié)果

圖3 TE過(guò)程故障10在SPCA下的檢測(cè)結(jié)果

圖4 TE過(guò)程故障10在SDPCA下的檢測(cè)結(jié)果

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

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