動態(tài)廣義主成分分析及其在故障子空間建模中的應用

馮曉偉，許劍鋒，何川

（火箭軍工程大學，陜西 西安 710025）

0 引言

大型工業(yè)設備在運行過程中的突發(fā)性故障容易造成經(jīng)濟損失甚至危及生命財產(chǎn)安全。在現(xiàn)代工業(yè)技術高速發(fā)展的今天，日趨復雜的系統(tǒng)結構使系統(tǒng)的安全性和可靠性面臨著前所未有的挑戰(zhàn)，也促使在線工業(yè)過程監(jiān)控和故障診斷方法越來越受到重視[1-2]。

目前，學者們針對工業(yè)過程的故障診斷技術做了大量研究，提出了多種不同的在線故障診斷方案。基于人工神經(jīng)網(wǎng)絡的方法[3-4]采用歷史數(shù)據(jù)訓練網(wǎng)絡模型來實現(xiàn)故障定位。以主成分分析（PCA,principal component analysis）[5-6]和偏最小二乘（PLS,partial least square）[7-8]為典型代表的多元統(tǒng)計分析方法通過歷史數(shù)據(jù)建立主子空間和殘差子空間，并在兩類子空間中分別構建T2和平方預測誤差（SPE,squared prediction error）統(tǒng)計量來實現(xiàn)過程監(jiān)控。在此基礎上，學者們又發(fā)展了核PCA（KPCA,kernel PCA）[9-10]和核PLS（KPLS,kernel PLS）[11]等方法將數(shù)據(jù)映射到高維空間來實現(xiàn)非線性工業(yè)過程監(jiān)控。傳統(tǒng)方法將工業(yè)過程視為靜態(tài)過程，但實際工業(yè)過程往往具有動態(tài)特性，即系統(tǒng)當前狀態(tài)與歷史狀態(tài)密切相關。因此，傳統(tǒng)靜態(tài)監(jiān)控方法在處理動態(tài)過程數(shù)據(jù)時容易出現(xiàn)較高的漏報率或誤報率。許多動態(tài)多元統(tǒng)計方法，如遞歸PCA（RPCA,recursive PCA）[12]、動態(tài)PCA（DPCA,dynamic PCA）[13-14]以及動態(tài)PLS（DPLS,dynamic PLS）[15-16]等方法的相繼提出，提升了動態(tài)工業(yè)過程的監(jiān)測準確度。

上述方法只考慮對系統(tǒng)狀態(tài)進行監(jiān)測，即確定系統(tǒng)是否發(fā)生了故障。在監(jiān)測到系統(tǒng)發(fā)生故障后，為了進一步確定故障發(fā)生的位置、類型和時間，一些故障診斷方法被相繼提出。其中，故障重構方法通過分離故障數(shù)據(jù)的正常工況和故障工況實現(xiàn)故障診斷，在實際工業(yè)過程監(jiān)控中得到了廣泛應用。基于重構的貢獻圖法[17]通過分析各變量的重構貢獻圖，可以定位到發(fā)生故障的傳感器，在較嚴重的故障拖尾效應下仍能表現(xiàn)出很好的性能。文獻[18]提出通過對標準化的故障數(shù)據(jù)直接做奇異值分解來提取故障子空間。基于閾值的故障子空間提取方法[19]在系統(tǒng)含有加性噪聲時比奇異值分解方法具有更小的重構誤差。文獻[20]通過相對變化分析來尋找故障方向，可實現(xiàn)多重故障重構。

基于奇異值分解的故障子空間提取方法直接針對包含正常工況和故障工況的原始數(shù)據(jù)進行分解，其不足之處是對數(shù)據(jù)特征分解不夠充分，難以得到精確的故障子空間。針對傳統(tǒng)故障子空間提取方法的不足，文獻[21]提出了一種基于廣義主成分分析（GPCA,generalized PCA）的故障子空間提取方法，通過剔除正常建模數(shù)據(jù)和故障建模數(shù)據(jù)之間的公共關系得到了更精確的故障子空間。針對動態(tài)工業(yè)過程的故障診斷問題，文獻[22]提出將數(shù)據(jù)進行時移并重組到更高維度來提取動態(tài)主成分，取得了良好效果。在此研究基礎上，本文提出一種動態(tài)廣義主成分分析（DGPCA,dynamic GPCA）方法并將其運用到故障子空間建模與故障診斷中。首先，分別將正常數(shù)據(jù)和故障數(shù)據(jù)重組到高維空間，并采用DGPCA 方法提取正常數(shù)據(jù)和故障數(shù)據(jù)之間的共性部分，實現(xiàn)對原始數(shù)據(jù)動態(tài)廣義主成分信息的提取；然后，剔除該共性部分得到純故障數(shù)據(jù)，并使用該純故障數(shù)據(jù)實現(xiàn)精確故障子空間提取和故障重構；最后，針對不同故障提取相應故障子空間建立故障子空間庫，并基于該故障子空間庫實現(xiàn)了在線故障診斷。所提方法實現(xiàn)了故障子空間提取流程的簡化和精度的提升，且直接構造故障子空間庫的方法降低了故障診斷的計算復雜度。

1 基于PCA 的故障子空間重構方法

令 X∈Rm×l為正常工況下的采樣數(shù)據(jù)集，其中m為變量個數(shù)，l為樣本長度。假設X已經(jīng)過標準化處理，即X各維數(shù)據(jù)均值為0，方差為1[23]。通過對X進行主成分分析，可得如下模型

其中，P∈Rm×r、T∈Rr×l分別為數(shù)據(jù)在主子空間的負載矩陣和得分矩陣，I∈Rm×m為單位矩陣，E∈R(m-r)×l為數(shù)據(jù)在殘差子空間的得分矩陣，r為主子空間維數(shù)。

當獲取新的采樣數(shù)據(jù)x后，將x在主/次子空間內進行分解可得

其中，Λ=diag(λ1,λ2,…,λr)，其對角元素由主子空間特征值 λi(i=1,2,…,r)組成。顯然，統(tǒng)計量 T2和SPE 分別用于監(jiān)控主子空間和殘差子空間的工況。如果滿足

在實際工業(yè)生產(chǎn)過程中，可能會發(fā)生多種不同的故障，而不同故障之間存在互不相同的特征信息。通過對故障數(shù)據(jù)的深入分析，區(qū)分不同故障的特征信息，可以實現(xiàn)故障診斷，即定位到具體故障及故障發(fā)生的時間。故障重構[18,20]是一種可實現(xiàn)故障定位與診斷的有效手段。該方法的核心思想是，當故障發(fā)生時，將故障數(shù)據(jù)分解為正常工況部分 x*和純故障工況部分 x-。顯然，如果能夠估計出純故障工況部分，則可以得到正常工況數(shù)據(jù)，進而實現(xiàn)故障重構，即

其中，Ξ是由正交向量組成的故障子空間代表故障方向； f 代表故障幅度。

如果故障子空間Ξ已知，那么故障重構的關鍵就是通過故障數(shù)據(jù)估計故障幅度f。在主子空間，要最大限度地消除x中故障信息的影響，故障幅度fT應滿足[25]

在進行故障重構后，使用重構數(shù)據(jù)重新計算2 個統(tǒng)計量。如果數(shù)據(jù)被正確地重構，那么重新計算的統(tǒng)計量應處于控制限以下。

假設工業(yè)過程中共有n類故障，且故障j(j=1,2,…,n)的故障子空間為Ξj。已獲得標準化的故障j的故障數(shù)據(jù)集Xj，x(k)為Xj中第k個樣本，即

由式(9)可知

因為 x*(k)為零均值數(shù)列，如果故障幅度 f(k)足夠大，則可以近似認為[18]

式(15)表明，Ξj具有與Xj相同的列空間[18]。于是，對數(shù)據(jù)列Xj做奇異值分解，得

其中，當確定故障子空間維度lj后即故障j的故障子空間。

2 動態(tài)廣義主成分分析方法

傳統(tǒng)靜態(tài)方法在包含動態(tài)因素的工業(yè)過程監(jiān)控中容易出現(xiàn)漏報或誤報的情況，而基于PCA 的故障重構方法因對故障分離不充分而未能對故障子空間精確建模。為此，本文提出了一種可在動態(tài)因素提取的基礎上實現(xiàn)精確故障子空間建模的DGPCA 方法。

假設 X∈Rm×l，Y(j)∈Rm×l分別代表用于訓練的正常樣本集和故障樣本集，其中 Y(j)代表第j(j=1,2,…,n)類故障的樣本集。令xk∈Rm，yk(j)∈Rm分別代表X和 Y(j)中第k個樣本數(shù)據(jù)。為了提取樣本數(shù)據(jù)中的動態(tài)信息，對樣本數(shù)據(jù)重組如下[22,26]

其中，s為延遲量，q為建模數(shù)據(jù)長度。

根據(jù)前面的分析，故障樣本應由兩部分組成，即正常工況部分和純故障工況部分。其中，正常工況部分的特征子空間應與正常樣本的特征子空間相關性最大。尋找該相關性，即尋找故障樣本集在正常樣本集X0中的最大投影。因此，提出如下目標函數(shù)

為了求解目標函數(shù)的最大值，定義拉格朗日函數(shù)

將式(21)對w求偏導，并令其為0，得

對式(23)左乘wT，可得

上述分析表明，要尋找的特征向量w本質上是對應矩陣束最大廣義特征值的廣義特征向量。

故障數(shù)據(jù)和正常數(shù)據(jù)的共性部分代表了故障數(shù)據(jù)中隱含的正常工況，是影響故障子空間建模精度的直接因素。DGPCA 方法實現(xiàn)了動態(tài)工業(yè)過程中對故障數(shù)據(jù)和正常數(shù)據(jù)共性部分的提取，為故障子空間的精確提取奠定了基礎條件。

3 基于動態(tài)廣義主成分分析的故障子空間建模

3.1 故障子空間提取

顯然，對純故障數(shù)據(jù)進行奇異值分解，即可得到更加準確的故障子空間。對做奇異值分解，有

當確定故障空間維度dd后，U(j)前dd個故障方向即故障j的故障子空間Ξ(j)。

在基于重構的故障診斷方法中，精確的故障重構是有效實現(xiàn)故障診斷的基礎。基于DGPCA 方法實現(xiàn)了對故障子空間的精確提取，為后續(xù)精確故障重構的實現(xiàn)提供了條件支持。

3.2 故障重構

為了判定故障重構效果，對正常工況數(shù)據(jù)集X0的期望矩陣做特征分解，構建直接DPCA 監(jiān)測模型[26]，即

取U0的前d0列構成主子空間P0，其余部分構成次子空間，并分別求 Y(j)在主/次子空間的統(tǒng)計量為

文獻[21,25]都是通過將故障數(shù)據(jù)分別投影到主子空間和殘差子空間來求故障幅度，然后分別進行故障重構。本文采用文獻[18]的聯(lián)合指標構建模式來計算統(tǒng)計量和控制限。聯(lián)合統(tǒng)計量計算方式為

其中，Φ 為對稱正定矩陣，有

依據(jù)文獻[18]可知，正常工業(yè)過程的聯(lián)合統(tǒng)計量應滿足

其中，λi為矩陣R0Φ 的第i個特征值，其重數(shù)為νi。

當獲取新的待監(jiān)測采樣數(shù)據(jù)y后，根據(jù)延遲量，聯(lián)合歷史數(shù)據(jù)采用式(18)將其擴展為高維數(shù)據(jù)在聯(lián)合指標限定下，其故障幅度計算式為

顯然，如果故障數(shù)據(jù)只包含單個故障j，則基于故障子空間Ξ(j)采用式(39)重構后的數(shù)據(jù)，其統(tǒng)計量應當處于控制限式(35)以下。反之，則不能將故障數(shù)據(jù)重構在控制限以下。

不同類型的故障的其故障子空間互不相同，因此可將其作為相應故障的特征信息。顯然，當且僅當故障類型與故障子空間相對應時才能成功將該故障數(shù)據(jù)重構到控制限以下。因此，基于故障子空間的故障重構方法為故障診斷提供了決策依據(jù)。

3.3 故障重構參數(shù)的確定

只有合理選擇上述延遲量s和幾種子空間維數(shù)d0、ds和dd等參數(shù)才能對故障子空間進行準確的提取。本節(jié)討論幾種維度參數(shù)的選擇方法。

延遲量s代表工業(yè)過程中的動態(tài)階數(shù)，d0代表工業(yè)過程中變量的非線性關系數(shù)量。兩者可通過文獻[22]提出的互相關分析法來確定，如圖1 所示。

圖1 s 和d0 計算流程

參數(shù)ds為的維度，代表故障建模數(shù)據(jù)和正常建模數(shù)據(jù)之間相關維度，令

4 基于動態(tài)廣義主成分分析的故障診斷

在實際工業(yè)過程中，可能發(fā)生的故障種類是多種多樣的，在檢測到系統(tǒng)故障后，需要進一步對故障進行定位，即故障診斷。本文故障診斷的基本思想為，首先通過DGPCA 方法對所有已知故障進行故障子空間提取，建立故障子空間庫；當發(fā)生故障后，嘗試采用故障子空間庫對該故障進行重構，通過對比重構前后統(tǒng)計量的變化，即可判斷出故障類型。根據(jù)同時發(fā)生的故障數(shù)量，可將故障診斷分為單故障診斷和多故障診斷。

4.1 單故障工業(yè)過程的故障診斷

在已知的n類故障中，對于每一類故障，都可以基于歷史故障數(shù)據(jù)提取其故障子空間，得到Ξ(1),Ξ(2),…,Ξ(n)。定義矩陣

代表n類故障組成的故障子空間庫。當獲取新的待監(jiān)測采樣數(shù)據(jù)y后，基于式(38)，采用式(18)將其擴展為高維數(shù)據(jù)?并在故障子空間庫中求其在每一類故障中的故障幅度為

進一步可求得由故障數(shù)據(jù)在各故障子空間重構后的統(tǒng)計量組成的統(tǒng)計量向量為

其中，diag(·)代表將矩陣的對角元素組合成向量。顯然，如果待檢測過程數(shù)據(jù)只含有第j類故障，則在重構后只有的統(tǒng)計量能恢復到控制限以下，而在其他故障子空間的故障重構則不能顯著改變統(tǒng)計量，即

因此，定義如下統(tǒng)計量變異向量

其中，n1 為所有元素均為1 的n維行向量，sφ 代表故障樣本數(shù)據(jù)在重構前基于式(33)求得的統(tǒng)計量，顯然，只有在正確的故障子空間進行故障重構，由于重構前后的統(tǒng)計量的顯著變化才能導致顯著增加，而在其他故障子空間重構的數(shù)據(jù)由于其統(tǒng)計量變化不顯著，導致的值近似為0。因此，在實際故障診斷過程中，可以將 ψc的各個元素作為故障指標單獨繪圖，實現(xiàn)故障診斷。ψc中對應位置的指標指示對應的故障信息。

如果系統(tǒng)同時發(fā)生多種故障，那么故障數(shù)據(jù)在任意單一故障子空間中都不能實現(xiàn)完全重構，因此重構數(shù)據(jù)的統(tǒng)計量不能降低到控制限以下。然而，由于在非故障子空間重構時的統(tǒng)計量變化較小，而在所有包含故障的子空間中重構的統(tǒng)計量變化仍然較顯著，此時 ψc中將會有多個指標出現(xiàn)顯著變化。因此，采用指標 ψc可實現(xiàn)多故障診斷。

討論1 工業(yè)過程故障診斷總體可以分為離線建模和在線監(jiān)測兩部分，而在線監(jiān)測算法的計算復雜度因其涉及硬件執(zhí)行效率而成為關注點之一。假設式(41)中矩陣Ξ的總列數(shù)為p，則式(42)的故障幅度的計算復雜度為

2(n+1)p；在式(42)計算結果的基礎上，式(43)的故障重構計算復雜度為 msnp2，而式(44)的統(tǒng)計量向量計算復雜度為 (ms+1)msn3。綜上，在忽略較小項后，算法的總計算復雜度約為

4.2 多故障工業(yè)過程的故障診斷

上述過程是針對單故障診斷分析的，在實際工業(yè)過程中，可能同時發(fā)生幾類故障。文獻[25]指出，在發(fā)生多故障情況時，將過程數(shù)據(jù)投影到任何一個單一故障子空間都不能實現(xiàn)完全重構，而應同時投影到由相應故障類型的故障子空間組合成的聯(lián)合故障子空間來進行故障重構。不難發(fā)現(xiàn)，雖然此時過程數(shù)據(jù)不能在單一故障子空間完全重構，但是只要在過程數(shù)據(jù)中所包含的故障類型的故障子空間重構時，重構數(shù)據(jù)的統(tǒng)計量都會有較顯著的減小，其減小幅度相對于在非對應故障類型的故障子空間重構時的微弱減小幅度而言十分明顯。因此，可以通過式(46)所求指標向量ψc中各元素的幅度來判斷過程數(shù)據(jù)所包含的故障類型，進而在這些故障類型的聯(lián)合故障子空間中進行故障重構。

假設所提取的k類故障的故障子空間分別為Ξ1～Ξk，定義如下矩陣

為當前故障的聯(lián)合故障子空間，注意此處Ξc不一定是列滿秩的。為了消除各故障子空間之間的線性關系，對Ξc做奇異值分解，得到

則取Uc的前r列Uc(:,1:r)為聯(lián)合故障子空間，仍記為Ξc，此處r為 Λc對角線中非零元素的個數(shù)。當獲取聯(lián)合故障子空間后，即可在聯(lián)合故障子空間中對數(shù)據(jù)進行故障重構

另一方面，指標 φc,1,φc,2,…,φc,k也同時用于多故障工業(yè)工程中的故障診斷。

討論2 已有故障庫需基于已知類型故障建立，且建立的故障庫只能對已知類型的故障進行區(qū)分。在實際工業(yè)過程中，由于實際工況的復雜性，可能會發(fā)生未知類型的新故障。根據(jù)前面的分析可知，在發(fā)生新故障后，由于未對該故障建模，因此系統(tǒng)只能監(jiān)測到該故障而無法區(qū)分其故障類型。針對此類情況，可以在監(jiān)測到未知類型故障后，先給出故障報警信號并對故障數(shù)據(jù)進行采集，在獲得足夠的故障數(shù)據(jù)集后，對其進行建模并更新故障庫，即可在后續(xù)監(jiān)控過程中實現(xiàn)對該故障的識別。

5 數(shù)值仿真實驗

本節(jié)通過數(shù)值仿真實驗來驗證本文所提方法在故障子空間建模和故障診斷的性能。實驗中對比了文獻[21]所提靜態(tài)廣義主成分分析算法在動態(tài)工業(yè)過程中的故障重構能力。

數(shù)值仿真模型為

其中，tk∈R2×100、uk∈R5×100為中間變量，tk代表2 個獨立的隨機變量，且服從標準正態(tài)分布；xk∈R5×100為傳感器輸出的數(shù)據(jù)點；G∈R5×2、C1,C2∈R5×5為建立過程狀態(tài)與傳感器測量的測量矩陣；e ,v∈R5×100均代表5 個獨立的高斯噪聲，均值為0，方差為0.01。測量矩陣G、C1、C2均隨機產(chǎn)生，在本文實驗中，其值分別為

對于建模的故障數(shù)據(jù)，共設置5 種故障，分別發(fā)生在5 個變量上，故障幅度均為5，即

5 種故障的故障子空間設置如表1 所示。

表1 建模故障數(shù)據(jù)的故障子空間設置

對于測試的單故障數(shù)據(jù)，設置為在30 采樣點以后發(fā)生故障2，故障幅度為7；對于測試的多故障數(shù)據(jù)，設置為在30 個采樣點以后同時發(fā)生故障2和故障5，故障幅度分別為6、9。在本文實驗中，延遲量s=1，d0=2，5 個故障的ds均為4，故障子空間維度dd分別為2、1、2、1、1。

故障2 在重構前后的統(tǒng)計量變化如圖2 所示。圖2 中包括重構前的統(tǒng)計量，以及分別采用故障2和故障3 的故障子空間對故障2 進行重構的結果。仿真結果表明，重構前的數(shù)據(jù)在30 個采樣點以后統(tǒng)計量明顯超出控制限，并在控制限以上振蕩，表明系統(tǒng)在30 個采樣點以后出現(xiàn)故障。采用故障2的故障子空間對故障數(shù)據(jù)進行重構后，重構數(shù)據(jù)的統(tǒng)計量都在控制限以下，說明重構后的數(shù)據(jù)已不再包含故障信息。然而，采用故障3 的故障子空間對故障2 的故障數(shù)據(jù)進行重構，由于故障子空間和發(fā)生的故障類型不一致，因此重構對于數(shù)據(jù)的統(tǒng)計量未產(chǎn)生明顯影響。

圖2 故障2 在重構前后的統(tǒng)計量變化

圖3 所示為針對包含故障2 和故障5 的多故障數(shù)據(jù)在重構前后的統(tǒng)計量變化，圖3 中包括重構前的統(tǒng)計量，分別采用故障2 和故障5 的故障子空間以及故障2 和故障5 的聯(lián)合故障子空間對故障數(shù)據(jù)重構的結果。仿真結果表明，對于包含故障2 和故障5 的多故障數(shù)據(jù)，采用故障2 和故障5 的故障子空間均能有效降低故障數(shù)據(jù)的統(tǒng)計量，但都不能將數(shù)據(jù)重構到控制限以下，采用故障2 和故障5 的聯(lián)合故障子空間則可以有效將故障重構數(shù)據(jù)重構到控制限以下。

圖3 多故障數(shù)據(jù)（故障2 和故障5）數(shù)據(jù)在重構前后的統(tǒng)計量變化

采用文獻[21]提出的廣義主成分分析方法對故障2 進行過程監(jiān)控的T2統(tǒng)計量和SPE 統(tǒng)計量結果分別如圖4 和圖5 所示。通過結果可以發(fā)現(xiàn)，由于沒有考慮動態(tài)因素的提取，傳統(tǒng)靜態(tài)廣義主成分分析方法得到的統(tǒng)計量均處于控制限以下，即不能實現(xiàn)動態(tài)過程下的過程監(jiān)控。顯然，基于式(17)和式(18)將數(shù)據(jù)重組到高維空間后，再采用該方法即可實現(xiàn)動態(tài)因素的提取。在采取動態(tài)因素提取的方式后，再采用GPCA方法重新計算相應統(tǒng)計量，如圖6 所示。此時可以發(fā)現(xiàn)，重新計算的SPE 統(tǒng)計量可實現(xiàn)對故障2 的監(jiān)測。

圖4 GPCA 算法中故障2 的T2 統(tǒng)計量

圖5 GPCA 算法中故障2 的SPE 統(tǒng)計量

圖6 動態(tài)因素提取后，GPCA 算法中故障2 的SPE 統(tǒng)計量

采用本文方法對故障2 的故障診斷結果如圖7所示。從圖7 中可以看到，在30 個采樣點后，統(tǒng)計量超過控制限，意味著此時系統(tǒng)發(fā)生了故障，而故障2 的統(tǒng)計量變化指標的明顯增大意味著此時發(fā)生的故障為故障2。

圖7 單故障過程（故障2）的故障診斷結果

采用本文方法對包含故障2 和故障5 的多故障過程故障診斷結果如圖8 所示。從圖8 中可以看到，在30 個采樣點后，統(tǒng)計量超過控制限，指示當前過程已發(fā)生故障。進一步，故障2 和故障5 的統(tǒng)計量變化指標明顯增大，表明此時系統(tǒng)同時發(fā)生的是故障2 和故障5，與故障預設情況一致。

圖8 多故障過程（故障2 和故障5）的故障診斷結果

6 結束語

通過動態(tài)廣義主成分分析方法，可以提取動態(tài)工業(yè)過程中正常數(shù)據(jù)和故障數(shù)據(jù)之間的公共部分，即故障數(shù)據(jù)中的正常工況信息，從而實現(xiàn)故障數(shù)據(jù)中正常工況和故障工況的分離。對獲得的只包含純故障信息的故障數(shù)據(jù)，采用奇異值分解方法可準確提取故障子空間，提高了故障子空間的建模精度。在實際工業(yè)過程監(jiān)控中，采用提取的故障子空間可對故障數(shù)據(jù)進行有效故障重構。進一步，可采用已知故障數(shù)據(jù)庫建立故障子空間庫，從而實現(xiàn)故障診斷。