基于稀疏性非負矩陣分解的故障監測方法

王帆，楊雅偉，譚帥，侍洪波

（華東理工大學化工過程先進控制和優化技術教育部重點實驗室，上海 200237）

引 言

隨著現代工業過程變得越來越復雜，許多工業過程，如化工過程，通常包含大量高維并且互相關聯的數據。在這種情況下，多元統計過程監控（MSPM）方法已經被提出，用來從過程數據中提取有用信息，并將其用來進行系統的故障監測。傳統的基于數據的多元統計過程監控方法主要包括主成分分析方法（PCA）、偏最小二乘方法（PLS）、獨立成分分析方法（ICA）。這些方法已經被成功地用于化工過程的故障監測。

主成分分析是過程監測與故障診斷研究應用最廣泛的一種方法，它通過將得到的樣本數據方差進行最優線性降維處理，解決了變量之間的相關性[1]。主成分分析方法提取的特征可以反映出原始數據的全局信息[2]。基于PCA 來進行故障監測的成果已有很多[3-5]。

然而，基于PCA 的故障監測方法，假定過程的潛變量服從多元高斯分布，這與實際的工業情況是不符的，實際過程的一部分變量并不服從高斯分布，這樣用PCA 對過程進行監控，勢必會導致監測性能不好[6]。而且經過PCA 處理之后，潛變量是互相正交且不相關的。最近出現了一種新的降維方法，叫做非負矩陣分解（NMF）[7]，它沒有對潛變量的性質進行假設，除了要求其非負。NMF 通過從原始數據中找到一個低秩矩陣逼近，可以產生對原始數據的局部特征的稀疏表示。NMF 在最近幾年吸引了越來越多的關注，在許多應用中都有良好的表現，比如數據挖掘和模式識別[8-9]。非負矩陣分解的方法要求非負性，這樣在進行分解時就只允許加而非減的組合，使非負矩陣分解方法可以學習局部的基于部分的表示[10]。每一個潛變量激活一個可見變量的子集，或者說是激活一個“部分”。所有的潛變量激活作用疊加在一起，產生一個整體[2]。從這個角度看，NMF 可以找到一個中間表示（基于部分的），介于數據的局部和全局結構之間[11]。

在化工過程故障監測領域，非負矩陣分解方法的應用很少。Li 等[11]提出了應用傳統的非負矩陣分解方法進行非高斯過程的故障監測，用仿真實驗證明了非負矩陣分解方法可以用于化工過程的故障監測，而且比基于主成分分析的方法有更廣泛的應用能力。隨后，Li 等[2]又改進了PNMF（投影非負矩陣分解），提出了GNMP（廣義非負矩陣投影），主要目的是減輕對原始數據的非負約束。但是，Li 等雖然在文中提到了NMF 可以增強對數據的壓縮性和可解釋性是由于其基于部分的稀疏表示，也提到了稀疏性是數據壓縮的關鍵屬性，但并沒有提出關于施加稀疏性的改進。因此，本文提出了采用SNMF（稀疏性非負矩陣分解）[12]方法來進行化工過程故障監測。由于在非負矩陣分解基礎上引入了稀疏編碼（sparse coding）方法，SNMF 可以得到對數據集更稀疏的表示。另外，對傳統的非負矩陣分解得到的原始數據低秩近似矩陣，即系數矩陣H，進行正交化處理，從而在降維時除去變量中的冗余信息，得到更少的投影方向。因為施加了稀疏性和正交化，信息集中到了更少的投影方向上，所以相比于傳統的非負矩陣分解方法具有優越性。同時，本文提出的稀疏性非負矩陣分解，其算法中的目標函數和更新策略不同于Li 等所使用的方法。最后用仿真實驗證明了在應用于故障監測時，基于稀疏性非負矩陣分解（SNMF）的方法具有可行性。

1 稀疏性非負矩陣分解

1.1 非負矩陣分解

非負矩陣分解是一種線性的、非負的近似數據表示。給定一個非負數據矩陣X?Rm×n，NMF 可以找到一個近似的分解，將X分解為兩個非負矩陣的乘積，如式(1)所示，

其中，W?Rm×k，H?Rk×n，k通常選擇為(m+n)k≤mn[7]。W是由基向量wi組成的基矩陣，H是由系數向量hi組成的系數矩陣。經過這種分解方法得到的基矩陣W和系數矩陣H都具有一定的稀疏性[7]。

這樣，非負矩陣分解問題就歸結為一個非線性優化問題，用一個目標函數來測度低秩近似的效果。Lee 等[10]定義了兩種簡單的目標函數，如式(2)、式(3)所示

Li 等[2,11]使用的是式(2)所示的目標函數，本文使用的是式(3)所示的目標函數。目標函數的不同，決定了相應的更新策略也是不同的。Lee 等[10]提出了一種迭代方法來解決這個問題，如式(4)、式(5)所示

這個方法叫做乘法更新（MU）算法，可以保證非負性，但是有可能會卡在非平穩點[13]。使用乘法更新算法解決優化問題，就可以得到W和H。基矩陣W保留了原始數據的空間關系和數據結構，系數矩陣H可以被看作是原始數據矩陣的低秩近似。從統計過程監測方法的角度來看，需要觀察矩陣H來反映工業過程的狀態[11]。因此，基于非負矩陣分解（NMF）的監測模型如式(6)所示

1.2 稀疏編碼

稀疏編碼的概念指的是一種代表方法。在這種代表方法中，少數單位就可以有效地代表典型數據向量，這就意味著，在表示時，大部分單位的價值接近于零，只有少部分單位作用明顯，價值是顯著的非零值[14]。經過非負矩陣分解后得到的矩陣本身具有稀疏性特點，因此，將稀疏編碼的方法與非負矩陣分解結合是很合理的。對非負矩陣分解施加一定的稀疏性約束，可以得到更稀疏的表示。

線性稀疏編碼模型的思想如下：假設觀察到的數據集是一些基向量的線性組合，基于這些基向量的系數代表了為作進一步分析的新的特征向量。稀疏編碼方法假定這些潛變量呈現出稀疏性，然后通過尋找一種特定形式的低熵編碼來最大化稀疏性，這種低熵編碼的概率密度的高峰出現在零周圍并且有嚴重的厚尾性[15]。

用 I(x,y) 來表示一個圖像，它可以表示為一些基函數的線性組合[12]，如式(9)所示

其中，fi(x,y)是基函數，a i是系數。如上所述，線性稀疏編碼可以表示為一個優化問題，通過最小化式(10)這個損失函數[12]

其中，λ是正常數，決定了第2 項相對于第1項的重要性。第1 項測度了編碼描述圖像的好壞，由式(11)來表示[15]

1.3 稀疏性非負矩陣分解

傳統的非負矩陣分解方法，只對矩陣元素施加了非負約束，然而，對矩陣元素施加稀疏性約束，就可以得到更稀疏的表示，即降維后得到更少的投影方向。因此，把稀疏編碼和非負矩陣分解結合在一起就可以得到稀疏性非負矩陣分解。進而，考慮到傳統的非負矩陣分解得到的原始數據低秩近似，即系數矩陣H，不具有正交性，這就意味著變量之間存在冗余信息，從而考慮對其進行正交化處理。這樣在降維時，可以將變量中的冗余信息除去，得到更少的投影方向。使信息集中到更少的投影方向上，與稀疏性的思想也是契合的。此時需要解決的優化問題如式(12)所示[12]

相應的乘法更新（MU）方法如式(13)～式(15)所示

這種方法是通過最小化所有Hkj的和來進行稀疏性的限制，從而對傳統的非負矩陣分解方法施加了稀疏性，將其改進成為稀疏性非負矩陣分解。對矩陣H的正交化處理體現在約束條件中的HHΤ=I，在更新算法求解過程中，體現在式(15)。沒有對矩陣H施加正交處理時，式(15)原本應該是，在迭代更新時以HXΤ代替WΤ，從而得到式(15)所示的更新迭代公式。這樣代替是因為X=WH，兩邊轉置得到XΤ=HΤWΤ，兩邊左乘矩陣H得到HXΤ=HHΤWΤ，約束條件HHΤ=I，所以在迭代更新時以HXΤ代替WΤ，這樣可以保證其正交性。

然后就可以得到如式(6)所示的監測模型。相應地，監測指標如式(7)和式(8)所示。

1.4 稀疏度比較

本文使用的稀疏性定量表示是基于L1 范數和L2 范數之間的關系[14]，如式(16)所示

其中，W?Rm×k，Wi是基矩陣W的列向量。稀疏度的值會有兩種極端情況：如果Wi中所有元素都相等且是非零值，則稀疏度為零；如果Wi中只有一個元素是非零值，其余元素值都為零，則稀疏度為1。

對于基矩陣W，其稀疏度由式(17)計算

使用隨機生成的矩陣進行計算驗證。用MATLAB 隨機產生100×10 的矩陣進行計算。NMF和SNMF 兩種方法條件相同，都取k為5，循環次數為1000 次。對矩陣W和H的初始化都采用隨機初始化方法。在SNMF 方法中，α的值選取為10。每個方法運行10 次，通過式(16)和式(17)來計算基矩陣W的稀疏度，最后選用10 次數據的平均結果。計算結果如表1所示。

表1 NMF 和SNMF 方法中W 的稀疏度Table 1 Sparseness of W in NMF and SNMF

如表1所示，在隨機計算中，經過SNMF 方法得到的基矩陣W的稀疏度大于傳統NMF 方法得到的稀疏度。通過稀疏度比較可以證明，由于施加了稀疏性的約束和進行了正交化處理，SNMF 可以得到比傳統NMF 更加稀疏的表示，從而具有了對數據更好的表示能力。

2 基于稀疏性非負矩陣分解的故障監測方法

基于上述提出的稀疏性非負矩陣分解，這里提出一種新的故障監測方法。

2.1 初始化

作為一種迭代算法，對W和H的初始化就顯得很重要。因為不同的初始化條件可能會導致不同的結果。正隨機數初始化和奇異值分解（SVD）是比較常見的初始化方法，其中SVD 是初始化基矩陣W的較好方法[16]。

本文采用PCA 中的負荷矩陣P的絕對值作為基矩陣W的初始化，這樣做方便了兩種算法的比較，而且PCA 選取主元的數目可以作為選擇基矩陣W維數大小的參考。系數矩陣H的初始化采用正隨機數的方法。

2.2 監測指標

如前文所示，基于稀疏性非負矩陣分解（SNMF）的監測模型如式(6)所示。相應地，監測指標如式(7)、式(8)所示。

2.3 控制限

核密度估計（KDE）被用來確定兩個監測指標N2和SPE 的控制上限。核密度估計是一類數據驅動的技術，用于密度函數的非參數估計。這是一個很有力的工具，可以從給定的樣本集中得到密度函數的經驗分布[17]。N2和SPE 兩個統計量都是單變量，所以很適合用KDE 這個方法來進行計算。

樣本集X={xi,i=1,2,…,n}，密度函數為P(x)，可以用式(18)表示

對單變量核估計的KDE 方程如式(19)所示[18]

2.4 故障監測流程

基于稀疏性非負矩陣分解（SNMF）方法的故障監測流程如下。

建模：

（1）初始化，正常樣本數據X經過PCA 后得到負荷矩陣P，用來初始化W。H用正隨機數方法來進行初始化；

（2）用稀疏性非負矩陣分解（SNMF）方法來計算得到W和H；

（3）根據式(7)和式(8)分別計算統計量N2和SPE；

（4）使用KDE 計算出兩個監測指標N2和SPE的控制上限。

監測：

（1）根據測試數據得到新的系數矩陣H；

（2）根據式(7)和式(8)分別計算新的N2和SPE；

（3）通過比較新的N2和SPE 與相應的控制上限，判斷其是否為故障。

3 仿真實驗

3.1 TE 過程

TE 過程可以人為設定21 種故障，已經被廣泛地用于評估和比較過程監測方法的有效性。TE 過程主要包括反應器、冷凝器、循環壓縮機、解吸塔和氣液分離器5 個單元。有4 種反應物和2 種產品，此外，還包含1 種副產品和1 種惰性氣體。該過程共有41 個測量變量和12 個操作變量。41 個測量變量中，22 個是連續的過程測量變量，19 個是非連續的過程測量變量。Downs 等[19]給出了關于TE 過程的詳細信息。TE 過程的工藝流程圖如圖1所示。

TE 過程可以仿真正常工作狀態和21 種故障情況。21 種故障類型如表2所示。

仿真數據的采樣間隔為3 min，每一個樣本數據集包括52 個變量。正常數據作為訓練數據用來建模。對每一個故障情況，進行960 次采樣，故障從第161 次采樣開始產生，這樣就得到了測試數據集。

這里分別將基于PCA 的監測方法、基于NMF的監測方法和基于SNMF 的監測方法，應用于TE過程，進行比較。在基于PCA 的方法中，選擇方差貢獻度為85%，此時，對應的k值為27。所以在基于NMF 和基于SNMF 的方法中，相應地將k值選取為27。另外，這兩種方法的循環次數都是1000次。在SNMF 方法中，α的值選取為10。這是因為，由迭代更新公式(13)～式(15)可以看出，當α的值選取為0 時，SNMF 退化為傳統的NMF；而當α的值選取過大時，不利于迭代算法的進行。α的取值過小和過大都對結果有不利的影響，所以經過嘗試，α的值選取為10。

圖1 TE 過程的工藝流程Fig.1 Process flow sheet of Tennessee Eastman process

表2 TE 過程21 種故障類型Table 2 Process faults for TE process

表3 PCA，NMF，SNMF 3 種方法在TE 過程的故障監測率Table 3 Fault detection rates of PCA,NMF and SNMF in TE process

3.2 結果討論

故障監測率作為評價監測好壞的標準。3 種方法的監測結果如表3所示。

由于故障3、9 和15 的測試數據在均值、方差和高階矩上沒有可觀察的變化，所以統計方法的故障監測率都很低[20]，這里沒有進行分析討論。對于故障18，由于NMF 和SNMF 這兩種方法的誤報率過高，所以不予討論。對于故障21，由于PCA 的監測指標SPE 的誤報率過高，所以也不予討論。對于其余的故障情況，基于NMF 和SNMF 的方法，誤報率始終在5%以內，而基于PCA 的方法的SPE指標，在故障8、12、17 和19 時，誤報率都稍微超過了5%。可以從表3中看出，對于故障1、2、6、7、8，3 種方法得到的結果基本相同，這是因為這些故障情況的故障程度很明顯，所以3 種方法都能取得很好的故障監測率。

對于故障4、5、10、16、19 和20，基于SNMF的方法都好于基于PCA 的方法。有些情況下，SNMF的結果不如PCA，這是因為SNMF 降維后的投影方向相當于PCA 的主元，因為稀疏性約束和正交化處理，基于SNMF 的方法將信息集中到了較少的投影方向上，可能會丟失一些與這些投影方向相關性差的信息。而PCA 的主元數較多，相比之下，包含更多的信息。但同時也注意到，在這些情況下，SNMF的結果與PCA 相差不多，說明SNMF 用很少的投影方向依然捕捉到了相關信息。

對于故障4、5、10、11、13、14、16、17、19和20，基于SNMF 的方法都好于基于NMF 的方法。故障4、11、16 和20 的N2指標，SNMF 相比于NMF有很明顯的提高（超過10%），在故障11 上提高了26.125%。最值得注意的是在故障4 上的表現，基于PCA 的方法，監測指標T2的故障監測率很低，而基于NMF 和SNMF 的方法，監測指標N2有較高的故障監測率。在故障4 的監測上，也可以很明顯地體現出基于SNMF的方法好于基本的NMF方法。3 種方法對故障4 的監測情況如圖2所示。

從整體上看，基于NMF 和SNMF 的方法好于基于PCA 的方法，而且由于加入了稀疏性的約束和正交化的處理，基于SNMF 的方法可以得到更稀疏的表示，因此其監測效果較基于傳統的NMF 的方法，有了一定程度的提高。

4 結 論

圖2 PCA，NMF 和SNMF 對故障4 的監測結果Fig.2 Monitoring results of PCA,NMF and SNMF for fault 4

相比于主成分分析（PCA），非負矩陣分解（NMF）是一種更普遍的方法，因為這種方法對潛變量沒有除了非負要求之外的任何假設。在非負性限制的基礎上，再加入稀疏性約束和正交化處理，可以更好地提高這種方法對數據集的表示能力。本文提出了使用稀疏性非負矩陣分解（SNMF）的方法來進行化工過程的故障監測，因為將稀疏編碼的概念引入到傳統的非負矩陣分解方法上，所以可以得到對數據集更稀疏的表示能力，另外，又對低秩近似矩陣進行正交化處理使得降維時除去了變量的冗余信息，將信息集中在更少的投影方向上，進而在故障監測時得到更好的效果。仿真實驗的結果表明，在TE 過程的監控上，基于稀疏性非負矩陣分解（SNMF）的故障監測方法好于傳統的非負矩陣分解（NMF）方法和主成分分析（PCA）方法。然而，非負矩陣分解的方法在近幾年還在不斷的發展中，以后會出現更多關于此方法的改進。另外，關于該方法應用于化工過程的故障監測時，還有一些值得考慮的問題，比如初始化問題和收斂性問題。

[1]Xia Luyue (夏陸岳),Pan Haitian (潘海天),Zhou Mengfei (周猛飛),Cai Yijun (蔡亦軍),Sun Xiaofang (孫小方).Process monitoring and fault diagnosis of propylene polymerization based on improved multiscale principle component analysis [J].CIESC Journal(化工學報),2011,62 (8):2312-2317

[2]Li Xiangbao,Yang Yupu,Zhang Weidong.Statistical process monitoringviageneralized non-negative matrix projection [J].Chemometrics and Intelligent Laboratory Systems,2013,121:15-25

[3]Ge Zhiqiang,Song Zhihuan.Distributed PCA model for plant-wide process monitoring [J].Industrial & Engineering Chemistry Research,2013,52:1947-1957

[4]Ge Zhiqiang,Song Zhihuan.Mixture Bayesian regularization method of PPCA for multimode process monitoring [J].AIChE Journal,2010,56:2838-2849

[5]Jia Mingxing,Xu Hengyuan,Liu Xiaofei,Wang Ning.The optimization of the kind and parameters of kernel function in KPCA for process monitoring [J].Computers and Chemical Engineering,2012,46:94-104

[6]Ge Zhiqiang (葛志強),Song Zhihuan (宋執環).New online monitoring method for multiple operating modes process [J].Journal of Chemical Industry and Engineering（China） (化工學報),2008,59 (1):135-141

[7]Lee D D,Seung H S.Learning the parts of objects by nonnegative matrix factorization [J].Nature,1999,401:788-791

[8]Shang Fanhua,Jiao L C,Wang Fei.Graph dual regularization non-negative matrix factorization for co-clustering [J].Pattern Recognition,2012,45:2237-2250

[9]Zeng Kun,Yu Jun,Li Cuihua,You Jane,Jin Taisong.Image clustering by hyper-graph regularized non-negative matrix factorization [J].Neurocomputing,2014,138:209-217

[10]Lee D D,Seung H S.Algorithms for non-negative matrix factorization //Advance of Neural Information Process System[C].MIT Press,Cambridge,MA,USA,2001

[11]Li Xiangbao,Yang Yupu,Zhang Weidong.Fault detection method for non-Gaussian processes based on non-negative matrix factorization [J].Asia-Pac.J.Chem.Eng,2013,8:362-370

[12]Liu Weixiang,Zheng Nanning,Lu Xiaofeng.Non-negative matrix factorization for visual coding//Proceedings of 2003 IEEE International Conference on Acoustics,Speech,and Signal Processing[C].2003,3:293-296

[13]Gonzalez E,Zhang Y.Accelerating the Lee-Seung algorithm for nonnegative matrix factorization[R].Technical Report TR-05-02,Rice University,2005

[14]Hoyer P.Non-negative matrix factorization with sparseness constraints [J].Journal of Machine Learning Research,2004,5:1457-1469

[15]Olshausen B A,Field D J.Emergence of simple-cell receptive field properties by learning a sparse code for natural images [J].Nature,1996,381:607-609

[16]Albright R,Cox J,Duling D,Langville A N,Meyer C D.Algorithms,initializations,and convergence for the nonnegative matrix factorization[R].NCSU Technical Report Math 81706,North Caroline State University,2006

[17]Chen Q,Wynne R J,Goulding P,Sandoz D.The application of principal component analysis and kernel density estimation to enhance process monitoring [J].Control Engineering Practice,2000,8:531-543

[18]Botev Z I,Kroese D P.The generalized cross entropy method,with applications to probability density estimation [J].Methodology and Computing in Applied Probability,2011,13:1-27

[19]Downs J J,Vogel E.A plant-wide industrial process control problem [J].Computers and Chemical Engineering,1993,17:245-255

[20]Yu Jianbo.Local and global principal component analysis for process monitoring [J].Journal of Process Control,2012,22:1358-1373