兩段法固定床甲醇制芳烴產物預測多元非線性回歸模型

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

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

引 言

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

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

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

1 實驗部分

1.1 實驗原料及反應設備

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

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

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

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

1.2 實驗設計

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

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

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

2 建模過程

2.1 數據準備及預處理

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

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

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

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

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

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

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

表2 MTA部分實驗數據Table 2 MTA partial experimental data

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

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

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

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

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

通過變換

2.3 模型評價指標

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

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

3 結果與討論

3.1 多元回歸模型及統計學檢驗

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

表3中列出了四個模型的統計學檢驗結果。可以看出四個模型的實驗組數和參數維度都分別為57、28，模型都滿足F>10F0.05，表明在考察的甲醇體積空速、一段反應溫度、二段反應溫度、反應壓力以及催化劑使用時長范圍內，本次建模結果是合適的，可以用于描述MTA 產物分布干基質量占比的變化過程。

表3 模型統計學檢驗Table 3 Model statistical tests

3.2 模型測試

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

表4 訓練集和測試集比較Table 4 Comparison of training and test sets

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

圖3 MTA產物分布奇偶圖Fig.3 Parity plot of MTA product distribution

圖4 總體預測樣本的奇偶圖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 軟件對芳烴回歸方程進行敏感性分析，研究芳烴制備條件中主要影響因素以及因素間的交互作用對芳烴質量百分數的影響，獲得各因素間交互影響的三維立體圖，如圖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 附近；而在空速增大的情況下，隨著壓力升高產物中芳烴占比遞增，在空速增加為0.5 h-1時，產物中最大芳烴占比對應的壓力峰值上移到0.8 MPa，此時對應的MTA 產物中芳烴占比值為62%。表明在高壓大空速交互作用影響下，芳烴的選擇性更高。這是由于在其他條件不變的情況下，大壓力本身有助于烯烴氣體向芳烴的轉變，但是小空速下物料停留時間更長，在大壓力下更易產生積碳而大量附著在催化劑的活性中心，使催化劑由于大量積碳增加而降低其活性。由圖5（b）可以看出，當一段反應溫度在390～510℃的范圍內時，相同壓力下，芳烴占比隨著溫度的升高而降低，在高壓情況下，降低的更明顯。這是因為高溫更容易加速積碳的產生，而導致催化劑上活性位點減少，減緩了芳烴生成的反應速率。由于空速不同，隨著裝置運行時間的增加，累計甲醇進料量會以不同的速度增加，因此E、F具有一定相關性。在分析交互影響時，由于工藝條件的物理意義限制，如果確定E、F中的任意一個因素，另一個因素也被確定在一定的范圍內。因此，本文選擇分析圖5（c）中壓力與裝置運行時間的交互影響。而在所限制的研究范圍內，芳烴占比隨反應裝置運行時間的增加無明顯變化，但壓力對芳烴占比的影響符合圖5（a）中的分析。

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

4 結 論

兩段法固定床反應器動力學具有復雜性，建立兩段法固定床產物分布的精確預測模型對固定床反應器的設計和優化具有重要的作用。本文建立的二次多項式回歸模型，主要結論如下。

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

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

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

符 號 說 明

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

MSE——均方誤差

p——壓力，MPa

T1——一段中心溫度，℃

T2——二段中心溫度，℃

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

w——干基質量占比，%

Yi——產物i的干基質量占比，%

βij——模型i的第j個參數

ρ2——決定性指標

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

下角標

calc——模擬計算值

exp——實驗數據