作为炼油生产过程的中间环节，加氢裂化流程受到上层调度命令和各项设备运行状况的限制，其运行状态经常发生改变。目前，在实际生产运行过程中，加氢裂化动态调整过程主要依靠现场操作人员凭借经验进行操作。为优化生产，可以利用加氢裂化历史数据，通过对历史已有的加氢裂化动态调整操作序列进行分析评价，选取历史最优的动态调整操作序列，建立最优动态调整智能决策系统指导生产。
历史最优加氢裂化动态调整操作序列的提取首先需要对动态调整区间进行识别。类似于流程工业的稳态判别，动态调整区间的识别是对工业过程从一种稳定运行模式向另一种稳定运行模式渐进式量变的动态过程的识别[1]，因此在加氢裂化调整过程识别的问题上，可以借鉴稳定操作状态(即稳态)的识别方法。
目前，稳态检测主要有组合检验法(CST法)[2]、滤波法[3]、滑动窗口法及F-检测法[4-5]等。陈文驰等提出一种融合自适应平滑技术的稳态检测方法，该方法可显著提高稳态检测处理的准确性[6]。文[7]提出对固定长度滑动窗口内数据通过最小二乘法得到二阶多项式系数，再根据反应系统稳态的指数来评价变量稳态。文[8]采用分段拟合、加权叠加的方法来检测稳态，将数据转换为多项式方程后，通过一阶导数和二阶导数的特性来判断系统稳态。文[9]提出了利用基于小波变换的多尺度表示来对测量数据进行有限连续逼近，通过分析该函数是否处于稳态来反映原信号的稳定情况。以上研究都是针对单一变量的稳定性进行分析的，对于多变量的情况，基本都是利用单变量的稳态检测结果来合成多变量的稳态结果。但是，加氢裂化过程是一个复杂体系，运行参数较多，各参数之间耦合严重，若是选择关键变量来进行过渡过程的识别，则不仅需要较为详细的机理知识而且关键变量数目也较多，一一进行分析会给计算带来很大负担。同时，实际过程数据中过渡状态前后数据特征不明显，容易被错误识别为稳态，从而对后续的分析造成较大影响[10]。
针对上述问题，本文提出一种基于主元导数特征聚类的加氢裂化动态调整区间识别方法，针对加氢裂化过程内部变量维度高、耦合强的特点，采用主元分析方法提取加氢裂化关键操作参数的主元，消除系统中冗余变量及各变量间随机噪声的影响；再基于滑动窗口的多项式拟合方法拟合主元数据，接着提取拟合数据的一阶导数作为聚类特征，采用导数特征来分析数据可以更加明显地反映数据的变化情况；设计基于密度峰值确定聚类初始点的K-means算法，进行聚类分析，从而识别出加氢裂化动态调整区间。采用中国某石化企业实际生产数据进行验证，结果表明本文所提方法可避免单个或几个变量误差的影响，有效识别动态调整区间，且不依赖先验知识，具有工业适用性。
1 问题描述 1.1 加氢裂化流程描述加氢裂化流程主要由加氢精制、加氢裂化、高低压分离、分馏系统这4个部分组成。加氢精制反应主要脱除原料油中的硫化物、氮化物和氧化物等非烃化合物。加氢裂化反应主要是使大分子裂解成小分子，使得产物中的氢含量提高，生成轻、中质产物，如图 1所示。加氢裂化主要处理的原料为：自常减压装置来的直馏轻蜡油、自罐区来的直馏轻蜡油、掺炼少量催化柴油。主要的产品为：轻石脑油、重石脑油、航空煤油、柴油、尾油。在反应过程中由于受到原料性质、催化剂活性等多种非线性因素的影响，加氢裂化流程进料流量、反应器的各个床层温度、急冷氢流量、分馏部分的汽提塔、分馏塔的塔顶温度、塔顶压力、塔顶回流量等需要调节的变量众多。
|
| 图 1 (网络版彩图)加氢裂化流程图 |
1.2 动态调整区间识别的问题描述
在实际的过程系统中，当系统受到外来干扰或者改变了设定值后, 原来的稳态就会遭到破坏, 系统中各组成部分的输入、输出量都相应发生变化，经过一段时间被控变量会变回原来的状态，这种从一个平衡状态到达另一个平衡状态的过程称之为过渡过程。
在理想情况下，对于只有一个变量x(t)的系统，若满足
| $ F\left( t \right) = \left\{ \begin{array}{l} {c_1}, \;\;\;\;\;\;\;\;\;\;t \in ({t_0}, {t_1});\\ {f_1}\left( t \right), \;\;\;\;{\rm{ }}\;t \in ({t_1}, {t_2});\\ {c_2}, \;\;\;\;\;\;\;\;\;\;t \in ({t_2}, {t_3});\\ {f_2}\left( t \right), \;\;\;\;{\rm{ }}\;t \in ({t_3}, {t_4});\\ {c_3}, \;\;\;\;\;\;\;\;\;\;t \in ({t_4}, {t_5});\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \end{array} \right. $ |
则系统在t∈(t1, t2)和t∈(t3, t4)分别存在两个不同的过渡状态，如图 2所示。其中：t0为系统初始时间，c1、c2和c3为常数。
|
| 图 2 理想情况下系统动态调整状态示意图 |
对于多变量系统X=[x1, x2, …, xn]，变量x1的过渡区间为T1，其中t∈(t11, t12)；变量x2的过渡区间为T2，其中t∈(t21, t22)；依次类推，变量xn的过渡区间为Tn，其中t∈(tn1, tn2)。记系统的过渡过程区间为T=T1∪T2…∪Tn。
在实际的工业过程中，对多变量复杂系统进行动态调整区间识别时，一般认为当系统的各个变量或者多数关键变量都处于过渡过程则系统处于过渡状态。
2 基于主元导数聚类特征的识别 2.1 总体思路针对加氢裂化流程长、变量多、耦合性强且实际过程数据中由于过渡状态前后数据特征不明显而容易被错误识别到稳定状态的问题，首先利用主元分析(PCA)方法将变量维度降下来；对降维之后的数据，利用滑动窗口，结合实际经验选取窗口长度，对这个时间窗口内的过程变量进行多项式拟合，再对拟合之后的各主元变量的曲线求取一阶导数。当工况在稳定运行状态和过渡过程状态之间切换时，各主元的一阶导数必然会发生较大的改变，采用导数特征来分析数据可以更加明显地反映数据的变化情况，克服实际过程数据过渡状态前后数据特征不明显、容易被错误识别到稳定状态的问题，为此将各主元变量的一阶导数作为聚类特征，利用基于密度峰值的K-means聚类方法对其进行聚类，从而区分加氢裂化动态调整过程。
2.2 基于PCA的主元分析对于加氢裂化这个复杂的多变量系统，利用主元分析的方法合并系统中的冗余变量从而降低系统的维度，同时这种方法也有一定的去噪效果。这种从原始空间到主元空间的映射不会改变数据的分布特性，通过变换可以提取对系统影响较大的变量。
加氢裂化流程生产过程，其历史数据可以用二维矩阵形式X∈Rn×z表示。其中: n和z分别表示过程的样本数和变量数。首先对数据进行标准化处理，消除过程变量测量值的量程差异对聚类结果的影响，标准化处理后的数据记为X1。
在信息损失最少的前提下，对原始数据进行压缩，提取数据特征，经过主元分析，历史运行数据被分解为
| $ \mathit{\boldsymbol{X}}_1^{\left( m \right)} = {\mathit{\boldsymbol{T}}^{\left( m \right)}}{\mathit{\boldsymbol{P}}^{\left( m \right){\rm{T}}}}. $ | (2) |
其中：T(m)为得分向量，P(m)为载荷向量，m(m≤z)为给定近似程度下所选取的主元个数。
2.3 主元数据导数特征的提取以第一主元x1为例来说明主元数据的特征提取过程。
选取长度为H的切割窗口对降维处理后的二维数据矩阵沿着采样方向进行切割，即将矩阵转置后沿横轴方向分割。变量H是根据经验选取的最佳拟合滤波窗口长度，将历史数据划分为一系列长度为H的窗口，在每个窗口内对数据进行多项式滤波。
设x1在各时刻的位置变量记录为
| $ {\mathit{p}_1}\left( t \right) = {\mathit{x}_1}\left( t \right). $ | (3) |
各变量的一阶导数S反映了各曲线在第t个采样时刻的变化快慢程度，用曲线各点的一阶导数能够描述曲线的变化趋势，因此对滤波后的主元变量求一阶导数，
| $ {\mathit{S}_1}\left( t \right) = {{\mathit{p'}}_1}\left( t \right). $ | (4) |
基于此，可以实现对加氢裂化调整过程更准确的识别。
2.4 基于主元数据导数特征的聚类K-means算法[11-12]是由MacQueen在1967年提出的一种聚类算法。由于传统的K-means算法的初始点是随机选择的，易造成聚类结果的随机性，Rodriguez等提出了基于快速搜索和密度峰值发现的聚类(clustering by fast search and find of density peaks，CFSFDP)，该聚类算法通过密度峰值找到聚类中心，再按照一定的顺序将数据点归类到距离其最近且密度比它自身密度高的点所属的类中[13-14]。利用CFSFDP确定K-means算法的初始聚类中心可以很好地降低K-means算法因随机选择初始聚类中心造成结果的随机性。
基于密度峰值的K-means算法的初始点确定分为3步：
步骤1 根据输入的主元导数样本矩阵和参数dc，计算每个数据点的局部密度ρ；
步骤2 计算每个数据点的高密度距离δ；
步骤3 使用决策图来确定聚类中心。
上述步骤中，局部密度ρ为一定距离内数据点的数量，即对于数据点i，其局部密度ρi[13-14]为
| $ {\rho _\mathit{i}} = \sum\limits_j {\chi \left( {{d_{ij}} - {d_{\rm{c}}}} \right)} . $ | (5) |
其中：当x＜0时，χ(x)=1；当x≥0时，χ(x)=0。dij是数据点i与数据点j之间的距离，而dc是人为的指定的一个距离。
高密度距离δ的定义为：对于数据点i，它与更高密度数据点的距离δi[13-14]为
| $ {\delta _\mathit{i}} = \mathop {\min }\limits_{j:{\rho _\mathit{j}} > {\rho _\mathit{i}}} {d_{ij}}. $ | (6) |
初始聚类中心的决策图可以在二维坐标系下建立，决策图中横坐标为局部密度ρ，纵坐标为高密度距离δ。不同的类会被低密度区域隔开，类中心会有较高的局部密度值ρ和较大的高密度距离δ，在决策图中，可以很明显地发现拥有高密度值和高密度距离的聚类中心。
基于密度峰值的K-means算法步骤如下：
步骤1 输入样本矩阵的主元数据的一阶导数矩阵X∈Rm×n(n表示样本数，m表示主元个数)；
步骤2 设定聚类数k，设定dc，根据式(5)、(6)来计算每个数据点的局部密度ρ以及数据点的高密度距离δ；
步骤3 使用决策图确定聚类中心；
步骤4 计算样本点到聚类中心的距离，将各点划分到与其距离最近的类中，形成初始分类；
步骤5 计算各类样本均值作为新的聚类中心；
步骤6 反复执行步骤4和5，直到聚类中心点都不再改变或者小于一定的误差范围。
3 加氢裂化工业过程的应用实例为验证所提方法的有效性，以国内某炼油厂加氢裂化反应系统为研究对象，所用的数据都来源于该加氢裂化流程现场，共提取了162个位号的数据, 通过结合专家经验与基于遗传算法的粗糙集属性约简后，得到主要表征加氢裂化流程动态调整过程的21个关键变量：加氢裂化流程入口流量，加氢精制反应器、加氢裂化反应器各床层温度，冷低分油流量，热低分油流量，脱硫化氢汽提塔和分馏塔的塔顶温度、塔顶压力，脱硫化氢汽提塔的塔底流量，分馏塔的中段回流量以及轻石脑油、重石脑油、航空煤油、柴油、尾油的抽出量。在数据采集过程中，加氢裂化流程正好进行了一次负荷调整，为此选取此次调整的数据作为分析的数据样本(自2016年2月12日至2016年2月13日，每5 min采集一个样本，共288组数据)。主要表征加氢裂化流程动态调整过程的21个关键变量的变化曲线如图 3所示。
|
| 图 3 加氢裂化动态调整过程关键变量变化曲线 |
图 3中横坐标表示采样序列，纵坐标表示各采样点对应的变量的实时值(标准化后)。图 3中，通过人工经验对加氢裂化动态调整区间进行识别得到本数据样本的动态调整区间为[90, 210]。
通过计算实例中各样本点的局部密度值ρ和高密度距离δ，决策图可以视为建立在二维坐标系下的数据表示，如图 4所示。不同的类总是被低密度区域隔开，因此选择类中密度最高点(图 4中已圈出)作为K-means算法的初始点。
|
| 图 4 CFSFDP法确定K-means初始聚类中心的决策图 |
图 5和6分别是利用CST稳态检测法和本文方法得到的识别结果。图中横坐标表示样本序列，纵坐标刻度值1表示稳定运行状态，2表示调整过程。从图 5和6中可以看出，本文方法比CST法对加氢裂化动态调整区间的识别具有更好的连续性。
|
| 图 5 CST法对调整过程的识别结果 |
|
| 图 6 本文方法对调整过程的识别结果 |
表 1是3种识别方法对本文案例的调整过程的识别结果。可以看出，本文提出的方法能够较为准确地识别出动态调整过程，尤其是在利用导数特征分析之后，对调整过程初期和末期的识别准确率更高。同时，利用本文方法在另外5组加氢裂化工业数据上进行验证，并与CST法进行对比，结果见表 2。表 2结果表明，本文方法可以有效识别加氢裂化调整过程的调整区间。
| 组号 | 本文方法准确率 | CST方法准确率 |
| 1 | 96.88 | 86.11 |
| 2 | 95.48 | 90.62 |
| 3 | 92.71 | 88.88 |
| 4 | 83.00 | 73.21 |
| 5 | 93.41 | 81.94 |
| 平均 | 92.30 | 84.15 |
4 结论
加氢裂化流程是一个高维复杂系统，而且各个变量对过程变化的表征不一。本文将PCA方法和多项式拟合结合起来对过程变量进行分析，在信息损失最少的前提下，对原始数据进行压缩；然后，对主元变化曲线求取一阶导数，采用导数特征来分析数据可以更加明显地反映数据的变化情况，可以避免数据中动态调整区间前后变化不明显的问题；采用K-means聚类方式对主元一阶导数进行聚类，可以直接提取出加氢裂化动态调整区间，同时基于密度峰值确定聚类初始点，可以避免初始点随机选择造成结果的随机性。最后，通过国内某石化企业加氢裂化流程的历史数据，对历史上数次调整过程的动态调整区间进行识别，并与经典稳态检测方法(CST法)的结果进行对比，验证了本文提出的方法的有效性和合理性。
| [1] |
张淑美, 王福利, 谭帅, 等.
多模态过程的全自动离线模态识别方法[J]. 自动化学报, 2016, 42(1): 60–80.
ZHANG S M, WANG F L, TAN S, et al. A fully automatic offline mode identification method for multi-mode processes[J]. Acta Automatica Sinica, 2016, 42(1): 60–80. (in Chinese) |
| [2] | NARASIMHAN S, MAH R S H, TAMHANE A C, et al. A composite statistical test for detecting changes of steady states[J]. AIChE Journal, 1986, 32(9): 1409–1418. |
| [3] |
李博, 陈丙珍, 胡惠琴.
稳态过程在线数据校正技术的工业实施[J]. 石油化工, 2000, 29(10): 768–771.
LI B, CHEN B Z, HU H Q. Industrial application of steady state online data reconciliation[J]. Petrochemical Technology, 2000, 29(10): 768–771. DOI:10.3321/j.issn:1000-8144.2000.10.011 (in Chinese) |
| [4] | LI G H. Inverse lag synchronization in chaotic systems[J]. Chaos Solitons and Fractals, 2009, 40(3): 1076–1080. DOI:10.1016/j.chaos.2007.08.062 |
| [5] | CAO S L, RHINEHART R R. Critical values for a steady-state identifier[J]. Journal of Process Control, 1997, 7(2): 149–152. DOI:10.1016/S0959-1524(96)00026-1 |
| [6] |
陈文驰, 刘飞.
一种改进的基于多项式滤波的稳态检测方法[J]. 控制工程, 2012, 19(2): 13–15, 20.
CHEN W C, LIU F. An improved steady state identification method based on polynomial filtering[J]. Control Engineering of China, 2012, 19(2): 13–15, 20. (in Chinese) |
| [7] | TAO L L, LI C C, KONG X D, et al. Steady-state identification with gross errors for industrial process units[C]//201210th World Congress on Intelligent Control and Automation (WCICA). New York, USA, 2012: 4151-4154. |
| [8] |
吕游, 刘吉臻, 赵文杰, 等.
基于分段曲线拟合的稳态检测方法[J]. 仪器仪表学报, 2012, 33(1): 194–200.
LÜ Y, LIU J Z, ZHAO W J, et al. Steady-state detecting method based on piecewise curve fitting[J]. Chinese Journal of Scientific Instrument, 2012, 33(1): 194–200. (in Chinese) |
| [9] | JIANG T W, CHEN B Z, HE X R, et al. Application of steady-state detection method based on wavelet transform[J]. Computers and Chemical Engineering, 2003, 27(4): 569–578. DOI:10.1016/S0098-1354(02)00235-1 |
| [10] |
季一丁. 多变量复杂系统的稳态检测和提取方法研究[D]. 杭州: 浙江大学, 2016. JI Y D. Research on steady state detection and extraction methods for multivariable complex system[D]. Hangzhou: Zhejiang University, 2016. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10335-1016185583.htm |
| [11] | MENG J L, SHANG H K, BIAN L. The application on intrusion detection based on K-means cluster algorithm[C]//International Forum on Information Technology and Applications. New York, USA, 2009, 1: 150-152. |
| [12] | CHOI Y, PARK C, KWEON I S. Accelerated K-means clustering using binary random projection[C]//Asian Conference on Computer Vision. Singapore, 2014: 257-272. |
| [13] | RODRIGUEZ A, LAIO A. Clustering by fast search and find of density peaks[J]. Science, 2014, 344(6191): 1492–1496. DOI:10.1126/science.1242072 |
| [14] | WANG S L, WANG D K, LI C Y, et al. Clustering by fast search and find of density peaks with data field[J]. Chinese Journal of Electronics, 2016, 25(3): 397–402. |