黃健，趙眾

（北京化工大學信息科學與技術(shù)學院，北京100029）

引 言

延遲焦化裝置加工的是劣質(zhì)重質(zhì)原料油,具有高溫操作的特點,裝置相互關聯(lián)嚴重[1]。在延遲焦化裝置中，焦化加熱爐屬于其中的關鍵設備，焦化加熱爐的燃料氣消耗量超過總耗能77%[2]，在衡量加熱爐操作先進性方面，熱效率是其中一項關鍵指標[3]，然而加熱爐熱效率難以在線測量，實時估計更加困難。影響其熱效率的因素有很多[4]，出口排煙的溫度、爐膛煙氣的氧含量、入爐的空氣量都會對熱效率產(chǎn)生影響。實際操作中熱效率主要依靠離線分析[5]得到，通過離線分析數(shù)據(jù)有正平衡法與反平衡法兩種方式計算[6]。正平衡法通過直接測量加熱爐輸入熱量和輸出熱量計算得到熱效率。而對于反平衡計算方法，則是通過測試和計算加熱爐各項熱損失。國內(nèi)外主要通過反平衡法計算。粒子濾波要追溯到20 世紀40 年代提出的Monte Carlo 方法，Hammersley 等［7］基于貝葉斯采樣估計提出順序重要性采樣濾波思想，其基本原理是尋找一組在狀態(tài)空間傳播的隨機粒子（樣本）描述系統(tǒng)的狀態(tài)，通過Monte Carlo 方法處理貝葉斯估計中的積分運算，從而得到系統(tǒng)狀態(tài)的最小均方差估計。當粒子數(shù)量趨于無窮時可以逼近服從任意概率分布的系統(tǒng)狀態(tài)。Gordon 等[8]通過重采樣技術(shù)的引入有效地解決了樣本權(quán)值退化問題，提出了粒子濾波算法。近年來，在現(xiàn)代信號處理、通信、人工智能、生物信息學、計算機視覺、自主移動機器人、目標跟蹤及統(tǒng)計學等領域不斷發(fā)展。粒子濾波的研究在中國最早見于2003 年，袁澤劍等[9]提出了一種性能優(yōu)于標準粒子濾波的高斯-厄米特粒子濾波器，莫以為等[10]把粒子濾波用于混合系統(tǒng)狀態(tài)監(jiān)測與診斷。隨后，粒子濾波得到了國內(nèi)學者廣泛的關注并被嘗試應用于目標跟蹤、組合導航和通信等領域，但需要較多的粒子數(shù)和計算時間仍然是亟待解決的問題。本文結(jié)合管式加熱爐物料，熱量平衡條件，構(gòu)造了加熱爐的狀態(tài)空間預測模型并提出了一種基于粒子濾波的聯(lián)合估計方法對模型預測值進行實時估計，并分析了方法的收斂性。在中國石化塔河煉化大型延遲焦化裝置的應用結(jié)果證實所提方法的可行性和有效性，為實現(xiàn)加熱爐熱效率實時優(yōu)化控制奠定了基礎。

1 工藝流程簡介及需求分析

中國石化塔河煉化延遲焦化裝置采用的是三爐六塔結(jié)構(gòu)，工藝主要包括三部分，第一部分是加熱爐焦炭塔、第二部分是分餾塔、第三部分是穩(wěn)定吸收部分。原料渣油經(jīng)一系列換熱器預熱至280℃后，進入分餾塔下段的換熱洗滌區(qū)，與來自焦炭塔頂?shù)?20℃高溫油氣接觸換熱，油氣中的重蠟油以上的組分被冷凝下來作為循環(huán)油與原料油一起流入分餾塔底，與原料油匯在一起成為350℃焦化油，塔底焦化油分四路進入加熱爐加熱至500℃。加熱后的高溫油氣，通過四通閥進入焦炭塔進行反應。高溫焦化油在焦炭塔內(nèi)發(fā)生裂解、縮合等一系列反應，生成反應油氣和焦炭。焦炭由下至上聚結(jié)在焦炭塔內(nèi)，反應油氣由焦炭塔頂逸出，進入分餾塔換熱塔板下部，進入分餾塔的焦化油氣與原料進行接觸換熱，循環(huán)油流入塔底，換熱后油氣上升進入分餾段，從下往上分餾出蠟油、柴油、粗汽油和富氣[11]。其中加熱爐作為加熱原料的龍頭部分，其在加熱介質(zhì)的過程中，加熱爐熱效率是加熱爐的一項重要的操作指標[12]。

2 加熱爐熱效率狀態(tài)空間模型

管式加熱爐是將爐管中通過的油料(或其他介質(zhì))加熱到所需的溫度,以滿足生產(chǎn)工藝的需要[13]。工業(yè)管式加熱爐是石化行業(yè)中必不可少的裝置。它通過燃料燃燒產(chǎn)生的火焰與高溫煙氣，對爐管內(nèi)高速流動的原料進行加熱，使其達到工藝規(guī)定的溫度。被加熱介質(zhì)在對流室內(nèi)與煙氣進行對流換熱，在輻射室內(nèi)通過火焰或者高溫煙氣進行輻射傳熱。管式加熱爐的結(jié)構(gòu)如圖1所示。

圖1 工業(yè)管式加熱爐Fig.1 Industrial tubular heating furnace

根據(jù)別洛康方法[14]及微元法，有：

穩(wěn)態(tài)時，煙氣對輻射室爐墻的熱射熱速率與輻射室爐墻對輻射管反射的熱速率相同，有：

根據(jù)以上傳熱方程及熱平衡條件，燃料氣的流量調(diào)節(jié)回路，對加熱爐煙氣理論氧含量[15]計算如下：

而對于加熱爐熱效率，目前有兩種方法計算得到：正平衡法和反平衡法?？筛鶕?jù)實際爐子熱負荷和實際燃料消耗量計算，即正平衡方法，也可根據(jù)各種損失由反平衡法計算。本文根據(jù)反平衡法計算[16]：

對于本文研究的中石化某煉油廠的渣油，Q=39330 kJ/kg，L=13.225，H=202.6 kJ/kg。圓筒爐及四壁排管的加熱爐q3取3，有廢熱回收系統(tǒng)在公式外額外增加0.5，根據(jù)以上條件將式(6)簡化：

其中對于化學不完全燃燒熱損失q2，在一氧化碳所占煙氣百分比為1%以下時根據(jù)經(jīng)驗公式（8）計算得到：

式中，CO為煙氣中的一氧化碳濃度。根據(jù)化驗室得到的數(shù)據(jù)得知常年煙氣中的一氧化碳濃度在0.01%左右，可以使用式(8)。

熱效率還與過剩空氣系數(shù)α和排煙熱損失q1有關。在工業(yè)生產(chǎn)中，燃料在化學平衡中所需空氣量即理論空氣量下是不可能完全燃燒的,因此需要過提供一定量的空氣即過?？諝?以保證燃料的完全燃燒[17]，這就是過?？諝庀禂?shù)。各種燃煤的工業(yè)爐窯和鍋爐排出的煙溫較高，燃料熱的相當大部分未加利用，此即排煙熱損失[18]。熱效率與過?？諝庀禂?shù)和排煙熱損失的關系如圖2所示。相關參數(shù)具體推導見附錄1[19-21]。

結(jié)合式(6)～式(10)，根據(jù)中石化某煉油廠一個月的化驗數(shù)據(jù)回歸得到熱效率擬合輸出為：

圖2 熱效率與過?？諝庀禂?shù)和排煙熱損失的關系Fig.2 Thermal efficiency versus excess air and flue gas heat loss

其中，O2為煙氣中氧含量，T為出口煙氣溫度。結(jié)合式(5)～式(11)選取出口介質(zhì)溫度Tof，爐膛煙氣溫度Tfg，煙氣氧含量Ofg，燃料氣進料Fg為狀態(tài)變量，爐膛煙氣出口溫度Tofg，燃料氣進料設定值Fsg，單位時間內(nèi)加熱爐空氣流入的體積Fa與進料介質(zhì)溫度Tif為輸入變量，W為狀態(tài)噪聲，V為觀測噪聲。得到如下狀態(tài)空間方程：

其中

3 基于粒子濾波聯(lián)合估計的加熱爐熱效率實時估計

3.1 Kalman濾波

實際工業(yè)過程由于建模偏差、噪聲及擾動的影響，會造成質(zhì)量指標的模型預測偏差，因此需要利用離線分析數(shù)據(jù)對模型預測進行實時反饋修正[22]。在線性離散系統(tǒng)中經(jīng)典Kalman 濾波（KF）是一種最優(yōu)估計方法[23]?？紤]如下的線性離散動態(tài)系統(tǒng)：

其中xk是k時刻的系統(tǒng)狀態(tài)，uk是k時刻對系統(tǒng)的控制量，wk和vk分別表示過程和測量的噪聲，假設成高斯白噪聲，協(xié)方差分別為Q和R。假設：

（1）過程噪聲wk和測量噪聲vk均是高斯白噪聲序列，且二者互不相關或者二者δ相關，即

式中，μw是wk的均值矩陣，μv是vk的均值矩陣，Qk為非負定方差矩陣，是系統(tǒng)過程噪聲wk的方差矩陣，Rk為正定方差矩陣，是系統(tǒng)測量噪聲vk的方差矩陣，Sk是wk和vk的互協(xié)方差矩陣，δkj是克羅尼克δ符號。

（2）系統(tǒng)初始狀態(tài)X0是某種已知分布或正態(tài)分布的隨機向量，其均值、協(xié)方差分別是:

（3）系統(tǒng)的過程噪聲wk和測量噪聲vk都與初始狀態(tài)X0互不相關，即

對于式(14)所示系統(tǒng)，滿足以上條件，其時間更新方程[24]如下：

3.2 粒子濾波

粒子濾波（PF）[25]是一種基于貝葉斯濾波思想，應用重要性采樣的概念，來獲得狀態(tài)向量后驗概率密度的序貫蒙特卡洛方法。考慮式（21）描述系統(tǒng)：

wk-1和vk為獨立同分布噪聲（非白噪聲），根據(jù)貝葉斯定理，假設xk服從一階馬爾科夫過程，且狀態(tài)序列xk和量測序列yk相互獨立，則狀態(tài)預測方程和更新方程分別如式（22）、式（23）所示：

其中p(yk|y1:k-1)為歸一化常數(shù)。在對后驗分布進行蒙特卡羅數(shù)值模擬計算時，直接利用未知的后驗概率密度是非常困難的。因此，根據(jù)貝葉斯重要性采樣定理[26-27]，引入一個容易采樣的參考分布q(x0:k|y1:k)，則后驗概率密度函數(shù)可近似為：

歸一化重要性權(quán)值，就可由式（24）計算得到狀態(tài)的概率分布，然而，隨著序貫重要性采樣算法的進行，一部分粒子會因為權(quán)值逐漸變小而失去作用，這對粒子濾波的進行非常不利。為了克服這種狀況，加入重采樣策略[28-29]，其基本思想是剔除權(quán)值小的粒子，同時對大權(quán)值粒子根據(jù)權(quán)值大小進行復制。

3.3 聯(lián)合估計

粒子濾波可以很好地克服有色噪聲對模型帶來的影響，但在一些過程中PF為了滿足估計精度要求，就需要更多的粒子數(shù)和更長的計算時間，這樣就大大削弱了實時估計的優(yōu)勢[30]。針對這種情況，本文參照狀態(tài)與未知參數(shù)聯(lián)合估計方法[31]，結(jié)合KF 和PF，提出了基于粒子濾波的狀態(tài)和校正系數(shù)聯(lián)合估計方法（JPF）來修正預測模型。考慮如下系統(tǒng)：

基于預測模型的狀態(tài)估計如下所示：

雖然對于式（27）可以利用KF 進行估計，但是這樣不能很好地處理非高斯噪聲和模型的大誤差。并且KF 穩(wěn)態(tài)濾波也不是常數(shù)，而應該是一個變值。故可以將K描述成：

其中K0由KF 得到，ξ3為噪聲。于是，可將系統(tǒng)擴維成：

從這個擴維系統(tǒng)中可以看出，濾波增益的動態(tài)方程中由于存在噪聲，整個增益輸出不夠平滑。為此，對其采用一階慣性濾波，即：

其中λ(0 <λ< 1)為濾波常數(shù)，整理得到最終的擴維系統(tǒng)如式(31)所示：

證明見附錄2[32-33]。由定理可見，當粒子數(shù)和時間趨于無窮時，該方法收斂于最優(yōu)濾波。

4 工業(yè)應用

本文建立了中石化某煉油廠如圖3所示的延遲焦化加熱爐裝置部分的狀態(tài)空間模型。熱效率輸出模型參數(shù)使用最小二乘法擬合得到，對所得到的模型進行模型預測驗證結(jié)果如圖4 所示，紅色點劃線部分為廠里化驗分析數(shù)據(jù)，藍色實線部分為狀態(tài)空間模型預測值。數(shù)據(jù)分析對比見表1，可以看到通過預測模型計算得到的熱效率數(shù)值與化驗分析得到的數(shù)值的誤差不大于0.3，誤差相比較分析計算得出的熱效率誤差不大于0.3%，完全可以作為加熱爐熱效率的實時輸出。將化驗分析值分別與聯(lián)合粒子濾波、粒子濾波和經(jīng)典Kalman濾波的估計值進行比較，效果如圖5所示，具體結(jié)果見表2。

圖3 延遲焦化流程圖Fig.3 Process flow diagram

圖4 模型計算對比Fig.4 Model prediction comparison with analysis data

表1 模型計算對比Table 1 Model prediction results

圖5 濾波效果圖Fig.5 Filter effect diagram

表2 各種濾波比較Table 2 RMSE comparisons with different filters

通過圖5 與表2 可以看出運用聯(lián)合粒子濾波比單獨使用Kalman濾波或單獨使用粒子濾波，標準誤差分別下降61.53%和70.18%。粒子濾波聯(lián)合估計的精度要比PF 和KF 高，這說明聯(lián)合粒子濾波在整體性能上是優(yōu)于粒子濾波和經(jīng)典Kalman 濾波的。但無論粒子濾波還是聯(lián)合粒子濾波都存在合適粒子數(shù)問題，具體分析結(jié)果見表3。

表3 粒子濾波與聯(lián)合粒子濾波對比Table 3 Comparison of particle filtering with joint particle filtering（PF/JPF）

通過表3 可以看出隨著粒子數(shù)的增加，無論是PF 還是JPF，均方根誤差（RMSE）都在減小，相對應的運行時間均在延長。但二者處于相同的均方根誤差（RMSE）時，使用JPF 比PF 的運行時間縮短20%，很好地克服了粒子濾波的需要更多的粒子數(shù)和更長的計算時間的問題。將模型參數(shù)由上層工程師站組態(tài)進DCS 中，建立軟測量界面，加熱爐熱效率的實時輸出如圖6 所示，為實現(xiàn)裝置的先進控制奠定了基礎。

圖6 熱效率在線輸出曲線Fig.6 Thermal efficiency online output

5 結(jié) 論

對于延遲焦化裝置，管式加熱爐消耗了能耗的大部分，但熱效率并不能實時測量，操作人員不能根據(jù)熱效率的實時變化進行相應的操作，以達到節(jié)能降耗的目的。本文結(jié)合管式加熱爐物料，熱量平衡條件，構(gòu)造了加熱爐的狀態(tài)空間方程并提出了一種基于粒子濾波的聯(lián)合方法對模型預測值進行實時估計，并對方法收斂性進行分析，工業(yè)應用結(jié)果證實所提方法的可行性和有效性。

符 號 說 明

A——爐管有效傳熱面積，m2

Ai,Ao,Am——分別為管內(nèi)、管外、平均表面積，m2

AR——輻射外表面積，m2

Cf,Cfg——分別為被加熱介質(zhì)焓值與煙氣焓值，kJ/kg

C0——煙氣氧含量的容量系數(shù)

Fa,Fif,Fofg——分別為單位時刻內(nèi)加熱爐空氣流入量、被加熱介質(zhì)流量以及煙氣出口流量，kg/h

Fg,Fsg——分別為燃料氣流量和流量設定值，kg/h

Hs——絕對黑體表面熱度，J/(K4·s)

hRC,kc——分別為煙氣傳熱系數(shù)與總的傳熱系數(shù)，W/(m2·K)

L——單位質(zhì)量燃料完全燃燒所需要的理論空氣量，是其比值

Ofg——爐膛內(nèi)煙氣氧含量，%

Q,H——分別為燃料低熱值與單位質(zhì)量空氣在預熱溫度下的熱焓，kJ/kg

QR——輻射室的傳熱速率，J/s

q1,q2,q3——分別為排煙熱損失、化學不完全燃燒熱損失與爐墻散熱損失，%

Ri,Ro——分別為加熱爐管內(nèi)與管外結(jié)構(gòu)熱阻，K/W

r——混合程度系數(shù)

Tg,Tt——分別為離開輻射室煙氣溫度及輻射外表面的平均溫度，K

Tif,Tof,Tfg,Tofg——分別為被加熱介質(zhì)的入爐溫度、出爐溫度、爐膛煙氣溫度以及爐出口煙氣溫度，K

U——平均表面熱強度，kJ/(m2·h)

Vf,Vfg——分別為進料爐管體積與煙道體積，m3

Vmf——在一定熱負荷下理論上所需要氧氣量，m3

α——過?？諝庀禂?shù)

εg,εr,εw——分別為煙氣、輻射管表面和輻射室爐墻的黑度

ρf,ρfg——分別為被加熱介質(zhì)與煙氣密度，kg/m3

τg——燃料氣回路的時間常數(shù)