许雄飞，刘鹏龙，张玮，许鑫，张侃，王俊文

（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——实验数据