兩段法固定床甲醇制芳烴產(chǎn)物預(yù)測多元非線性回歸模型

許雄飛，劉鵬龍，張瑋，許鑫，張侃，王俊文

（1 太原理工大學(xué)化學(xué)化工學(xué)院，山西 太原 030024； 2 中國科學(xué)院山西煤炭化學(xué)研究所，山西 太原 030001）

引 言

芳烴是重要的有機化工原料。其中，輕質(zhì)芳烴（苯、甲苯、二甲苯，即BTX）市場需求量巨大，與工業(yè)中的“三烯”（乙烯、丙烯、丁二烯）有著同等重要的地位[1]，在能源、交通、紡織服裝和飲料包裝等領(lǐng)域都有廣泛應(yīng)用[2-5]。針對我國“富煤貧油少氣”的能源現(xiàn)狀，煤基甲醇制芳烴（MTA）技術(shù)不僅可以降低芳烴生產(chǎn)成本，而且可以降低我國的對外石油依存度[6-8]。芳烴中BTX 的收率直接影響到MTA 工藝的經(jīng)濟價值，是重要的評價指標(biāo)之一。因此，準(zhǔn)確預(yù)測和優(yōu)化包括BTX 在內(nèi)的產(chǎn)物分布是一個重要課題。

目前，對MTA 產(chǎn)物分布預(yù)測的研究主要采用熱力學(xué)計算與動力學(xué)建模的方法。文獻[9-11]分析了熱力學(xué)計算平衡時產(chǎn)物分布隨溫度、壓力等工藝條件的變化規(guī)律，研究了工藝條件對MTA 反應(yīng)產(chǎn)物的熱力學(xué)規(guī)律。但是熱力學(xué)平衡計算與反應(yīng)路徑無關(guān)，且對反應(yīng)體系做了一定的簡化處理[9-10]，而真實的反應(yīng)需要考慮催化劑及反應(yīng)設(shè)備等問題，故實驗結(jié)果與其計算結(jié)果相比有較大誤差。因此，熱力學(xué)模型只能從能量變化的角度對工藝條件優(yōu)化做指導(dǎo)，不能用于產(chǎn)物分布的精確預(yù)測。施麗麗等[12]和徐亞榮等[13]構(gòu)建了五集總穩(wěn)態(tài)動力學(xué)模型，將MTA 反應(yīng)體系劃分為含氧化合物、烯烴、芳烴、烷烴以及C1五個集總，分析了產(chǎn)物分布隨溫度和空速的變化規(guī)律。Li 等[14-15]則建立失活動力學(xué)模型，將催化劑活性考慮進模型中，并將它作為影響因素，根據(jù)不同的反應(yīng)體系劃分，研究MTA 的失活動力學(xué)行為。上述兩種動力學(xué)模型具有明確的物理意義，但該數(shù)學(xué)模型的假設(shè)前提是等溫反應(yīng)，模型的適用范圍會受到簡化和假設(shè)的限制；并且MTA 反應(yīng)機理復(fù)雜[16-23]，建立精確的動力學(xué)模型需要對不同溫度下的平衡常數(shù)k進行參數(shù)估計[24-28]，而動力學(xué)參數(shù)估計涉及對常微分方程或偏微分方程組的反復(fù)迭代求解，求解速度較為緩慢[29]。在化工研究中，上述基于實驗數(shù)據(jù)的動力學(xué)模型仍然是預(yù)測MTA 產(chǎn)物分布的主要方法，但該方法未考慮因素之間的交互作用，而數(shù)據(jù)驅(qū)動的多元回歸模型對非線性高維度的數(shù)據(jù)具有很強的擬合能力[30]，可以研究各因素之間的交互影響，能夠有效解決動力學(xué)模型在求解速度、精度和通用性方面的不足。

鑒于此，本文提出一種多元非線性回歸分析的數(shù)據(jù)驅(qū)動建模方法，用于MTA 產(chǎn)物預(yù)測。以兩段式固定床MTA 裝置為研究對象，采用實驗設(shè)計的方法開展實驗，以實驗過程中6 個工藝條件（壓力、甲醇體積空速、一段中心溫度、二段中心溫度、裝置運行時間、累計甲醇進料量）為模型輸入，以4 種產(chǎn)物（C1、烷烴、烯烴、芳烴）的干基質(zhì)量占比為模型輸出，研究各個因素及其交互作用對產(chǎn)物分布的影響。最后對所建模型進行性能測試。

1 實驗部分

1.1 實驗原料及反應(yīng)設(shè)備

實驗原料主要包括催化劑與甲醇。催化劑為中國科學(xué)院山西煤炭化學(xué)研究所制備，其中一段反應(yīng)器所用催化劑為HZSM-5 催化劑，二段反應(yīng)器所用催化劑為Zn 改性的ZSM-5 催化劑。甲醇為工業(yè)甲醇，是山西焦煤五麟公司的產(chǎn)品。

實驗裝置采用如圖1所示的兩段式固定床甲醇制芳烴小試裝置。甲醇通過神州微科公司的平流泵預(yù)熱氣化后分別進入一段反應(yīng)器和二段反應(yīng)器，催化劑裝填在反應(yīng)器的中間部分，兩個反應(yīng)器分別設(shè)有6個測溫點來采集反應(yīng)的中心溫度。經(jīng)過兩段反應(yīng)器的反應(yīng)，高溫產(chǎn)物經(jīng)冷卻循環(huán)系統(tǒng)進行冷凝，通過氣液分離罐后液相產(chǎn)物進入儲罐，由儲罐出口取樣分析。氣相產(chǎn)物通過濕式流量計計量后收集分析。

圖1 兩段式固定床甲醇制芳烴小試裝置Fig.1 Small test unit for methanol to aromatic in two-stage fixed bed

兩段法固定床甲醇制芳烴反應(yīng)裝置的兩段反應(yīng)器規(guī)格相同，內(nèi)徑32 mm，長1000 mm。實驗前，分別將200 ml 的HZSM-5 催化劑和200 ml 的Zn 改性的ZSM-5 催化劑置于一二段反應(yīng)器的中部，上下兩端用惰性瓷環(huán)填充。反應(yīng)所得產(chǎn)物采用離線分析方法，氣相組成及含量分析采用GC-920 型氣相色譜儀，色譜柱為：5A 分子篩[熱導(dǎo)（TCD）檢測器]，碳分子篩[氫火焰離子（FID）檢測器]及改性氧化鋁（FID 檢測器），三根色譜柱所測氣相組分含量通過甲烷進行關(guān)聯(lián)，歸一化處理后得到氣相產(chǎn)物全組成。油相產(chǎn)物采用GC-950 型氣相色譜儀進行分析，色譜柱為Rtx-1 型毛細管柱(柱長60 m×內(nèi)徑0.25 mm×膜厚0.5μm)，F(xiàn)ID 檢測器。水相產(chǎn)物分析采用GC-950 型氣相色譜儀、GDX-405 填充柱、熱導(dǎo)檢測器（TCD）檢測。

1.2 實驗設(shè)計

甲醇制芳烴反應(yīng)機理極其復(fù)雜，一段反應(yīng)器主要是甲醇制烯烴，二段反應(yīng)器里是未反應(yīng)完的甲醇和烯烴在催化劑作用下轉(zhuǎn)化為以芳烴為主的產(chǎn)物。溫度、空速和壓力作為影響MTA 反應(yīng)的主要因素，會影響到產(chǎn)品的分布。在兩段式固定床MTA 工藝中，這些工藝條件分別是：壓力（p）、甲醇體積空速（τ）、一段中心溫度（T1）、二段中心溫度（T2）。根據(jù)先驗知識，這些工藝條件的變動范圍如表1所示。

表1 兩段式固定床MTA工藝條件的變動范圍Table 1 Change scope of process conditions of MTA in two-stage fixed bed

傳統(tǒng)化工實驗用的是簡單對比法，即改變一個因素而其他因素固定不變，逐步找到每個因素的最優(yōu)值，從而得到最佳工藝條件。這種實驗方法優(yōu)勢是實驗次數(shù)少，但是試驗點不具有代表性，無法分清因素的主次，未考慮因素之間的交互作用。而全面試驗法是將所有水平一一搭配起來，從而找到最好的生產(chǎn)條件[31]。但是這種方法實驗次數(shù)太多，特別費時。為了合理安排試驗，減少試驗次數(shù)，但又能獲得范圍更廣的數(shù)據(jù)，本文采用Design-Expert 軟件進行自定義四因素五水平的實驗設(shè)計。

2 建模過程

2.1 數(shù)據(jù)準(zhǔn)備及預(yù)處理

MTA 反應(yīng)機理極其復(fù)雜[16-23]，主要有甲醇脫水生成二甲醚反應(yīng)、烯烴的生成反應(yīng)以及烯烴發(fā)生聚合氫轉(zhuǎn)移、環(huán)化脫氫等二次反應(yīng)生成芳烴和氫轉(zhuǎn)移生成烷烴等反應(yīng)類型。故本文擬建立一個基于數(shù)據(jù)驅(qū)動的模型，采用多元非線性回歸分析方法實現(xiàn)對MTA 主要產(chǎn)物分布的預(yù)測。模型的預(yù)測精度及其效率取決于數(shù)據(jù)的質(zhì)量，因此需要對實驗數(shù)據(jù)進行預(yù)處理。

首先，根據(jù)收集的氣液相產(chǎn)品分析得到MTA 的主要產(chǎn)物。進而參考施麗麗等[12-13]的集總劃分，將產(chǎn)物劃分為：含氧化合物MDOH（液相中的甲醇和二甲醚）、C1類產(chǎn)物（氣相中的H2、CH4、CO 和CO2等副產(chǎn)品）、烷烴（氣相烷烴和油相中的C+5以上烷烴）、烯烴（氣相中的烯烴）、芳烴（油相中的芳烴）。張寶珠[9]在MTA 反應(yīng)熱力學(xué)計算結(jié)果中得知甲醇生成二甲醚反應(yīng)為可逆反應(yīng)，能夠快速達到平衡，因此，本文將含氧化合物MDOH 直接記作甲醇。其次，對兩段式固定床MTA 反應(yīng)進行質(zhì)量守恒計算，從而對數(shù)據(jù)進行篩選。實驗中的測量誤差以及一些系統(tǒng)誤差假設(shè)，通過進出反應(yīng)體系的物質(zhì)質(zhì)量得到平衡常數(shù)之后，本文選擇誤差絕對值在10%以內(nèi)的71組數(shù)據(jù)進行MTA產(chǎn)物預(yù)測研究。

MTA 產(chǎn)物分布的預(yù)測模型有6 個輸入變量，分別是壓力（A）、甲醇體積空速（B）、一段中心溫度（C）、二段中心溫度（D）、裝置的運行時間（E）和累計甲醇的進料量（F）。相對應(yīng)的輸出為4個變量，分別是MTA 主要產(chǎn)物分布C1、烷烴、烯烴、芳烴的干基質(zhì)量占比，即Y1～Y4。下面是輸入輸出變量的說明。

對于模型輸入變量，首先選擇實驗設(shè)計中的四個工藝參數(shù)作為模型的部分輸入變量（A～D）。此外，整個實驗過程是在兩段式固定床小試反應(yīng)裝置上進行的，總運行時間超過1663 h，中間換了7 次催化劑。隨著反應(yīng)時間的變化，催化劑的活性會發(fā)生變化，這會極大地影響到不同產(chǎn)物的選擇性。Aguayo 等[27]在之前的研究中也提到，由于積炭導(dǎo)致催化劑失活的結(jié)果，除了表現(xiàn)在甲醇轉(zhuǎn)化率降低之外，還有產(chǎn)物選擇性的改變。如圖2所示，是在工藝條件為：τ=0.3 ml/min、p=0.4 MPa、T1=450℃、T2=480℃的工況下的實驗結(jié)果，反應(yīng)了MTA 產(chǎn)物分布隨裝置運行時間的變化過程。圖中可以明顯看出，C1和烷烴隨著裝置運行時間的延長而減少，烯烴的變化不明顯，芳烴的占比由于二甲苯和C+9重芳烴的占比增加而增加。說明隨著裝置運行時間的增加，甲醇芳構(gòu)化的程度有所加深，二甲苯和C+9重芳烴的選擇性升高，可能是由于上層催化劑部分失活，甲醇及其中間產(chǎn)物在催化劑中間部分繼續(xù)芳構(gòu)化反應(yīng)導(dǎo)致的。因此，本文在MTA 產(chǎn)物分布的預(yù)測模型研究中考慮催化劑的活性因素，將兩段式固定床MTA 裝置的運行時間（TOS）和累計甲醇進料量（M）作為輸入變量E和F。

圖2 兩段法固定床MTA產(chǎn)物分布隨裝置運行時間的變化Fig.2 Distribution of MTA products in two-stage fixed bed with device running time

輸出變量選取產(chǎn)物的主要成分的干基質(zhì)量占比，由圖2 可得，在裝置運行時間長達227 h 的情況下，產(chǎn)物中甲醇的干基質(zhì)量占比仍然小于1%。因此，催化劑的轉(zhuǎn)化率很高，產(chǎn)物中的甲醇可以忽略不計。所以本文將預(yù)測重點放在C1、烷烴、烯烴、芳烴四種主要產(chǎn)物上，即輸出變量Y1～Y4。

綜上所述，6 元多項式回歸模型的輸入、輸出變量及部分實驗數(shù)據(jù)如表2所示。

表2 MTA部分實驗數(shù)據(jù)Table 2 MTA partial experimental data

為了防止輸入數(shù)據(jù)的量綱差異影響建模效果，在進行建模前對輸入數(shù)據(jù)X進行了歸一化處理，如式(1)所示。

2.2 考慮催化劑活性的多元回歸模型

多元回歸模型（MVR）是多輸入單輸出(MISO)回歸建模的著名方法。它是研究兩個或者兩個以上自變量與一個因變量之間的關(guān)系，反映一個因變量受多個自變量影響而產(chǎn)生相應(yīng)變化的規(guī)律，是一種基于多個變量之間線性以及非線性數(shù)學(xué)模型的統(tǒng)計方法[32]。并且，MVR 屬于數(shù)據(jù)驅(qū)動模型的范疇，不需要詳細的過程機理知識，具有很強的泛化能力[30]。

在本研究中，采用的是多元回歸中的二次多項式非線性回歸模型。選擇80%的實驗數(shù)據(jù)作為訓(xùn)練集進行模型構(gòu)建，20%的實驗數(shù)據(jù)用來測試模型的準(zhǔn)確性。為了使測試數(shù)據(jù)具有代表性，需考慮反應(yīng)中間催化劑置換的過程，因此在每次催化劑使用過程中選擇兩個數(shù)據(jù)點組成測試集。建立以下四個六元二次回歸方程，響應(yīng)值分別為C1類、烷烴、烯烴、芳烴的干基質(zhì)量占比（Y1～Y4），自變量為T1(A)、T2(B)、τ(C)、p(D)、TOS(E)、M(F)，模型通式如式（2）所示。

式中，Yi表示第i類產(chǎn)物的干基質(zhì)量占比，1≤i≤4，分別對應(yīng)于MTA 主要產(chǎn)物分布C1、烷烴、烯烴、芳烴的干基質(zhì)量占比；模型中包含6個自變量以及它們的交叉項和平方項作為影響產(chǎn)物分布的因素，對于每種產(chǎn)物Yi的模型，需要對βi,0～βi,27共28 個參數(shù)進行參數(shù)估計；εi為隨機誤差。

通過變換

2.3 模型評價指標(biāo)

為了更加直觀地評價模型的建模性能，本文引入了統(tǒng)計學(xué)檢驗指標(biāo)F比檢驗、決定性指標(biāo)ρ2和均方誤差MSE 作為評價指標(biāo)。一般認為決定性指標(biāo)ρ2>0.9，F(xiàn)>10F0.05時，該模型是合適的，其置信區(qū)間為95%[33]。而MSE 作為衡量樣本離散程度的指標(biāo)，越接近于0說明精度越高。

其中，yi為實驗值；f(xi)為模型計算的預(yù)測值；n為樣本數(shù)量；m為參數(shù)維度。

3 結(jié)果與討論

3.1 多元回歸模型及統(tǒng)計學(xué)檢驗

為驗證所提方法的有效性，采用Design-Expert軟件對訓(xùn)練集實驗數(shù)據(jù)進行擬合，計算平臺為：Intel（R）i7-10700F，16G RAM，Windows 10。建模過程共耗時4 s，相比傳統(tǒng)動力學(xué)模型參數(shù)估計動輒上百秒的時間[29]，本方法在求解速度上有很大的優(yōu)勢。本次回歸得到四個二階非線性多元回歸模型，并用于預(yù)測MTA的產(chǎn)物分布（Y1～Y4）。

表3中列出了四個模型的統(tǒng)計學(xué)檢驗結(jié)果?？梢钥闯鏊膫€模型的實驗組數(shù)和參數(shù)維度都分別為57、28，模型都滿足F>10F0.05，表明在考察的甲醇體積空速、一段反應(yīng)溫度、二段反應(yīng)溫度、反應(yīng)壓力以及催化劑使用時長范圍內(nèi)，本次建模結(jié)果是合適的，可以用于描述MTA 產(chǎn)物分布干基質(zhì)量占比的變化過程。

表3 模型統(tǒng)計學(xué)檢驗Table 3 Model statistical tests

3.2 模型測試

為了驗證模型的可靠性，用建模過程中未使用的14組測試數(shù)據(jù)（占總數(shù)據(jù)的20%）進行模型分析。表4 中列出了4 個模型及其總體在訓(xùn)練集和測試集上的表現(xiàn)，包括各自的決定性指標(biāo)ρ2和均方誤差MSE。由表4可知，所建4個模型的訓(xùn)練集的ρ2都大于0.98。相比已有的MTA 動力學(xué)模型[12-13]精度更高。另外，在測試集中C1類回歸模型ρ2較小，其他模型的ρ2都大于0.94。說明4 個模型各自在訓(xùn)練集上的表現(xiàn)高于測試集，但是總體上在測試集上決定性指標(biāo)ρ2= 0.9576，MSE=0.0037，表明建立的4 個6元二次非線性回歸模型具有很高的預(yù)測能力。

表4 訓(xùn)練集和測試集比較Table 4 Comparison of training and test sets

為了更直觀地分析模型在數(shù)據(jù)集上的表現(xiàn)，將模型計算結(jié)果與實驗數(shù)據(jù)對比。圖3 為MTA 產(chǎn)物分布干基質(zhì)量占比（Y1～Y4）的奇偶圖，圖4 為總體預(yù)測樣本干基質(zhì)量占比的奇偶圖。正方形和三角形分別代表訓(xùn)練集和測試集數(shù)據(jù)。兩側(cè)的虛線代表20%以內(nèi)的相對誤差。由圖可以看出，大部分訓(xùn)練集和測試集數(shù)據(jù)均勻分布在y=x線的兩側(cè)，意味著模型預(yù)測值和實驗值取得了良好的一致性，可以用于預(yù)測兩段式固定床MTA產(chǎn)物分布。

圖3 MTA產(chǎn)物分布奇偶圖Fig.3 Parity plot of MTA product distribution

圖4 總體預(yù)測樣本的奇偶圖Fig.4 Parity plot of overall predicted sample

3.3 芳烴模型的敏感性分析

對芳烴回歸方程進行方差分析如表5所示。

從表5中可以得到，在制備條件的影響因素中，B、C以及交互作用AB、AC、AE、AF顯著，二次項中A2、D2也較顯著，其他項的顯著性均不明顯（P>0.05）。利用Design Expert 軟件對芳烴回歸方程進行敏感性分析，研究芳烴制備條件中主要影響因素以及因素間的交互作用對芳烴質(zhì)量百分數(shù)的影響，獲得各因素間交互影響的三維立體圖，如圖5所示。

圖5 各變量間的交互影響Fig.5 Effect of the interaction among the various variables

表5 芳烴模型方差分析Table 5 Analysis of variance for aromatic model

圖5（a）為壓力與甲醇空速的交互影響，在小空速的情況下芳烴選擇性隨著壓力的增加先升高后降低，峰值在0.2 MPa 附近；而在空速增大的情況下，隨著壓力升高產(chǎn)物中芳烴占比遞增，在空速增加為0.5 h-1時，產(chǎn)物中最大芳烴占比對應(yīng)的壓力峰值上移到0.8 MPa，此時對應(yīng)的MTA 產(chǎn)物中芳烴占比值為62%。表明在高壓大空速交互作用影響下，芳烴的選擇性更高。這是由于在其他條件不變的情況下，大壓力本身有助于烯烴氣體向芳烴的轉(zhuǎn)變，但是小空速下物料停留時間更長，在大壓力下更易產(chǎn)生積碳而大量附著在催化劑的活性中心，使催化劑由于大量積碳增加而降低其活性。由圖5（b）可以看出，當(dāng)一段反應(yīng)溫度在390～510℃的范圍內(nèi)時，相同壓力下，芳烴占比隨著溫度的升高而降低，在高壓情況下，降低的更明顯。這是因為高溫更容易加速積碳的產(chǎn)生，而導(dǎo)致催化劑上活性位點減少，減緩了芳烴生成的反應(yīng)速率。由于空速不同，隨著裝置運行時間的增加，累計甲醇進料量會以不同的速度增加，因此E、F具有一定相關(guān)性。在分析交互影響時，由于工藝條件的物理意義限制，如果確定E、F中的任意一個因素，另一個因素也被確定在一定的范圍內(nèi)。因此，本文選擇分析圖5（c）中壓力與裝置運行時間的交互影響。而在所限制的研究范圍內(nèi)，芳烴占比隨反應(yīng)裝置運行時間的增加無明顯變化，但壓力對芳烴占比的影響符合圖5（a）中的分析。

在二階非線性多元回歸模型的分析中，操作條件對芳烴選擇性的影響與張寶珠[9]的研究結(jié)論一致，但是卻較之前有所擴充，得到在甲醇空速增大的情況下，產(chǎn)物中最大芳烴占比對應(yīng)的壓力峰值也隨之增高。因此，在不同的操作條件下，二階非線性多元回歸模型對兩段法固定床MTA 小試裝置的產(chǎn)物預(yù)測具有良好的顯著性。

4 結(jié) 論

兩段法固定床反應(yīng)器動力學(xué)具有復(fù)雜性，建立兩段法固定床產(chǎn)物分布的精確預(yù)測模型對固定床反應(yīng)器的設(shè)計和優(yōu)化具有重要的作用。本文建立的二次多項式回歸模型，主要結(jié)論如下。

（1）模型可以用于兩段法固定床MTA 反應(yīng)裝置的產(chǎn)物分布預(yù)測。測試樣本的總體決定性指標(biāo)ρ2為0.9576，均方誤差MSE為0.0037，具有產(chǎn)物預(yù)測精確度高、計算量小以及泛化性強的特點，并且模型計算值與實驗結(jié)果具有良好的一致性，符合兩段法固定床反應(yīng)器設(shè)計要求和發(fā)展方向。

（2）模型對反應(yīng)器操作條件的變化具有良好的顯著性。在研究的工藝條件范圍內(nèi)，小空速的情況下芳烴選擇性隨著壓力的增加先升高后降低，峰值在0.2 MPa 附近；而在空速增大的情況下，產(chǎn)物中最大芳烴占比對應(yīng)的壓力峰值也隨之增高。在空速增加為0.5 h-1時，產(chǎn)物中最大芳烴占比對應(yīng)的壓力峰值上移到0.8 MPa，此時對應(yīng)的MTA 產(chǎn)物中芳烴占比值為62%。表明在高壓大空速交互作用影響下，芳烴的選擇性更高。

（3）該建模方法可用于非線性高維度的復(fù)雜工藝，為實際工廠的反應(yīng)系統(tǒng)放大研究提供了思路。

符 號 說 明

M——MTA裝置累計甲醇進料量，g

MSE——均方誤差

p——壓力，MPa

T1——一段中心溫度，℃

T2——二段中心溫度，℃

TOS——MTA裝置運行時間，min

w——干基質(zhì)量占比，%

Yi——產(chǎn)物i的干基質(zhì)量占比，%

βij——模型i的第j個參數(shù)

ρ2——決定性指標(biāo)

τ——甲醇體積空速，h-1

下角標(biāo)

calc——模擬計算值

exp——實驗數(shù)據(jù)