朱雄卓，張瀚文，楊春節(jié)

（浙江大學控制科學與工程學院，浙江杭州310027）

引 言

鋼鐵工業(yè)是我國現代工業(yè)和國民經濟重要支柱產業(yè)，也是國民經濟密不可分的重要組成部分。而高爐煉鐵是整個鋼鐵工業(yè)的核心環(huán)節(jié)，也是鋼鐵制造過程中能質流轉換的關鍵工序，它的能耗占鋼鐵生產總能耗的比例高達70%。實現高爐煉鐵過程安全運行，是鋼鐵工業(yè)深度節(jié)能減排和提質增效的重中之重。在高爐煉鐵過程中，由于原料品質波動、人員的誤操作、設備故障等原因，高爐爐況異常時有發(fā)生。高爐爐況一旦發(fā)生異常，往往會引起燃料比增加、休風檢修時間增長、鐵水質量不達標等問題，不僅會造成資源和設備的重大損失，降低高爐爐齡，甚至有可能引發(fā)事故造成人員傷亡。因此，提高高爐煉鐵過程的安全性，尤其是減少高爐事故和異常爐況出現的頻率，對減少能源消耗、保證設備和人員的安全以及提高鋼鐵生產流程的經濟效益有著十分重要的意義。而一個準確的異常監(jiān)測模型可以及時地為操作者提供報警，提前對爐況進行調整，避免危險情況的發(fā)生[1-4]。

作為鋼鐵生產的重要環(huán)節(jié)，高爐煉鐵主要是通過內部一系列復雜的物理化學反應生產鐵水。高爐是主要的反應堆，是整個煉鐵過程的核心。高爐的輸入和輸出通?？梢苑譃閮刹糠?。第一部分由含鐵材料（例如鐵礦石、燒結礦、球團礦）、焦炭和助熔劑組成，它們從頂部進入高爐。輸入的第二部分由干燥的熱空氣、富氧燃料（例如焦油或煤粉）和濕氣組成，這些氣體從爐底吹入爐中。同樣，作為輸出的一部分，液態(tài)鐵水和爐渣從爐底流出，而煤氣則從爐頂收集[5-6]。然而，大型高爐數據有著數據非線性、非高斯分布和時變等特點，這使得高爐系統故障監(jiān)測與診斷成為極具挑戰(zhàn)性的課題，成為當今世界冶金科技研究的前沿熱點、難點[7-12]。

本文聚焦高爐數據的非高斯性，而針對一般的數據非高斯分布的監(jiān)測問題，傳統的基于PCA 和PLS的監(jiān)測方法可能不能很好地發(fā)揮作用。在構建過程監(jiān)控的T2和SPE 統計量時，需要在假設潛在變量為高斯分布的前提下計算其控制限。為了解決這個問題，Li 等[13-14]提出了ICA 方法，通過尋找統計上獨立且非高斯分布的隱成分完成故障監(jiān)測問題；后來Kano 等[15]又開發(fā)了基于統計學習的過程監(jiān)測的統一框架，將基于PCA 和基于ICA 的監(jiān)測方法結合起來；隨著算法的發(fā)展，Ge 等[16]提出了一種基于ICA-PCA 的兩步監(jiān)測方法來進行故障檢測和識別。與此同時，也有其他學者Thissen 等[17]、Yu 等[18]、Chen等[19]采用高斯混合模型對非高斯監(jiān)測問題進行解決。

目前高爐異常爐況的監(jiān)測方法主要可分為兩類：基于專家系統的方法和數據驅動的方法。其中，數據驅動的方法又可以分為基于機器學習的方法和基于多元統計分析的方法，如圖1所示。

早在1986 年，日本的NKK 公司就將專家系統引入了高爐，日本川崎水島廠1987 年開發(fā)了Go-Stop系統，通過采用模糊控制的思想，基于專家經驗構建規(guī)則庫，結合復合參數對高爐運行狀況進行綜合判斷[20]。李芳[21]基于統計學和模糊數學，提出了高爐故障的特征參數提取方法，并給出了異常爐況的診斷方法。

圖1 常見高爐故障監(jiān)測方法Fig.1 General blast furnace fault monitoring methods

而基于機器學習的高爐煉鐵過程監(jiān)測方法主要是利用高爐正常爐況以及各種異常下的樣本數據訓練模型，從而達到實時監(jiān)測的目的。在高爐的異常爐況監(jiān)測和診斷方法中，支持向量機(SVM)和人工神經網絡(ANN)取得了較多的研究成果。Tian等[22]在2010 年提出了一種基于Bagging 思想的支持向量機集成框架，利用此框架提出了一種新的異常爐況診斷方法；為了提高診斷的速度，Liu 等[23]提出了一種基于成本意識的最小二乘支持向量機多分類方法用于高爐的異常診斷；趙明[24]基于BP 人工神經網絡建立了異常爐況的預報模型。

但是由于高爐系統復雜，知識積累緩慢，專家系統只能針對性地解決部分問題；而機器學習方法因為存在歷史故障數據難以獲得、標簽不一定準確等問題，也難取得很好的效果。相比之下基于統計分析的方法有兩個更大的優(yōu)勢：

（1）和那些基于模型的方法相比，基于數據驅動的方法需要更少的機理知識和過程的因果關系。

（2）與專家系統、SVM 和其他機器學習方法相比，它不需要故障的先驗信息，對故障樣本也沒有要求。

因此近些年來基于統計分析的高爐異常監(jiān)測方法得到了更多的發(fā)展。其中，Vanhatalo[25]針對高爐的多變量過程監(jiān)測問題，利用投影實現了多變量監(jiān)測信息的降維，進行了主元模型在實際過程中的測試并設計了在線算法；Zhang 等為了解決高爐運行過程中的熱風爐切換問題提出來一種兩階段PCA 的監(jiān)測方法[26]，同時將該方法拓展到存在多種故障的多座高爐數據集上[27]；而進一步考慮到高爐數據呈現的非高斯特性和時變特性，Zhou 等[28]采用凸包代替了傳統需要滿足高斯分布的T2統計量，提出了一種基于滑窗主元凸包的PCA 算法，實現在線的高爐過程監(jiān)測；為了解決同樣非高斯的問題，孫夢園[29]提出了一種基于改進ICA 算法的異常爐況監(jiān)測方法，通過結合多元指數加權平均的方法，適應了高爐的時變特性，也一定程度上提高了監(jiān)測的準確性。隨后，Shang 等提出了一種基于主趨勢邏輯回歸的方法用于異常爐況的監(jiān)測和分離[30]，同時還針對高爐故障診斷提出了一種遞推變元統計分析方法，并采用了一種指標切換策略剔除了熱風爐切換帶來的影響[31]。

1 高爐煉鐵過程數據分析

1.1 高爐煉鐵數據簡介

本次實驗數據采集于華南某大型鋼鐵集團的2650 m3高爐。存儲于數據庫中的數據是反映高爐運行狀態(tài)的關鍵性指標。爐體的數據平均10 s采樣一次，含有35個變量，如表1所示。

根據操作員的經驗，熱風壓力、全壓差、理論燃燒溫度和鼓風動能是反映高爐運行狀態(tài)的關鍵性指標。因此選擇2018 年7 月1 日一整天共7819 個樣本，其中這四個關鍵指標如圖2所示。

從圖中可以看出，有部分變量每間隔一段時間會出現一個峰狀的擾動，這是由于高爐爐體熱風爐切換導致的，高爐往往有不止一個熱風爐，當一個熱風爐在給高爐加熱時，另外的熱風爐需要進行蓄熱，然后在一段時間后再進行切換，以保證高爐不間斷運行。

在對煉鐵廠實地進行調研和閱讀多篇文獻后，總結目前高爐異常爐況監(jiān)測面臨的問題主要有以下兩點：

表1 數據集的變量列表Table 1 A list of variables for the dataset

圖2 關鍵指標一天內變化情況Fig.2 General blast furnace fault monitoring methods

（1）高爐狀態(tài)和數據具有時變特性。

（2）高爐數據具有主要由熱風爐切換的峰狀干擾帶來非高斯分布的特點。

為了進一步說明上述存在的問題，利用實際數據進行驗證。

1.2 高爐煉鐵數據簡介

為了檢驗數據樣本的非高斯特性，本次測試選取了2018年7月1日一整天共7819個樣本的4個典型變量，繪制其數據分布直方圖如圖3所示。

從4個典型變量熱風壓力、全壓差、理論燃燒溫度和鼓風動能的直方圖分布粗略來看，變量并不符合高斯分布的要求，為了進一步對問題進行說明，對該四個變量采用KS（Kolmogorov-Smirnov）檢驗中的正態(tài)分布檢驗方法進行檢驗，結果如表2 所示。事實上，高爐樣本的非高斯分布很大程度上就是由于熱風爐切換所導致的。

表2 KS正態(tài)分布檢驗Table 2 A list of variables for the dataset

1.3 高爐煉鐵數據時變特點

為了研究高爐過程數據的時變特性，選取2018年7 月1 日至2018 年7 月20 日時隔20 天的兩組數據，樣本數分別為7819 個和8467 個，針對這兩組數據的四個典型變量計算均值、標準差，相關數據如表3所示。

表3 兩組數據均值與標準差Table 3 Data mean and standard deviation of the twogroups

PCA 降維分解結果和數據的均值、標準差密切相關。從表3 中可以看出，僅僅20 天，變量的均值和標準差變化就十分明顯，如果在PCA 建模過程中，不對模型進行更新，舊模型顯然無法適應新數據的工作狀態(tài)，會與真實狀態(tài)產生較大的偏差。因此，針對數據具有的時變特性，模型必須也進行相應的處理。

2 高爐煉鐵異常監(jiān)測算法

針對上述存在的2 個難點，本文提出了一種基于高斯混合模型的MWPCA 異常監(jiān)測算法，通過利用高斯混合模型的方法將基于高斯分布假設的T2統計量拓展到非高斯分布的數據集上，然后加入滑窗機制使模型可以適應數據的時變特性。

2.1 PCA監(jiān)測算法簡介

圖3 關鍵指標分布直方圖Fig.3 Distribution histogram of key index

主成分分析（PCA，principal component analysis）是一種以保留數據信息量為目的的降維算法，由Pearson[32]在1901 年提出。隨后被Hotelling[33]和Jackson[34]應用于多元統計領域，然后逐漸被引入到過程監(jiān)測領域。

PCA 在執(zhí)行前要先對樣本進行標準化，假設標準化后的數據為X(N×M)，其中N 為樣本的數目，M為變量數目，需要先求取協方差矩陣：

通過對協方差矩陣使用SVD 分解可以得到特征 值(λ1,λ2,…,λA)，以 及 其 對 應 的 特 征 向 量(p1,p2,…,pA)，其中A(A＜M)為保留主元的個數，通常由累計方差貢獻率方法得到。在得到特征向量后，主成分可以表示為：

其中P(N×A)為由特征向量組成的負載矩陣，單獨一個主成分可以表示為：

當建立好PCA 模型后，通常用T2統計量和SPE統計量來監(jiān)測過程是否發(fā)生了異常，其中T2統計量的計算方式為：

其中，t為主成分向量，S = diag(λ1,…,λA)為前A個主成分所對應的特征值的對角矩陣。

當數據樣本近似高斯分布時，T2統計量可視為服從F分布，其控制限可以根據式（5）計算：

其中，N 為樣本的數目，A 為保留的主成分的個數，α 為置信度，Fα(A,N-A)為在置信度為α、自由度為A和N-A情況下F分布的臨界值。

SPE統計量計算方式：

其中，e為重構x后的殘差。其控制限為：

2.2 高斯混合模型簡介

高斯混合模型（Gaussian mixture model,GMM），是一種被廣泛應用的聚類算法，它可以將數據聚類成若干個基于高斯概率密度函數的成分，因此它幾乎可以擬合逼近各種分布[35-36]。

對于數據集X={x1,x2,…,xN}中的一個采樣xi，其高斯混合分布的密度函數為：

高斯混合模型的目的就是估計模型參數K、πk、θk=(μk,σ)，但是參數K 一般提前給定，因此需要求解πk和θk=(μk,σ)。傳統的參數估計往往使用極大似然估計，但對于高斯混合模型，在取對數之后難以對其進行求偏導，因此求解高斯混合模型需要使用迭代求解的EM（expectation maximization）算法。

求解高斯混合模型的主要步驟如圖4所示。

圖4 GMM模型迭代流程圖Fig.4 Iteration flow chart of GMM model

即確定所有高斯密度函數之后，可以計算出現xi的概率，以及第k個高斯成分的貢獻率。在計算完所有的zik后，再反向更新每一個高斯成分：

如此反復迭代，直到滿足收斂條件為止。最終可以得到每個高斯成分的權重，以及高斯成分的均值、方差和協方差。

通常高斯混合模型的高斯元個數可以根據AIC和BIC指標進行選取，但是由于高爐數據維度高，數據量大，指標使用的效果并不理想。因此在本文中高斯元個數是根據離線訓練時的使用效果以及程序所需的運行時間進行綜合考慮進行選取的。一般k 的選取值可以在取得較為理想的效果后，最終選為k+1或k+2，以適應高爐可能存在的新分布。

2.3 基于高斯混合模型的T2統計量

由于熱風爐切換等原因的存在，導致數據具有非高斯分布的特點，傳統的T2統計量在對高爐過程進行監(jiān)測時會有誤報的情況產生，在先前的研究當中，學者們往往選擇區(qū)分正常爐況和熱風爐切換的情況，再進行監(jiān)測。在本文中，提出了一種基于高斯混合模型的T2統計量改進方法，使PCA 模型可以在不區(qū)分是否為熱風爐切換的情況下進行監(jiān)測。

選取高爐數據中2018 年10 月1 日至2018 年10月10日10天的80000個數據樣本，為了更加直觀地感受樣本分布和T2統計量范圍，將前兩個主成分在2 維空間進行展現，并繪制了置信度為99.99%、99.9%和99%下T2統計量覆蓋的范圍，置信度越高T2統計量包括的范圍就越大，如圖5所示。

圖5 樣本與T2統計量覆蓋范圍Fig.5 Samples and the coverage of T2 statistics

由于是2 維空間，相比于真正用于監(jiān)測的8 維空間，必然缺失了比較多的信息，但是仍然可以發(fā)現數據的分布是有偏向性的，而不是真正符合高斯分布，T2統計量對樣本分布的擬合是存在比較嚴重的問題的?；谶@個問題以及高斯混合模型強大的分布模擬能力，對T2統計量進行了改進。

首先對原始數據降維后的主成分使用高斯混合模型進行求解，此時選擇的高斯元個數k 為6，將高斯混合模型解算的結果從2 維空間進行展現，如圖6所示。

其中，每一個點代表一個樣本，不同顏色的點代表他們歸屬于不同高斯成分，除去樣本點和代表不同置信區(qū)間下T2統計量范圍的虛線，剩下橢圓形的淺灰色陰影部分為高斯混合模型得到的覆蓋范圍，相比于T2統計量覆蓋的范圍，高斯混合模型的結果顯然更加準確一些，不僅基本上覆蓋了所有樣本，而且一些沒有樣本的空間不被包含在其中。對于這種通過高斯混合模型劃分區(qū)域再計算是否超限的T2統計量改進方法，將其稱為GMM-T2統計量。

圖6 樣本和GMM-T2統計量覆蓋范圍Fig.6 Samples and the coverage of GMM-T2 statistics

2.4 基于GMM-T2統計量的MWPCA算法

在實際監(jiān)測過程中，基于GMM-T2統計量的MWPCA 監(jiān)測算法可以分為離線建模、在線監(jiān)測和模型更新三部分，具體如下所示。

2.4.1 離線建模

（1）收集監(jiān)測過程數據，將其進行標準化后的數據記為Xtrain(N × m)，樣本數為N，變量數為m。計算數據的協方差矩陣Rtrain，利用SVD 分解得到特征向量和特征值，然后根據累計方差貢獻率選取主成分個數A，進而得到負載矩陣Ptrain以及得分矩陣Ttrain，其中Ttrain為Xtrain在Ptrain方向上的投影。

（2）確定高斯混合模型的高斯成分個數為k，對使用PCA 得到的主成分使用高斯混合模型進行求解，得到每個高斯成分的均值μi和協方差Σi，并求出每個樣本的GMM-T2統計量以及每個高斯成分的GMM-T2統計量上限GMM - T2UCL：

其中，0 ≤i ≤k，即第i 個高斯成分，GMM - T2ij為隸屬于第i個高斯成分的第j個樣本的GMM-T2統計量，tij為隸屬于第i 個高斯成分的第j 個樣本的主成分，bi為第i 個高斯成分包含的樣本數，每個高斯成分的GMM-T2統計量上限即為該成分包含的樣本中，GMM-T2統計量的最大值。

同時也計算出SPE統計量的上限值SPEUCL。

2.4.2 在線監(jiān)測 當有新樣本時，同樣進行標準化，將標準化后的數據記為Xtest，然后計算其主成分t，利用高斯混合模型判斷其歸屬于第i個高斯成分，然后計算其對應的GMM - T2統計量和SPE 統計量：

其中，μi和Σi分別是離線訓練時第i個高斯成分均值和協方差。將計算出來的統計量與GMM -和SPEUCL進行比較，如果有一個指標超限，則可判斷系統出現異常；若兩個指標均未超限，則系統處于正常運行狀態(tài)。

2.4.3 模型更新 當算法判斷樣本為正常樣本時，就將該樣本加入訓練集中，然后將訓練集中最早的一個樣本剔除，保持訓練集數量不變，當判斷為正常樣本的新樣本數到達W 個時，后臺使用新的訓練集重新執(zhí)行離線建模，完成模型的更新后替代當前用于在線監(jiān)測的模型。

3 高爐數據集應用實例

3.1 算法克服非高斯性檢驗

為了檢驗基于GMM-T2統計量的MWPCA 算法的應用效果，選取2018 年11 月5 日至2018 年11 月15 日10 天的80000 個數據樣本作為訓練集；2018年11月26日的5000個樣本為測試集。其中測試集包含了一個故障，為了證明算法的有效性，將普通T2統計量和GMM-T2統計量進行對比，監(jiān)測結果如圖7所示。

如圖7 所示，對于故障，SPE 統計量、T2統計量和GMM - T2統計量表現類似，都發(fā)出了報警，其中SPE統計量在1360時刻報警，T2統計量在1359時刻報警，GMM - T2統計量在1360 時刻報警，GMM -T2統計量報警晚于T2統計量1個時刻。

其中T2統計量由于熱風爐切換的原因有很多誤報的情況出現，而GMM - T2統計量則沒有受到任何干擾。此外，T2統計量和GMM - T2統計量在3100 時刻左右有一個較大的分歧，即鼓風動能有一個微小的波動導致T2統計量超限了，通過觀察其他過程變量發(fā)現并沒有這個波動的出現，結合操作者的報告也并不存在故障現象，因此可以判斷這是一個正常的波動。針對這個正常波動，T2統計量產生了大規(guī)模誤報，而GMM - T2統計量則只有個別誤報。由此可見，GMM - T2統計量不僅可以克服熱風爐切換擾動，并且在一定程度上可以避免正常波動帶來的誤報，擁有更強的魯棒性。

圖7 高斯性檢驗的監(jiān)測結果Fig.7 Monitoring results of Gaussian test

3.2 算法克服時變性檢驗

為了測試基于GMM - T2統計量的MWPCA 算法對時變數據的適應能力，選取2018年10月1日至2018 年10 月10 日10 天的80000 個數據樣本作為起始訓練集；2018年11月5日的4000個樣本為最終測試集。未加入滑窗時的監(jiān)測效果如圖8 所示，加入滑窗后的監(jiān)測效果如圖9所示。

相比于圖8未加入滑窗時兩個統計量大幅度超限的異常情況，圖9 中當算法加入滑窗來更新樣本之后，SPE 統計量明顯不會出現因為模型未更新導致的大規(guī)模誤報情況，GMM - T2統計量不僅沒有大規(guī)模誤報，而且上限值相對變得更小了，這說明監(jiān)測范圍也變得更加準確了。

圖8 加入滑窗前的監(jiān)測結果Fig.8 Monitoring results of moving windows were not added

圖9 加入滑窗后的監(jiān)測結果Fig.9 Monitoring results after moving windows were added

3.3 算法對比

由于GMM - T2統計量是PCA 監(jiān)測算法傳統T2統計量的改進，而且其在判斷統計量的控制限時，取值為隸屬于當前高斯成分的所有樣本的GMM -T2統計量的最大值，這與傳統T2統計量根據置信度決定控制限是有區(qū)別的。因此為了對比更加全面，也加入了傳統T2統計量不根據置信度，而直接取T2統計量最大值作為控制限時的監(jiān)測效果，這里簡稱為MAX-T2統計量，很顯然相比于主元凸包[28]簡單將樣本點相連包圍的形式，MAX-T2所覆蓋的正常范圍更廣。此外，也比較了GMM - T2統計量與專門克服熱風爐擾動的兩階段PCA 算法和常用于克服非線性的KPCA 算法，可以進一步表明GMM - T2統計量在克服非高斯、非線性和熱風爐擾動上的獨特優(yōu)勢。其中，統計量選取的置信區(qū)間的置信度均為99.99%。

故障監(jiān)測的監(jiān)測效果主要體現于監(jiān)測的誤報率、故障的報警起始時刻和報警持續(xù)時間（可以判斷是否有漏報的情況產生）上，因此下面的對比主要表現在以上3個指標里。

對于3.1節(jié)中的測試集，結果對比如表4所示。

表4 不同算法對比Table 4 Comparison of different algorithms

其中時刻指的是樣本數，即第1859 時刻為第1859 個樣本，持續(xù)186 個時刻即持續(xù)186 個樣本的時間，其中每個樣本的采樣間隔為10 s，由于操作員的交班記錄里并不會明確給出故障的具體時間，因此希望通過故障報警的持續(xù)時間從側面看出算法報警的完整性。

對于該測試集，每個統計量的表現都不太一樣，對于T2統計量，很明顯存在誤報的情況出現，因此需要進行改進；但是如果直接簡單地利用MAXT2統計量，可以看出其效果是不如GMM - T2統計量的，由于MAX-T2統計量是直接取了最大T2統計量值作為控制限，因此其覆蓋的范圍很大且沒有針對性，這從其報警時刻的延后和持續(xù)時間的縮短上都可以有所體現；而主元凸包認定的正常樣本范圍是直接由樣本連接得到的，這個范圍會小于MAX-T2統計量的范圍，但是準確度會比一般的T2統計量好，所以它的誤報率介于兩者之間，具體數值與樣本的預處理息息相關；兩階段PCA 的T2統計量由于其有一個報警持續(xù)時間限W，只有當報警持續(xù)時間超過W 時才會視為發(fā)生故障，因此報警時刻的延后和持續(xù)時間的縮短是不可避免的。至于KPCA 監(jiān)測算法由于處理了樣本的非線性，使得誤報率有所降低，但是效果并不是特別理想。而GMM - T2統計量不僅在誤報率上得到了大幅降低，而且報警時刻也僅僅只有1 個時刻的落后，由于高斯混合模型更加具有針對性的擬合方式，也避免了MAX-T2統計量簡單粗暴方式所帶來的覆蓋范圍過大問題，同時相比于專門針對熱風爐切換的兩階段PCA 方法，GMM - T2統計量在報警時刻和持續(xù)時間上都有著無可比擬的優(yōu)勢。

4 結 論

異常監(jiān)測是高爐煉鐵過程中的一項重要任務，可靠的早期檢測可以為操作員提供足夠的時間來采取控制措施，以使過程恢復正常并防止可能發(fā)生的事故。在本文中，對傳統的基于PCA 的監(jiān)測方法進行了兩方面的改進，以應對過程數據的非高斯性和時變特性，追求誤報率和故障報警時間上的平衡。其中，基于GMM - T2統計量的MWPCA 算法，在實際應用的過程中表現優(yōu)異，對高爐異常監(jiān)測帶來了顯著的幫助。