基于因果分析的燒結生產狀態預測模型

李浩然，邱彤

（1 清華大學化學工程系，北京100084； 2 工業大數據系統與應用北京市重點實驗室，北京100084）

引 言

燒結單元是高爐煉鐵系統的重要一環，其生產水平的高低直接決定了供給后續高爐生產的燒結礦的質量[1-2]。經驗表明，生產的燒結礦中鐵含量（TFe）降低1%，高爐燃料比率增加2%，而高爐產量下降增加3%[3]，因此燒結單元生產的燒結礦質量與煉鐵過程的生產效益緊密相關，維護燒結機的健康工作狀態至關重要。通過在某鋼鐵企業調研過程發現，在實際生產過程中的燒結機狀態維護主要是通過手動調節來完成的，具有隨意性和滯后性的缺陷。為了改善這一問題，需要針對燒結過程建立模型，對生產狀態進行提前判斷，從而為現場操作工的操作提供可靠指導。

典型的燒結過程由配料、混合、進料、點火、卸料、破碎和篩選等步驟組成。首先各種鐵礦石、燒結返礦、焦炭和熔劑等根據設計方案配比混合，然后進料在移動小車上形成均勻燒結床層。接下來點火器點燃床的表面，在此期間移動小車下方的鼓風機通過抽風的方式在床層下方的風箱內產生負壓，由于床層內部含有焦粉和煤粉，點火之后床層在燃料、空氣和合適的溫度下繼續發生自燃，隨著臺車逐漸向燒結機尾端移動，燃燒不斷向下發展[4]。最終使得原料礦粉逐步形成具有一定粒度的燒結礦，作為鐵料進入到后續的高爐煉鐵生產過程中。

燒結過程具有以下兩個特點。

（1）時滯性。一方面，燒結系統中過程變量的變化與其對其下游變量產生影響之間存在時間間隔，稱為機制時滯。機制時滯隨著傳播路徑的不同而會有所差異，體現出多時間尺度特征。另一方面，對于生產過程中的同一批原料，由于不同變量測量的時間不同，因此會產生技術時滯。在針對同一批料的變量測量中，傳感器分別在t0、t1、t2和t3時刻測量了焦粉比例、一混加水量、點火溫度和終點溫度這四個變量，這四個時刻之間的不同時間間隔即為不同的技術時滯。這兩類時滯的存在使得燒結系統變量之間的影響呈現明顯的滯后。

（2）非線性。在燒結過程中涉及的化學反應和物理變化通常處于動態平衡中。在這種狀態下，燒結床內部發生一系列化學和物理變化的過程中能量和質量保持平衡，系統處于穩態[5]。燒結過程中除了發生焦炭燃燒、石灰石分解、白云石分解、鐵氧體系變化和含硫成分脫除等一系列化學變化外，還存在二維三相復雜傳熱傳質關系，因此燒結系統變量之間呈現明顯的非線性關系特征。

由于燒結過程時滯性和非線性這兩大特點，雖然現有一些利用機理模型對于燒結過程進行模擬的方法，如姜麗娟[6]建立了基于燒結過程的傳熱、傳質關系的二維模型，以此來獲得量熱法中漏風前的氣體溫度值準確率達91.3%。Zhou 等[7]提出了預測氣體濃度的方法，Cheng 等[8]提出了預測燒結礦的物理質量指標的能力的方法。但是由于燒結過程來料不穩定，生產模式也經常隨之變化，因此建立詳細的機理模型不僅工作量大，而且往往難以適應生產狀態變化。與此對應的是，近年來數據驅動方法得到了學界的關注和研究，如陳偉等[9]使用BP 神經網絡預測燒結礦化學組成，尹學等[10]使用GA-PSOBP 神經網絡預測燒結終點等，劉加達等[11]分別使用BPNN 和RNN 預測了燒結礦質量，易正明等[12]使用RBF-BP 混合神經網絡預測了燒結煙氣氮氧化物含量。數據驅動的方法將系統看作黑箱模型，在一些機理難以詳盡模擬的系統建模中取得了良好的效果。

因此本文采用結合機理分析和黑箱建模優點的方法，首先利用因果性分析的方法獲得燒結過程中的部分機理知識，即燒結變量之間影響的因果關系和作用時間窗口，然后將此知識應用到黑箱模型建模中，建立燒結系統生產狀態預測模型。

1 燒結過程因果性分析

1.1 建模變量選取

在燒結過程中，風箱負壓能夠維持空氣從料層表面向下穿過料層到達主風道的壓力差，使得料層內部發生充分燃燒，因此風箱負壓與穿過料層的氣流量、燃燒的程度、燒結的質量緊密相關。風箱廢氣溫度是空氣經過燃燒完全穿過料層后的溫度，由于在燒結過程中無法直接測量料層的溫度，因此風箱內的廢氣溫度作為能夠間接反映料層溫度的變量顯得尤為重要。燒結終點控制是燒結現場操作工調節設備參數的重要目標[13-15]，燒結終點是指燃燒深度剛好完全穿透料層時的位置，若燒結終點過于靠后，則會導致料層沒有燒透，影響燒結礦的質量，而若燒結終點過于提前，則會導致從燒結終點到燒結機機尾的這段距離內臺車篦條溫度持續過高，損壞設備，同時導致燒結礦在到達機尾時溫度偏低，不利于維持較好的燒結礦質量。由于以上狀態變量對燒結機生產狀態的重要指示意義，本文選取第14 號風箱（對應于燒結機中部）和第22號風箱（對應于燒結機尾部）的負壓和廢氣溫度以及燒結終點的位置和溫度共6 個狀態變量作為研究目標。

為了準確預測以上6 個關鍵狀態變量（state variables,SVs）的變化趨勢，需要分析這6 個變量變化背后的機理。通常來說狀態變量的變化與操作變量（operating variables,OVs）和狀態變量歷史值有關。本文中采用因果分析的手段來分析這種機理，具體來說，一方面通過自相關函數來分析狀態變量的自相關性窗口，另一方面通過收斂交叉映射方法來分析操作變量對于狀態變量的因果影響及其時間窗口。

1.2 狀態變量自相關性分析

燒結生產狀態具有一定的自相關性，即前一段時間的系統狀態變量會對之后的狀態變量產生明顯的因果作用，因此有必要對這種自相關性進行分析。自相關函數[16-17]（autocorrelation function, ACF）是用來衡量同一序列中不同時間間隔數據間的相關性隨時間間隔的變化的指標，其計算方法如式(1)所示：

其中，τ是滯后期，μx是序列均值。

選取六個系統狀態變量：14 號風箱負壓（OV1）、22 號風箱負壓（OV2）、14 號風箱廢氣溫度（OV3）、22 號風箱廢氣溫度（OV4）、燃燒終點位置（OV5）、燃燒終點溫度（OV6）進行自相關函數分析。結果如圖1 所示。可以發現隨著時滯期延長，所有六個變量的ACF 值均從1 開始逐步下降，取閾值水平為0.6，對應的平均自相關窗口長度為27 min。

1.3 操作變量對狀態變量的因果性分析

變量間因果關系能夠揭示系統中各因素之間的因果性，從而能夠幫助明確系統中某些變量發生變化背后的機理。變量因果性分析已有許多方法研究，如Bauer 等[18]提出了互相關分析（crosscorrelation analysis, CCA），Granger[19]提出了格蘭杰因果關系（granger causality, GC）檢驗，Schreiber[20]基于信息熵的概念提出了傳遞熵（transfer entropy,TE），但這些方法存在計算量大[21-22]或只適用于線性或弱非線性系統[23-27]等缺點，不能確保對強非線性系統中的耦合變量進行因果關系的正確識別，而Sugihara 等[28]提出的收斂交叉映射（convergent cross-mapping, CCM）可以克服此問題。McCracken等[29]證明了CCM 算法可用于雙向因果性推斷，M?nster 等[30]驗證了其在噪聲數據中的適應性，Luo等[31]進一步發展了CCM 算法以適應化工過程的需要。因此本文采用CCM 方法對于操作變量對狀態變量的因果性進行研究。

本文中所使用的CCM方法主要步驟如下：

（1）讀取數據并歸一化 從數據集中讀取兩個時間序列變量X =[x1,x2,…,xt]和Y =[y1,y2,…,yt]，并按照以下規則進行歸一化：

（2）重構流形 按照延時嵌入原則分別重構X和Y的吸引子流形：

其中i = 1,2,…,N -(E - 1)τ。E 是嵌入維度，τ是時間滯后。

（3）流形預測 在Mx上尋找Mx,i的E + 1 個由近 到 遠 近 鄰 點 Mx,i1,Mx,i2,…,Mx,iE+1，計 算Mx,i1,Mx,i2,…,Mx,iE+1離Mx,i的歐式距離，得到權重u：

圖1 6個狀態變量的自相關函數值Fig.1 Autocorrelation values of six state variables

在My上尋找對應的E + 1個點My,i1,My,i2,…,My,iE+1，計算使用Mx,i預測My,i的預測值：

（4）計算收斂交叉映射能力值

當考慮時滯時，CCM 可以通過數據平移的方法實現[32]，即在計算收斂交叉映射能力值時將真實值平移一個滯后期（lag），然后與預測值計算相關系數。

按照時滯算法取CMS 值的最大值作為有時滯情況下的操作變量對狀態變量的CMS 值，總結如表1所示。

分析表1 中的數據可以發現1 號主風機頻率和2 號主風機頻率對于14 號風箱和22 號風箱負壓的CMS 值均在0.6 左右，明顯高于系統中其他變量的CMS 值，說明風機頻率與風箱負壓具有強因果關系。這一點也能夠從機理上得到解釋，即主抽風機抽風的頻率決定了風箱負壓。

表1 操作變量對狀態變量的CMS值Table 1 CMS values of OVs towards SVs

除上述數據之外，表1 中的其他數據均位于0～0.4的范圍之內，說明在燒結系統中的數據因果性并不強。這是因為對于實際的工業系統，一方面各過程變量相互之間具有強耦合特征，另一方面由于變量之間的關系復雜，往往變量之間的關系不是直接相關，所以CMS 值會隨著傳播路徑的變長而減小。基于以上兩點，將燒結系統中變量因果性的閾值取為0.1，即認為在CMS＞0.1時變量之間存在明顯的因果關系。

在此閾值水平下，可以認為對于OV1～OV 3構成其原因的變量個數為9 個，對于OV6 構成其原因的變量個數為8 個，對于OV4～OV 5 構成其原因的變量個數為5 個。對于14 號風箱負壓和廢氣溫度、22 號風箱負壓這三個狀態變量，表1 中所選擇的9個操作變量均對其構成顯著因果關系，而對于22號風箱廢氣溫度，除點火溫度、臺車速度、圓輥速度、七輥速度CMS 值略低于0.1，其余操作變量均顯著。同理對于燒結終點位置，除了料層厚度、B 排點火強度、臺車速度、圓輥速度略低，其余也均構成顯著因果關系。對于燒結終點溫度，除七輥速度外，其余因果關系均顯著。從冶金實踐來看，料層厚度、圓輥速度和七輥速度與燒結床層的性質相關，點火強度和點火溫度與初始燃燒狀態相關，臺車速度、主風機頻率分別對應燃燒時間和風箱負壓，均與燃燒相關。綜上結合因果分析和機理分析，可認為本文中所選擇的9 個操作變量，對于6 個狀態變量均構成較好的因果解釋性，可以用于后續的預測。

設定樣本量為5000，取滯后期為-40～40，嵌入維度E=3，嵌入時滯τ=1，選取9 個操作變量對于狀態變量進行CCM 分析。圖2 是A 排點火強度（SV2）、點火溫度（SV4）、圓輥速度（SV6）和1號主風機頻率（SV8）對14 號風箱負壓和廢氣溫度的CMS值時滯曲線，可以發現CMS 值在20～40 min 的時滯期內均較高，此后逐漸降低。這一現象說明大部分的操作變量對狀態變量的因果關系存在明顯的滯后，滯后期為20～40 min，因此可將20～40 min 作為本文選取的操作變量對狀態變量影響的因果窗口。結合本數據集對應的冶煉生產實際情況來看，從進料到燒結完成的時間大致為40 min，與此處的因果性窗口較好對應。

圖2 CCM算法時滯分析結果Fig.2 CCM time lag analysis result

2 燒結生產狀態預測模型

2.1 模型結構

燒結狀態預測模型如圖3所示，共分為兩層：因果分析層和神經網絡層。因果分析層是通過數據因果性分析的方法獲得關鍵操作變量集、自相關窗口和因果影響窗口。具體來說，通過收斂交叉映射的方法從過程操作變量中選取與狀態變量有顯著因果關系的操作變量構成關鍵操作變量集，如圖3中OV1～OV9 所示，這些操作變量中蘊含了未來狀態變量的變化趨勢，可以作為后續神經網絡的輸入。除此之外，收斂交叉映射還可以得到操作變量集對狀態變量影響的時間窗口，如圖3中OV對應的黑色和深灰色色塊所示，這一窗口包含了變量間作用的時滯信息。自相關函數分析能夠對狀態變量的自相關性進行分析，從而獲得狀態變量的自相關窗口，如圖3 中SV 對應的黑色和深灰色色塊所示，這一窗口包含了變量自身作用的時滯。

神經網絡層中接收這些信息作為學習到的知識作為輸入，然后進一步學習輸出關鍵狀態變量的預測值。神經網絡層采取具有全連接層的人工神經網絡作為基本模型，采用開源的pytorch 庫構建神經網絡。優化器選擇torch.optim.Adam()，學習率是10-3。損失函數選擇torch.nn. L1Loss()，該損失函數對應的是平均絕對誤差（mean absolute error,MAE），數學表達如式(10)所示：

對于隱藏層，采用的激活函數是線性整流函數（ReLU 函數），該函數在x＜0 時的取值為0，在x＞0 時的取值為x，其數學表達如式(11)所示：

2.2 模型預測結果及精度分析

使用此模型對6 個狀態變量預測結果如圖4 所示，可以發現對于6個狀態變量，使用具有因果分析層和神經網絡層的生產狀態預測模型均能較準確地預測其變化趨勢。

將模型的預測精度總結見表2，可以發現對于狀態變量的預測平均相對誤差介于0.5%～3.1%之間，其中對于第14號風箱負壓的預測平均相對誤差是0.51%，而對于第14 號風箱廢氣溫度的預測平均相對誤差是3.07%，較精準地預測了狀態變量的變化趨勢，從而為現場操作工提前預知燒結機生產狀態變化，及時調整操作變量提供了參考。

3 結 論

圖3 燒結生產狀態預測模型框架Fig.3 Scheme of sintering production state prediction model

本文融合因果性機理和神經網絡黑箱模型，建立了燒結生產狀態預測模型。因果分析層使用自相關函數分析和收斂交叉映射的方法分別獲得了燒結狀態變量的自相關窗口、關鍵操作變量集及其對狀態變量的因果性影響窗口。神經網絡層采用誤差反向傳播神經網絡，接受因果分析層的輸入并對6 個關鍵狀態變量進行預測。經過工業數據測試，該模型能夠在3.1%的誤差范圍內對燒結生產狀態實現較精準的預測，從而能夠根據過程數據初步預判燒結機生產狀態的變化趨勢。為后續進一步研究燒結機生產狀態調整策略，實現對現場操作工的指導提供了基礎。

圖4 燒結生產狀態預測結果Fig.4 Prediction results of sintering production state

表2 燒結生產狀態預測模型預測誤差Table 2 Prediction errors of sintering production state prediction model

符 號 說 明

CMS——收斂交叉映射能力值

MAE——平均絕對誤差

OV——操作變量

SV——狀態變量