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

馮曉偉，許劍鋒，何川

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由式(9)可知

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

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

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

2 動態廣義主成分分析方法

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

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

其中，s為延遲量，q為建模數據長度。

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

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

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

對式(23)左乘wT，可得

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

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

3 基于動態廣義主成分分析的故障子空間建模

3.1 故障子空間提取

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

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

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

3.2 故障重構

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

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

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

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

依據文獻[18]可知，正常工業過程的聯合統計量應滿足

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

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

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

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

3.3 故障重構參數的確定

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

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

圖1 s 和d0 計算流程

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

4 基于動態廣義主成分分析的故障診斷

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

4.1 單故障工業過程的故障診斷

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

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

進一步可求得由故障數據在各故障子空間重構后的統計量組成的統計量向量為

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

因此，定義如下統計量變異向量

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

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

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

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

4.2 多故障工業過程的故障診斷

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

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

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

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

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

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

5 數值仿真實驗

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

數值仿真模型為

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

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

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

表1 建模故障數據的故障子空間設置

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

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

圖2 故障2 在重構前后的統計量變化

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

圖3 多故障數據（故障2 和故障5）數據在重構前后的統計量變化

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

圖4 GPCA 算法中故障2 的T2 統計量

圖5 GPCA 算法中故障2 的SPE 統計量

圖6 動態因素提取后，GPCA 算法中故障2 的SPE 統計量

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

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

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

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

6 結束語

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