杜 建 張 磊 李繼庚 洪蒙納，3 滿 奕，4，*

（1. 華南理工大學制漿造紙工程國家重點實驗室，廣東廣州，510640；2. 廣東省節(jié)能中心，廣東廣州，510030；3. 中新國際聯(lián)合研究院，廣東廣州，510555；4. 人工智能與數(shù)字經(jīng)濟廣東省實驗室（廣州），廣東廣州，510335）

隨著過程工業(yè)生產(chǎn)市場競爭日益激烈，提高生產(chǎn)效率，降低生產(chǎn)損耗，是過程工業(yè)生產(chǎn)企業(yè)亟待解決的問題。其中，生產(chǎn)事故是造成企業(yè)產(chǎn)量及生產(chǎn)效率不能有效提高的重要原因之一。造紙工業(yè)是一種市場前景廣闊的過程工業(yè)［1-3］，但受漿料種類、設備故障、工藝環(huán)境等眾多不穩(wěn)定性因素的影響，造紙機在生產(chǎn)過程中易發(fā)生紙幅斷裂的故障。Lahti 等［4］通過拉力測試定量地分析出定量、纖維取向、緊度和厚度等紙張性能變化對紙張纖維結(jié)合強度的影響。Krasnoshlyk等［5］使用X 射線成像技術成功分析出低密度紙張（LDP）中裂縫擴展引起的紙張結(jié)構(gòu)變化。然而這些方法只適用于離線檢測，無法應用于實時造紙生產(chǎn)過程中。紙廠對斷紙故障的提前預防大多采用定期巡查的方法，但這樣易導致維修過剩，造成不必要的能源資源浪費。

近年來，利用傳感器監(jiān)測的狀態(tài)參數(shù)數(shù)據(jù)對設備運行狀況進行監(jiān)測逐漸成為研究的熱點［6-8］，即構(gòu)建設備健康運行階段數(shù)據(jù)集。實時運行數(shù)據(jù)狀態(tài)與健康基準數(shù)據(jù)集之間的差值為故障信息，當故障信息超過某個設定閾值時，模型發(fā)出預警信號。Wang 等［9］通過計算軸承振動信號與已知健康狀態(tài)的多個統(tǒng)計偏差來預估設備故障信息，從而了解軸承的健康狀態(tài)。Zhang等［10］提出一種基于運行工況識別技術和高斯混合模型（GMM）的在線健康狀態(tài)評估方法，該方法對風力發(fā)電機組的健康狀態(tài)評估效果良好。然而，在選取狀態(tài)特征變量時，上述幾種方法高度依賴人工經(jīng)驗，但過程工業(yè)故障與現(xiàn)場實際工藝環(huán)境相關性較大，只依賴人工經(jīng)驗無法提取與過程工業(yè)故障相關性強的特征變量。

針對以上方法在過程工業(yè)生產(chǎn)應用中存在的不足，筆者提出了一種基于GMM-馬氏距離（MD）組合算法的過程工業(yè)故障預測模型。GMM 是一種無監(jiān)督學習方式的數(shù)據(jù)聚類算法，常被用于識別最能代表特征空間數(shù)據(jù)的子集［11］。MD是一種被廣泛應用于健康狀態(tài)評估和故障檢測領域的距離度量方法［12］，可用于度量空間中數(shù)據(jù)點與數(shù)據(jù)分布之間基于樣本分布的距離。該過程工業(yè)故障預測模型具有如下特點：

（1）通過滑動窗口以提取時序數(shù)據(jù)的平均值和方差特征，之后對所有特征進行聚類分析，判斷異常類別數(shù)據(jù)，并借助相關性分析篩選與異常類別數(shù)據(jù)有關聯(lián)的特征；

（2）針對不同參數(shù)與斷紙故障的關聯(lián)關系和耦合現(xiàn)象，基于GMM 和MD 的組合算法構(gòu)建健康指標，充分考慮各參數(shù)之間的相關性。

1 過程工業(yè)故障預測模型

本研究所提出的過程工業(yè)故障預測模型主要包括數(shù)據(jù)準備、特征工程、距離度量和健康評估4 個部分，具體流程圖如圖1所示。

圖1 過程工業(yè)故障預測流程圖Fig.1 Flow chart of process industry fault prediction

第1階段，分析狀態(tài)變量間的相關性，去除冗余變量與無關變量并完成輸入變量選擇，然后針對數(shù)據(jù)采集系統(tǒng)中存在的一些問題對原始數(shù)據(jù)進行清洗；第2階段，通過滑動窗口處理保留的變量特征，然后從平均值和方差2個維度提取特征信息，并利用基于聚類的方法對特征信息進行合理篩選，以獲得核心特征變量；第3階段，依據(jù)不同的故障原因?qū)⒑诵奶卣髯兞縿澐譃椴煌淖犹卣髯兞考瑢γ總€子特征變量集使用GMM-MD 組合算法以生成健康指標，從而實現(xiàn)多種參數(shù)實時在線監(jiān)測；第4階段，利用核密度估計確定健康指標的閾值，以判斷參數(shù)的故障信息和評估生產(chǎn)過程的健康狀況。本研究利用國內(nèi)某造紙廠實時生產(chǎn)數(shù)據(jù)對該過程工業(yè)故障預測模型進行驗證。

1.1 數(shù)據(jù)準備

對生產(chǎn)過程中的故障進行研究，首先需要從工業(yè)實時采集系統(tǒng)中選取在生產(chǎn)過程中與生產(chǎn)故障相關且影響較大的狀態(tài)特征參數(shù)，構(gòu)建狀態(tài)特征向量。相關系數(shù)是目前較為常用的一種相關性分析方法，即假設存在變量x和y，則它們之間的相關系數(shù)rxy可通過式（1）計算獲得［13］：

式中，n是樣本的數(shù)量，xi和yi分別是x和y第i個值，分別是x和y的平均值。r值表示 2 個變

xy量之間的相關程度，取值在［-1，1］之間，rxy的絕對值越接近1，2個變量的線性相關性越強。

1.2 特征工程

1.2.1 特征提取

多變量時間序列是按多個變量的時間序列順序形成的一組數(shù)據(jù)，可視為按照時間索引的多個單變量序列的組合，蘊含了故障的發(fā)生規(guī)律和趨勢，在數(shù)學上可表示為二維矩陣。本研究利用滑動窗口截取時間序列數(shù)據(jù)，將原始數(shù)據(jù)截取成許多按照時間索引的數(shù)據(jù)片段。為挖掘出數(shù)據(jù)片段的潛在故障信息，需要將原始數(shù)據(jù)的不穩(wěn)定特性對生產(chǎn)過程健康狀態(tài)評估的影響降到最低，加強關鍵特征對評估的影響力［14］。

1.2.2 特征加強

本研究主要從集中趨勢和離散趨勢2個角度分析故障信息，利用滑動窗口截取數(shù)據(jù)片段后，對窗口內(nèi)每個變量分別提取平均值和方差2 個數(shù)學統(tǒng)計特征，這樣原來n維特征變成2n維特征。為更加準確提取與生產(chǎn)故障有相關性的特征，本研究提出一種基于聚類方法對特征信息進行合理篩選以獲得核心特征變量的方法。首先，對2n維特征變量分別使用K-Means聚類算法，判斷出異常類別；然后，提取生產(chǎn)故障前一段時間內(nèi)的數(shù)據(jù)，對出現(xiàn)少數(shù)類別的數(shù)據(jù)點進行標記；最后，利用點二列相關性系數(shù)對提取的數(shù)據(jù)與新生成的標簽進行相關性分析。

1.2.2.1 聚類數(shù)判斷

使用聚類算法前，首先需給定聚類數(shù)k值。為了確定最佳的k值，可利用聚類簇內(nèi)的稠密程度和聚類簇間的離散程度來評估聚類效果。本研究采用手肘法則和輪廓系數(shù)來確定最佳的k值。

手肘法則的核心思想是：隨著k值的增大，樣本劃分更加精細，每個簇的聚合程度逐漸提高，則誤差平方和（SSE）逐漸下降。并且，當k值小于真實最佳聚類數(shù)時，k值的增大會大幅提高每個簇的聚合程度，故SSE的下降幅度很大；而當k值達到真實最佳聚類數(shù)時，再增大k值所得到的聚合程度回報迅速變小，因而SSE的下降幅度驟減；之后隨著k值的繼續(xù)增大，SSE趨于平緩，即SSE和k值的關系圖呈手肘形狀，此肘部對應的k值就是數(shù)據(jù)的真實最佳聚類數(shù)［15］。SSE可通過式（2）計算得到：

式中，p代表第i個類組Li中的數(shù)據(jù)對象，qi代表第i個類組中所有數(shù)據(jù)對象的平均值。

聚類完成后，類別內(nèi)部數(shù)據(jù)的協(xié)方差越小，類別之間的協(xié)方差越大，即類本身越緊密，類與類之間越分散，聚類效果越好，相應地輪廓系數(shù)（S（i））的值介于［-1，1］，其絕對值越趨近于1，聚類簇內(nèi)的稠密程度和聚類簇間的離散程度越優(yōu)。S（i）可按式（3）表示［16］：

式中，a(i)表示i向量到同一簇內(nèi)其他點不相似程度的平均值，b(i)表示i向量到其他簇的平均不相似程度的最小值。

1.2.2.2 相關性分析

點二列相關法[17]是可用于測量連續(xù)變量和二分類變量之間關聯(lián)程度的統(tǒng)計方法，二列相關系數(shù)（r）的取值在[-1，1]之間，其絕對值越接近1，變量間相關程度越高。r的計算如式（4）所示：

式中，二分類變量y的取值可用0 和1 來表示，并設取1的比率為p，取0的比率為q；表示連續(xù)變量中與二分類變量p對應數(shù)據(jù)的平均值，表示連續(xù)變量中與二分類變量q對應數(shù)據(jù)的平均值，Sx表示連續(xù)變量的標準差。

1.3 算法描述

基于GMM-MD 組合算法的過程工業(yè)故障預測模型的流程圖如圖2 所示。該模型主要利用GMM 構(gòu)建健康基準數(shù)據(jù)集并利用最大期望（EM）算法進行參數(shù)估計［18］，然后利用MD生成健康指標，以評估紙機生產(chǎn)過程的健康程度。

圖2 GMM-MD組合算法流程圖Fig.2 Flow chart of the GMM-MD combinational algorithm

GMM-MD組合算法的主要步驟如下［18］。

輸入：樣本集D={x1，x2，…，xm}；高斯混合成分個數(shù)G。

Step1.初始化高斯混合分布的模型參數(shù){(αi，μi，Σi)∣1 ≤i≤G}，其中αi為選擇第i個混合成分的概率，μi是n維均值向量，Σi是n×n的協(xié)方差矩陣。

Step2.根據(jù)貝葉斯定理生成式（5），并由此計算xj由各混合成分生成的后驗概率 （γj）i，即γji=pM(zj=i∣xj)(1 ≤i≤G) （EM 算 法 的 E 步 ）， 其 中 ，zj∈{1，2，…，G}表示生成樣本xj的高斯混合成分，P表示樣本數(shù)據(jù)的概率密度函數(shù)。

Step4.當EM 算法滿足停止條件（例如已達到最大迭代輪數(shù)，或者如式（9）所示的似然函數(shù)LL(D)增長很少，甚至不再增長）時，根據(jù)及高斯混合分布確定簇劃分。

Step5.輸出：簇劃分。

Step6.生成健康基準數(shù)據(jù)集。

Step7.計算采樣點與健康基準數(shù)據(jù)集之間的MD。

式中，μ和S分別為健康基準數(shù)據(jù)集的均值和協(xié)方差矩陣，x為采樣點的值。

Step8.本研究的健康指標（HI）可按式（11）計算：

1.4 健康評估

1.4.1 確定閾值

核密度估計是概率論用來估計隨機變量的概率密度函數(shù)，屬于非參數(shù)檢驗方法之一，是一種從數(shù)據(jù)樣本本身出發(fā)研究數(shù)據(jù)分布特征的方法［19］。本研究采用核密度估計確定報警閾值。假設x1，x2，…，xi，…，xN為獨立同分布的N個樣本點，則概率密度函數(shù)如式（12）所示：

式中，h是一個超參數(shù)（稱為帶寬），N表示樣本總數(shù)，Kh表示核函數(shù)。當帶寬大小固定后，基于不同核函數(shù)的估計結(jié)果幾乎沒有差異［20］。高斯核函數(shù)由于具有良好的數(shù)學計算特性而被廣泛應用于實際問題中。因而，本研究中Kh選取高斯核函數(shù)。h大小采用經(jīng)驗法并按式（13）計算［21］，其中，?為樣本標準差。

1.4.2 評價指標

如式（14）和式（15）所示，為了評估基于GMM-MD組合算法的過程工業(yè)故障預測模型的效果，故障預測精準率（Pew）和召回率（Rew）被分別描述為：

2 案例分析

本研究以國內(nèi)某生活用紙企業(yè)楊克烘缸紙機作為研究對象，利用其2個月的實時過程數(shù)據(jù)（期間斷紙故障發(fā)生次數(shù)為330次）作為訓練數(shù)據(jù)，然后又提取了2 個月的數(shù)據(jù)（期間斷紙故障發(fā)生次數(shù)為299 次）對基于GMM-MD 組合算法的過程工業(yè)故障預測模型進行驗證，采樣間隔為30 s，共351362組數(shù)據(jù)。經(jīng)相關系數(shù)分析后，選取卷紙缸速度、斷紙刮刀加載壓力、上漿泵出口壓力等16 個特征參數(shù)。利用滑動窗口截取數(shù)據(jù)片段，帶寬大小為10 min 內(nèi)20 個數(shù)據(jù)，步長為1。經(jīng)過提取變量平均值和方差2 個特征后，原來16維特征變成32維特征。

利用K-Means聚類算法分別對所有特征進行聚類分析，最佳聚類數(shù)由肘部法則和計算不同k值下輪廓系數(shù)判斷。以上漿泵出口壓力均值為例（見圖3），當k值小于 3 時，SSE隨著k值增大驟降；當k值大于 5時，SSE下降幅度隨著k值的繼續(xù)增大趨于平緩；而當k值為3～5 時，輪廓系數(shù)在k值為4 時最大。綜合考慮，將上漿泵出口壓力均值聚類成四類時效果最佳。

圖3 k值及對應的SSE和輪廓系數(shù)Fig.3 k Values and the corresponding SSE and silhouette coefficients

過程工業(yè)生產(chǎn)過程中變量數(shù)據(jù)波動與生產(chǎn)故障間存在滯后性，并且由于過程工業(yè)生產(chǎn)過程的敏感性和過程復雜性，變量數(shù)據(jù)波動與生產(chǎn)故障間的滯后時間通常沒有明顯的對應關系，因此無法有效、準確地標記異常波動數(shù)據(jù)。鑒于上述情況，本研究提取出斷紙故障前1 h 的數(shù)據(jù)，對出現(xiàn)異常類別的時間點進行標記，再進行相關性分析。認為當k值≤5時，數(shù)量比例低于10%的類別被視為異常類別。

由相關性分析得到的相關性結(jié)果可知，與斷紙故障相關性最大的是上漿泵出口壓力均值、抄前池液位均值、清潔刮刀加載壓力方差、起皺刮刀加載壓力方差、斷紙刮刀加載壓力方差這5個核心特征變量。這5個核心特征變量出現(xiàn)異常均可能導致斷紙故障的發(fā)生，因此依據(jù)斷紙原因?qū)⑦@5個核心特征變量劃分成2個子特征變量集，供后續(xù)研究分析。

分別對子特征變量集進行歸一化處理，然后使用GMM-MD組合算法得到相應的健康指標。本研究利用MD 的數(shù)據(jù)分布確定健康指標的報警閾值，其中，由上漿泵出口壓力均值和抄前池液位均值生成MD1，對應的健康指標為HI1；清潔刮刀加載壓力方差、斷紙刮刀加載壓力方差和起皺刮刀加載壓力方差3個特征生成MD2，對應的健康指標為HI2。

圖4為MD1和MD2的概率分布和累積概率分布曲線。根據(jù)90%的置信度可分別設定故障判斷閾值為1.921 和 4.632，由此計算HI1 閾值為 0.419，HI2 閾值為0.317。測試結(jié)果表明，GMM-MD 組合算法對斷紙故障的Pew為76.82%，Rew為72.50%，可較好地跟蹤造紙過程設備運行狀態(tài)的變化過程，起到過程工業(yè)故障提前預警作用。

圖4 概率分布及累積概率分布曲線Fig.4 Probability distribution and cumulative probability distribution curves

3 結(jié) 論

為了使過程工業(yè)生產(chǎn)企業(yè)能夠提高生產(chǎn)效率和降低生產(chǎn)成本，本研究介紹了一種基于高斯混合模型（GMM）和馬氏距離（MD）組合算法的過程工業(yè)故障預測模型，該模型可構(gòu)建健康指標以評估生產(chǎn)過程的健康狀況。利用國內(nèi)某造紙企業(yè)實時生產(chǎn)數(shù)據(jù)對模型進行了驗證；結(jié)果表明，模型的故障預測精準率為76.82%，召回率為72.50%，表明該模型可較好地跟蹤造紙生產(chǎn)過程中設備的運行狀態(tài)，起到過程工業(yè)故障提前預警作用。