基于趋势分析的间歇过程异常工况超早期报警研究

胡瑾秋， 郭 放， 张来斌

(中国石油大学 油气资源与工程国家重点实验室 机械与储运工程学院， 北京 102249)

间歇过程具有灵活多变的特点，适用于小批量、高附加值、多品种的工业要求，在精细化工生产中得到广泛应用[1]。间歇过程反应复杂，具有多时段特性，过程参数相关关系随操作进程呈周期性的阶跃变化，仅通过分布式控制系统(Distributed control system, DCS)监测参数是否超出阈值范围并不能有效地判断系统运行情况[2]。参数的变化趋势是判断阶跃变化中工况是否异常的重要指标，阶跃变化过程中，若参数变化速率控制不当，极易出现飞温、飞压等危险情况，造成严重后果[3]。

目前国内外对于化工间歇过程的研究多集中于多时段工况监测[4]，1994年首次提出针对间歇过程的不同子时段，采用多向主元分析(Multiple principal component analysis, MPCA)逐个建立监测模型。郭小萍等[5]提出基于滑动窗口的多时段主元分析(Principal component analysis, PCA)建模方法，解决了短期内难以取得理想数据的问题。常玉清等[6]提出一种改进的多时段MPCA间歇过程监测方法，根据不同子时段主成分特性的不同对间歇过程进行划分，再利用MPCA方法实现各个子时段的在线监测，但该方法需要对各个时段的划分有一定的先验知识。王静等[7]提出一种主元分析-多方向高斯混合模型(Principal component analysis-multiple Gaussian mixture model, PCA-MGMM)，先采用PCA对数据降维，在低维空间中用GMM对得分向量训练，将局部概率指标融合为全局概率监控指标，该方法不需要模型的先验知识。这些方法需要将过程数据划分成若干个阶段逐一建模，并假设各个阶段服从高斯分布[8]，不仅计算复杂，而且在实际过程监测中实时性差。为了更好地实时监测工况，袁彬等[9]在3σ准则的基础上提出一种浮动报警阈值的方法，根据系统的危险状况，建立一级和二级浮动报警值，与固定报警阈值相比，能够及时有效地发出警报。朱海等[10]将3σ准则应用于自适应小波阈值的选取，能够有效地区分海洋背景噪声和磁目标信号，剔除异常值，使小波阈值降噪效果更好。但是3σ准则成立的条件仍是假设数据服从高斯分布，不适于阶跃变化过程。最小二乘法[11-13]作为常用的数据拟合方法，可以简便地求得未知数据，并使求得数据与实际数据之间误差的平方和最小。笔者提出一种基于最小二乘法的拟合-微分-再微分的趋势分析方法，通过拟合-微分之后过程参数的一阶变化率进行间歇过程多时段工况识别，针对阶跃变化过程，利用二次微分对参数一阶变化率的变化快慢分析，结合滑动窗算法，实现实时的间歇过程异常工况超早期报警。

1 基于最小二乘的趋势分析原理

最小二乘法(Least-squares，LS)是根据“使偏差平方和最小”的原则来选取拟合曲线y=f(x)。利用最小二乘法解决问题需要两个步骤：根据已知数据确定f(x)的具体形式；按照最小二乘的原则求得最小二乘解。经常采用的最小二乘多项式方程形式为：

f(x)=anxn+an-1xn-1+…+a1x+a0

(1)

式(1)为n次多项式拟合形式。最小二乘法的目的就是要确定ai，并且是各个点的偏差δi平方和最小。将m对数据(xi,yi)代入式(1)中，得到如下方程：

(2)

将式(2)转换成矩阵的形式，求得矩阵的唯一的一组最优近似解，使得偏差δi的平方和最小，从而求得最小二乘拟合多项式。

2 基于趋势分析的超早期报警方法

2.1 工况识别

假设滑动窗口长度L是曲线拟合需要的数据长度，步长D是一次拟合数据更新时窗口添加和丢弃的数据长度。滑动窗口的选取对于最小二乘拟合的准确度非常关键。滑动窗长度太小，拟合曲线与实际趋势不匹配；长度太大，计算量大，拟合速率慢。由于间歇过程中参数变化频繁，故选取较小窗口长度。取当前时刻前L个数据作为拟合所需数据，记为X0∈RL0×m，m为过程变量个数，L0为初始拟合数据长度，每隔步长D对拟合数据进行一次更新Xk∈RLk×m，k为数据的更新次数。Lk为第k次更新训练数据的长度。其中：

(3)

选取历史数据中正常样本数据根据公式(2)进行最小二乘曲线拟合，求出长度为L的数据的最佳拟合函数，并进行一阶、二阶微分计算得到参数的一阶和二阶微分方程。一阶方程式反映参数变化率υ的情况，用来识别工况处于的阶段；二阶方程式反映参数的二阶变化率γ的情况，用来判断参数在某一工况下是否出现过快或过慢变化。

间歇过程工况大致可分为3个阶段：上升阶段、反应阶段和下降阶段。假设上升阶段中参数X的值为α，一阶变化率为υ。根据参数和一阶变化率的数值区间，判别工况的运行阶段。3个阶段的参数X和其变化率υ的取值范围如表1所示。

2.2 参数超早期报警

所谓超早期报警，即在异常发生的早期就发出警报，提前时间应在5 min以上。确定工况处于上升阶段或下降阶段后，采取拟合、微分、再微分的方法，计算参数X和参数的二阶变化率γ的取值区间。每个阶段分为2种趋势：稳定趋势和不稳定趋势。不同趋势下参数X和参数的二阶变化率γ的取值范围如表2所示。

表1 参数X和一阶变化率(υ)的取值区间Table 1 The range of values for parameters X and first order derivative (υ)

表2 参数X和二阶变化率(γ)的取值区间Table 2 The range of values for parameters X and its second order derivative (γ)

数据拟合之后通过与表2中的取值范围比较，判断是否处于相应的区间内，若在区间范围内，则参数变化正常，否则，参数变化异常。整个工况结束后，将3个阶段的数据拟合，并进行一阶、二阶微分求导，得出表1和表2中的各个参数的取值区间。随后将最大值和最小值保存到样本库中，并与相应的样本库中的数据比较，选取样本库中的最大值和最小值作为最新的控制限。根据表1判断间歇过程工况阶段，识别为上升阶段或者下降阶段，进行变化率监测。利用最小二乘法对当前T时刻长度为L的实时数据拟合，微分计算后与表2中的取值范围比较，判断是否处于相应的区间内，若在区间范围内，则参数变化正常；否则，参数变化异常，发出报警信号。

3 实施步骤

针对化工间歇过程工况识别以及阶跃变化报警监测的问题，提出基于最小二乘的趋势分析方法，步骤如下：

步骤1：选取间歇过程中存在阶跃变化的监测变量的10组不同时长的正常历史数据作为训练样本。间歇过程中参数变化范围大，剔除数据中离群点会忽略部分正常生产过程信息，降低报警准确度，因此，不对历史数据进行预处理。依据公式(1)对历史数据利用最小二乘法进行多项式拟合，求得拟合后的方程f(x)，并画出参数拟合曲线；

步骤2：根据工艺操作要求和步骤1中得到的拟合曲线趋势，对间歇过程的各个工段进行划分，分为上升阶段、反应阶段和下降阶段，并确定各个阶段中参数X的取值区间；

步骤3：对历史数据进行一阶微分求导，画出参数一阶变化率的曲线。根据步骤2中划分的各个工段的参数数值，确定各个工段的参数一阶变化率υ的取值区间；

步骤4：选取间歇过程中工艺操作较危险的阶跃变化阶段，对该段的数据进行再微分处理，画出参数的二阶变化率γ的曲线。根据二阶变化率的大小，将该阶跃变化阶段划分成稳定阶段和不稳定阶段。稳定阶段为参数在这个阶段保持较稳定的速率上升或下降，不稳定阶段为参数在这个阶段呈现出不太稳定的上升或下降。结合该工段的参数取值，确定稳定阶段和不稳定阶段的参数二阶变化率γ取值区间；

步骤5：由于间歇过程参数波动较大，窗体大小和步长不宜过大。选取窗体大小为400～500 s，步长为5～10 s。对在线数据实时进行拟合-微分-再微分，将在线数据的一阶变化率和步骤3求得的一阶变化率取值区间比较，识别出间歇过程的具体工段。将在线数据的二阶变化率与步骤4得出的二阶变化率取值区间比较，当参数超出取值区间时，则发出警报。

4 案例分析

4.1 工艺介绍

采用常见的间歇生产装置——聚丙烯装置的聚丙烯聚合过程作为研究对象，使用现场采集数据进行间歇过程早期异常情况识别。聚合釜是聚丙烯装置中的重要部分，它使丙烯单体在催化剂、氢气、活化剂和催化助剂的作用下自身发生加成聚合生成聚丙烯[14-15]。聚合釜反应过程中的主要工艺监测参数有聚合釜上温、聚合釜中温、聚合釜下温和聚合釜压力。聚合反应属于剧烈反应，冷却系统需保持良好水平，若温度或压力上升过快，可能引发爆炸。以聚合釜上温为例进行说明。

4.2 趋势分析报警监测

(1)分别选取历史数据中连续10组样本求取平均值，利用最小二乘法对数据进行多项式拟合，并对样本进行工段划分，图1为3个阶段中聚合釜上温的拟合曲线图，图中的两个采样点之间隔有20个点。

图1 聚合釜上温3个阶段最小二乘拟合曲线Fig.1 The curve of the polymerization reactor temperature in three stages by least squares

(2)对历史数据曲线进行一阶微分计算，根据已经划分的3个阶段的参数取值，分别确定3个阶段的一阶变化率υ，如图2所示。

图2 聚合釜上温和一阶变化率曲线Fig.2 The curve of the polymerization reactor temperature vs the first order derivative

根据图1中的阶段划分和图2中一阶变化率的变化趋势，确定一阶变化率的取值区间，如表3所示。

表3 聚合釜上温和一阶变化率取值区间Table 3 The range of values for polymerization reactor upper temperature and its first order derivative

(3)由于升温过程较危险，因此选取聚合釜上温的升温阶段重点进行监测。利用再微分方法计算参数二阶变化率，并根据二阶变化率的值将升温阶段划分成稳定升温阶段和不稳定升温阶段，如图3所示。通过二阶变化率曲线趋势，确定聚合釜上温和二阶变化率γ的取值区间，如表4所示。

图3 聚合釜上温二阶变化率曲线Fig.3 The curve of the polymerization reactor upper temperature vs the second order derivative

TrendXγStable22-29,56-60-0.07-0.07Unstable29-56-0.15-0.15

(4)由于升温过程用时较短，滑动窗的长度也应较短。设置滑动窗长度为400 s，步长为5 s。选取一组聚合釜上温升温过快的异常数据，共511个采样点(约42 min)，一般升温过程持续50 min以上。前400 s的聚合釜上温温度及一阶变化率曲线如图4和图5所示。由图4、图5可知，该段故障数据处于聚合过程中升温阶段。

图4 聚合釜上温数据曲线Fig.4 The curve of the polymerization reactor upper temperature

图5 聚合釜上温一阶变化率曲线Fig.5 The curve of the polymerization reactor upper temperature vs the first order derivative

监测过程中，一阶变化率曲线在第805 s时超出控制限，但瞬间恢复正常，此时的聚合釜上温为32℃。之后在第1080 s，聚合釜上温达到39℃时，对照表4，该时刻处于不稳定升温阶段，且开始连续超出报警限，如图6所示。利用再微分得到该段数据的二阶变化率曲线，如图7所示。图7中，在第1080 s二阶变化率同样超出报警线，说明在该采样点处聚合釜上温存在升温速率过快的危险。

图6 聚合釜上温异常情况一阶变化率曲线Fig.6 The curve of the polymerization reactor upper temperature abnormal situation vs the first order derivative

图7 聚合釜上温异常情况二阶变化率曲线Fig.7 The curve of the polymerization reactor upper temperature abnormal situation vs the second order derivative

4.3 结果分析

(1)依据现场操作经验，操作人员只能在聚合釜上温升到60℃时才能发现发生升温过快的异常情况。笔者所提出的趋势分析方法能够在第1080 s时监测出异常的发生，并及时发出报警，为现场操作提供了可靠的安全指导。

(2)为了验证所提方法的有效性，采用3σ法与所提方法进行对比。根据聚合釜上温升温阶段的10组历史数据，由3σ法求得，聚合釜上温服从均值为42.36，标准差为12.47的正态分布N(42.36，12.472)。利用3σ法对故障数据进行报警监测，采用[μ-2σ,μ+2σ]的报警阈值区间，但故障数据一直处于报警阈值内，在整个升温阶段都未能发现异常情况，如图8所示。由3σ法和所提方法的比较可以得出，3σ法在聚合釜上温出现升温过快的早期无法监测出异常情况并报警，而所提趋势分析方法在第1080 s时准确监测出异常的发生，提前24 min 34 s 发出报警。

图8 3σ法异常监测曲线Fig.8 The curve of the abnormal monitoring by the 3σ method

5 结 论

(1)针对现有的化工间歇过程报警研究没有考虑在阶跃变化过程中参数变化是否过快的问题，提出基于最小二乘的拟合-微分-再微分的趋势分析方法，能够准确识别间歇过程工况阶段，重点监测间歇过程中的阶跃变化阶段，实现参数变化异常早期报警。

(2)采用最小二乘法对过程数据进行拟合，通过拟合-微分之后过程参数的变化趋势进行间歇过程多时段工况识别，利用再微分对参数变化率的变化快慢分析，结合滑动窗算法，与报警阈值比较后，判断参数变化是否异常。

(3)将所提方法应用于聚丙烯聚合过程，准确识别出故障数据处于聚合过程升温阶段，相比于3σ法，能够提前24 min 34 s监测出异常情况的发生并发出报警。

符号说明：

a——拟合多项式系数；

D——滑动窗步长，s；

k——滑动窗更新次数；

L——滑动窗长度，s；

m——过程变量个数；

X——过程参数；

α——过程参数值；

β——参数一阶变化率值；

ε——参数二阶变化率值；

ν——参数一阶变化率；

γ——参数二阶变化率；

μ——正态分布均值；

σ——正态分布标准差。

[1] PENG Kaixiang, LI Qianqian, ZHANG Kai, et al. Quality-related process monitoring for dynamic non-Gaussian batch process with multi-phase using a new data-driven method[J].Neurocomputing, 2016,214: 317-328.

[2] 胡瑾秋， 张来斌， 王安琪. 基于格兰杰因果关系检验的炼化系统故障根原因诊断方法[J].石油学报(石油加工), 2016, 32(6): 1266-1272.(HU Jinqiu, ZHANG Laibin, WANG Anqi. Fault root-cause diagnosis based on Granger causality test for petrochemical process system[J].Acta Petrolei Sinica(Petroleum Processing Section), 2016, 32(6): 1266-1272.)

[3] 马曦， 张来斌， 胡瑾秋， 等. 基于IRML的油气加工系统多层次故障传播模型研究[J].石油学报(石油加工), 2015, 31(5): 1193-1202.(MA Xi, ZHANG Laibin, HU Jinqiu, et al. Hierachical fault propagation model for petroleum refining system based on IRML[J].Acta Petrolei Sinica(Petroleum Processing Section), 2015, 31(5): 1193-1202.)

[4] 蔡华伟. 聚合工艺重点参数的监控及安全控制[J].化工管理, 2014, (15): 186-188, 134.(CAI Huawei. Monitoring and safety control of key parameters in polymerization process[J].Chemical Management, 2014, (15):186-188, 134.)

[5] 郭小萍， 陆宁云， 高福荣， 等. 间歇过程滑动窗口子时段PCA建模和在线监测[J].控制与决策, 2005, 20(9): 1034-1037.(GUO Xiaoping, LU Ningyun, GAO Furong, et al. PCA modeling and on-line monitoring of sub-period of sliding window in batch process[J].Control and Decision, 2005, 20(9): 1034-1037.)

[6] 常玉清， 王姝， 谭帅， 等. 基于多时段MPCA模型的间歇过程监测方法研究[J].自动化学报, 2010, 36(9): 1312-1320.(CHANG Yuqing, WANG Shu, TAN Shuai, et al. Research on the monitoring method of batch process based on multi-period MPCA model[J]. Journal of Automation, 2010, 36(9): 1312-1320.)

[7] 王静， 胡益， 侍洪波. 基于GMM的间歇过程故障检测[J].自动化学报, 2015, 41(5): 899-905.(WANG Jing, HU Yi, SHI Hongbo. Fault detection based on GMM for batch processes[J].Journal of Automation, 2015, 41(5): 899-905.)

[8] 肖应旺， 徐保国. 基于ICA-MPCA的间歇过程监测方法[J].仪器仪表学报, 2009, 30(5): 990-996.(XIAO Yingwang, XU Baoguo. Monitoring method for batch process based on ICA-MPCA[J].Chinese Journal of Scientific Instrument, 2009, 30(5): 990-996.)

[9] 袁彬， 刘才学， 黄文楼， 等. 基于浮动报警阈值的燃料元件包壳破损监测方法的分析与研究[J].核动力工程, 2010, 31(1): 102-106.(YUAN Bin, LIU Caixue, HUANG Wenlou, et al. Analysis and research on fuel cladding failure monitoring based on floating alarm threshold[J].Nuclear Power Engineering, 2010, 31(1): 102-106.)

[10] 朱海， 高胜峰， 蔡鹏， 等. 基于拉依达准则的自适应小波阈值选取方法[J].海洋技术学报, 2016, 35(4): 50-54.(ZHU Hai, GAO Shengfeng, CAI Peng, et al. Study on the self-adaptive wavelet threshold selection method based on the Pauta criterion[J].Journal of Ocean Technology, 2016, 35(4): 50-54.)

[11] 田垅， 刘宗田. 最小二乘法分段直线拟合[J].计算机科学, 2012, 39(S1): 482-484.(TIAN Long, LIU Zongtian. Least-square method of piecewise linear fitting[J] .Computer Science, 2012, 39(S1): 482-484.)

[12] MARCO A, MARTNEZ J J.Polynomial least squarefitting in the Bernstein basis[J].Linear Algebra and Its Applications, 2010, 433(7): 1254-1264.

[13] 陈桂秀. 用程序求解最小二乘拟合多项式的系数[J].青海师范大学学报(自然科学版), 2010, 26(3): 14-17.(CHEN Guixiu. Solve the least square curve fitting polynomial coefficient with program[J].Journal of Qinghai Normal University (Natural Science), 2010, 26(3): 14-17.)

[14] 蔡华伟. 聚合工艺重点参数的监控及安全控制[J].化工管理, 2014, (15): 186-188, 134.(CAI Huawei. Monitoring and safety control of key parameters in polymerization process[J].Chemical Management, 2014, (15): 186-188, 134.)

[15] 林观炎， 乔建江， 吴达岭. 聚丙烯生产工艺危险性分析及安全措施[J].安防科技, 2012, (1): 24-27.(LIN Guanyan, QIAO Jianjiang, WU Daling. Hazard analysis and safety measures of polypropylene production process[J].Safety and Environmental Engineering, 2012, (1): 24-27.)