基于ARIMA-LSTM模型的核电厂运行事件变化分析及预测

侯秦脉 朱 伟 邹 象 刘时贤 吴彦农

（生态环境部核与辐射安全中心 北京 102445）

《中国的核安全》首次发布于2019年，该书是中国政府发表的第一部核安全白皮书，白皮书全面介绍了中国核安全事业的发展历程，阐述中国核安全的基本原则和政策主张，分享中国核安全监管的理念和实践，阐明中国推进全球核安全治理进程的决心和行动［1］。截至2019年6月，已安全稳定运行累计300余堆年，未发生过国际核与放射事件分级表（International Nuclear Event Scale，INES）2级及以上的事件或事故，且0级偏差和1级异常事件发生率呈下降趋势。该报告提供了1991～2018年我国运行核电厂的运行事件数量，同时《核安全年报》也提供了2019～2020年的运行事件数量［2］，此领域研究目前多为对运行事件本身的分析和一般的趋势性描述。郑丽馨等［3］指出运行事件对了解机组性能、运行经验趋势的评价有重要的贡献。胡攀［4］认为对低级别事件趋势分析能够跟踪管理异常趋势并提高异常趋势的识别效率。刘时贤等［5］从误碰方面分析了运行事件对核电厂的安全性和经济性造成影响。修炳林等［6］指出2015～2016年新机组投运高峰带来运行事件的增加。但缺乏有效的工具对运行事件数量变化趋势及预测进行研究，因此，需要研究一套方法进行核电厂运行事件数量的趋势分析和预测。

差分自回归移动平均模型（Autoregressive Mobile Average Model，ARIMA）是经典的时间序列预测方法，武云云等［7］探讨了运用ARIMA模型在核电站外围环境放射性水平动态分析和预测，为核电站外围环境影响评价和核事故放射性污染提供参考。胡衍坤等［8］指出ARIMA模型能够较好地体现时间序列数据中的线性特征，但是单一的ARIMA模型对非线性变化难以充分有效地处理，结合网络结构模型能有效提高河流水质预测精度7%。在深度学习算法中，长短期记忆神经网络模型（Long and Short-Term Memory，LSTM）由于其特殊的网络结构，在处理事件序列问题时，比传统的神经网络更快更容易收敛到最优解，非常适合处理核电厂运行事件指标这种时序数据。因此，本文建立ARIMA与LSTM组合模型进行核电厂运行事件预测的研究，通过Mann-Kendall检验法（简称M-K检验法）探讨运行事件变化特性以及预测的合理性，并结合白皮书中我国运行核电厂运行事件数据，分析运行事件变化趋势，为核安全监管提供数据支持。

1 研究方法描述

1.1 M-K检验法

M-K检验法是一种非参数检验方法［9］，对样本分布无要求，能有效排除个别异常值对时间序列的影响，该方法广泛应用于大气和地质等行业时间序列突变点检测研究中［10］。在突变分析方面，M-K检验能够查找出数据异常点，确定突变区域。

根据核电厂运行事件数时间序列的秩序列为：

在时间序列随机独立的假定下，定义统计量：

式中：E（Sk）、Var（Sk）分别是秩序列Sk的均值和方差。

根据UBk=−UFk，同样方法计算反序列得到统计量UBk。依据UBk和UFk可得到UB曲线和UF曲线，差值大于零则呈下降趋势，小于零呈上升趋势，差值等于零为相交点［11−12］，如果相交点在0.05显著性水平下，则该点是突变点。

1.2 ARIMA模型

ARIMA是一种广泛用于时间序列分析和预测的模型，对处理线性时间序列有较好的优势［13］。非平稳时间序列通过该方法差分后转化为平稳时间序列，模型对因变量滞后值以及其随机误差进行回归，融合了两种模型：自回归模型（Autoregressive，AR）和移动平均模型（Moving Average，MA）。组间随机变量具有的依存关系，既包括本身的变化规律同时也包括了外部影响因素，该方法有较高的预测精度，有助于解释预测变动的规律。

首选需要确保预测时，ARIMA模型的时间序列是平稳的或近似平稳的，自相关系数只有时间间隔一个变量，即时间变化不影响均数和方差。基于同质性是非平稳时间序列的普遍特征，可以通过差分处理转化为平稳时间序列，差分次数就是其次的阶数。那么完整的差分自回归移动平均过程ARIMA（p，d，q）模型可以用滞后算子B表示为［14］：

式中：p是自回归滞后阶；d是非平稳序列转化为平稳序列的差分次数；q是移动平均阶数，偏自相关函数值显著不为零的最高阶数；ei是i时刻的偏差；Фi是自回归部分的系数；B是滞后算子；θi是移动平均部分的系数［15］。

ARIMA模型建立方法如下：

1）模型识别，验证时间序列数据平稳性，非平稳数据需进行差分处理，利用自相关函数图（Autocorrelation Function，ACF）和偏自相关函数图（Partial Autocorrelation Function，PACF）进行分析，确定自回归滞后阶数、差分次数和移动平均阶数。

2）参数估计和模型检验。用特定时间序列的数据，估计模型的参数并进行检验［16］，判断该模型是否恰当，如不恰当，返回步骤1）。

3）趋势预测应用，并对设定的将来某个时间段给出预测值。

1.3 LSTM模型

LSTM是一种递归神经网络（Recursive Neural Network，RNNs）的改进［17］，可以解决长时间周期依赖的预测，并且对于有限数据样本深入分析其长期趋势，解决RNNs距离感知能力对梯度的依赖，在各行业时间序列预测上LSTM有着广泛的应用［18］。

LSTM模型在递归神经网络的基础上添加了遗忘、输入、输出等三个门层，解决了有限时间和数据样本的问题［19］，通过梯度下降得到模型的最佳误差参数。LSTM的第一步是确定遗忘门层遗忘哪些信息，借助Sigmoid层决定前一时间周期的节点数据。读取上一层节点数据输出ht−1和输入数据Xt，并为单元状态Ct−1输出0~1的值ft，0和1分别表示“不允许任何量通过”和“允许任何量通过”的逻辑意义［20］。ft通过Sigmoid函数的输出值得到，分配给Ct−1决定前一个状态的值是否丢弃，ft为0表示“遗忘”，1表示“保留”，Wf是权重矩阵，bf是偏置项，其表达式为：

下一步是在输入层中决定输入哪些信息。Sigmoid函数决定将更新哪些值，添加到新候选值Ct的tanh函数输入值（−1~1）组成一个向量，并乘以Sigmoid函数值，得到输出值。其表达式为：

最后需要通过输出门层输出函数的参数。首先运行一个Sigmoid层，将得到输出参数信息ot，包含决定输出的单元状态信息［21］，最终得到将要输出的部分ht，其表达式为：

1.4 组合预测模型

核电厂运行事件数据具有线性特征和非线性特征，将运行事件数据记为At，分解为线性和非线性两个部分，分别用Bt表示线性部分和Ct表示非线性部分，表达式如下：

第一步，时间序列模型建模；对时间数据序列At进行建模，然后将训练数据进行拟合并进行预测，通过自相关图和偏自相关图建立ARIMA（p，d，q）模型，输出时间数据序列Bt'，其数据序列描述线性部分，得到偏差et的时间数据序列。

第二步，LSTM残差预测；使用LSTM模型对偏差的时间数据序列et进行拟合。对数据序列At中非线性部分数据信息进行挖掘，识别Ct中与et的异常关系，通过对et进行LSTM处理，得到预测结果Ct'。

第三步，模型组合；通过ARIMA（p，d，q）模型表达线性部分和LSTM模型表达非线性部分建立混合模型，得到混合预测模型表达式。

因此，将刻画数据序列线性关系的ARIMA（p，d，q）模型和描述数据序列非线性关系的LSTM模型的各自优点有机结合，能够显示混合模型的优越性和提高预测的准确度。

2 实验数据和变化特征分析

2.1 实验数据及趋势分析

实验运行事件数量数据来源于2019年9月国务院新闻办公室发表的《中国的核安全》白皮书，白皮书提供了1991～2018年我国运行核电厂的运行事件数量。国家核安全年报提供了2019～2020年的运行事件数量，将用作组合模型预测值的对比。

采用线性回归的方法对1991～2018年运行事件数量时间序列进行线性趋势分析（图1），运行事件数量倾向率为线性方程的斜率，在0.95置信限下运行事件数量趋势函数为f（x）=69.252−0.017 5x。由图1可以看出，1）运行事件数量与时间呈负线性相关趋势。2）年度运行事件波动性变化明显，波动范围在8~93，1993年、2002年和2016年达到波峰，1991年、2000年和2009年达到波谷，经历了枯-丰-正常-枯-丰-正常-枯-丰阶段，最大值与最小值相差10倍。

图1 1991～2018年运行事件数量及趋势Fig.1 Number and trend of operation events from 1991 to 2018

2.2 运行事件数量时间序列突变分析

采用M-K方法检验突变点，对1991～2018年我国运行核电厂的运行事件数量进行检验，检验结果见图2。结果表明：年度运行事件的突变检验曲线存在交汇点，图2中时间序列M-K突变检验在0.05显著性水平下，UF统计量与UB统计量从1995～2013年无交点即无突变点，这一时期我国核电机组数量以年均不到1台机组缓慢发展，运行事件数量无显著突变。

图2 运行事件数量M-K检验结果Fig.2 M-K test results of the number of operation events

而在1991年、1994年、2014年、2016年和2017年存在交点，说明1991～1994年以及2014～2018年为运行事件序列显著突变年份，查阅核电机组数量从1991～1994年核电机组从无到有的摸索阶段，同时2014～2018年我国核电机组数量年均增长4.2台，对于核电行业来说这两个阶段都是快速发展期，处在新人员积累经验的阶段，可以说明突变点伴随着行业的高速发展，同时也说明机组数量年均增长1台不会引起运行事件数量突变，而机组数量年均增长4台会引起国内运行事件数量的突变，应对方法将是在高速发展阶段前提前储备人才，加强人员职业能力培养，核安全监管队伍规模也应该做出相应调整。

3 运行事件时间序列预测

从文献［8］中可以看出，客观事物具有线性和非线性特征，本文§2趋势分析和突变分析反映了运行事件也具有线性和非线性的双重特征。对于运行事件的线性关系，通过R软件调用Tseries和Forecast程序，实现ARIMA模型对运行事件数量的预测。对于运行事件的非线性关系，通过R软件调用基于Tensorflow框架的keras程序，实现LSTM模型的残差预测。

采用组合模型按照5个步骤对我国核电厂运行事件数量建模并预测：

1）平稳化处理1991～2018年运行事件数量时间序列。首先通过ADF检测检测未差分时间序列是否平稳，未差分的p-value=0.22，大于0.05，可知此时间序列为非平稳序列；经过一阶差分后p-value达到0.05，可知1阶差分后时间序列平稳化，确定数值1作为ARIMA模型参数d的值。

2）ARIMA模型参数p、q定阶。通过一阶差分后运行事件数量时间序列的自相关图和偏自相关图，如图3所示。图3（a）中说明，有3个自相关系数显著不为0，参数q可以推断取3或2（考虑3是偶然出现的因素，则取2是显著不为0的阶数），图3（b）中说明，有2个偏自相关系数显著不为0，参数p可以推断取2或1（考虑2是偶然出现的因素，则取1是显著不为0的阶数）。

图3 一阶差分后的自相关图ACF(a)和偏自相关图PACF(b)Fig.3 ACF(a)and PACF(b)maps of first order difference

根据参数p、q初步定阶结果，ARIMA模型可暂定为ARIMA（3，1，2）、ARIMA（3，1，1）、ARIMA（2，1，2）、ARIMA（2，1，1），等待赤池检验确定最终模型参数。由赤池信息准则（Akaike Information Criterion，AIC）可知，AIC结果大小决定着模型优劣，数值越小表明模型越佳，模型参数估计见表1。从表1中AIC值可知，ARIMA（2，1，2）为运行事件预测的最佳模型。

表1 ARIMA模型参数估计Table 1 ARIMA model parameter estimation

3）白噪声检验。通过计算Ljung-Box检验统计量，用于检验给定时间序列中独立的零假设。检验ARIMA（2，1，2）模型残差的白噪声，残差的白噪声检验统计量结果p-value=0.998 3，大于0.05，说明模型残差显著性检验没有通过，确定了残差无自相关，ARIMA（2，1，2）模型能够较好地拟合核电厂运行事件数量时间序列。

4）ARIMA模型预测。本文采用ARIMA（2，1，2）模型对1991～2018年运行事件数数量进行预测，并预测了在95%置信限值下的2019～2021年运行事件数量预测结果，如表2所示。

表2 2019～2021年运行事件数量预测Table 2 Predicted values of operation events from 2019 to 2021

5）LSTM残差预测，合并得到组合模型预测值。采用LSTM模型对ARIMA模型预测运行事件数量残差进行训练和预测，训练集和测试集比例分别占数据集的2/3和1/3，对于训练集控制训练迭代50次，训练损失随迭代次数快速下降，在20次后趋于平稳，说明拟合良好。计算得到运行事件数量残差的预测结果，将ARIMA运行事件数量预测结果与LSTM运行事件数量残差预测结果相加得到组合模型预测结果，组合模型预测结果如图4所示。ARIMA预测值相对真实值滞后的问题在时域分析是普遍存在的，如果运行事件是一个平滑序列的预测，出现滞后的概率就会偏低。如果在新的时刻运行事件出现了一个跳变被认为是客观规律，不能判断为异常，那么预测将会滞后。从图4看出，2019～2021年运行事件趋势稳中向下，有利于减轻运行板块的核安全监管压力，实现在不增员的前提下突出监管重点，科学协调审评和监督的资源向在建板块倾斜，有利于板块轮动和提升监管能力。

图4 组合模型运行事件数量预测结果Fig.4 Predicted values of combined model

4 模型评价

根据ARIMA模型和ARIMA-LSTM组合模型预测的运行事件数量结果，分别计算模型预测运行事件数量的均方误差（Mean-Square Error，MSE）、根误差（Root Mean Square Error，RMSE）和平均绝对百分比误差（Mean Absolute Percentage Error，MAPE）（表3），可用于量化模型预测的准确性。

由表3可以看出，ARIMA（2，1，2）模型运行事件预测值的MSE、RMSE、MAPE分别为234.92、15.32和36.5%，ARIMA-LSTM组合模型的MSE、RMSE、MAPE分 别 为223.94、14.96和33.6%。ARIMA（2，1，2）-LSTM组合模型比ARIMA（2，1，2）模型运行事件数量的预测精度提高约3%。可以看出，ARIMA模型能够较好地体现时间序列数据中的线性特征，但是单一的ARIMA模型对非线性变化难以充分有效的处理，胡衍坤等结合网络结构模型能有效提高河流水质预测精度7%，与本文运行事件数量的预测精度提高能力相近，在具有更多固定成分的不稳定时间序列上表明LSTM做得更好，随着未来运行数据的增加，预测精度将会更加显著。

表3 模型评价参数Table 3 Model evaluation parameters

根据国家核安全局发布的《国家核安全局2019年报》和《国家核安全局2020年报》，可以得到2019～2020年的核电厂运行事件数量分别为31起和20起，与组合模型预测值的2019年27.64和2020年21.08非常接近，证明该组合模型能够精确地预测核电厂运行事件数量，能够很好地描述运行事件数量时间序列的趋势。有利于提前调配在建板块和运行板块的核安全监管审评力量，依据预测趋势变化大小，在趋势下降时调配资源到在建板块，在趋势上升时适当从在建板块调配资源到运行板块，前瞻性地调配有限监管力量，实现板块间合理轮动。

5 结语

本文研究利用M-K检验法对《核安全白皮书》中1991～2018年核电厂运行事件数量进行了趋势和突变分析，利用ARIMA（2，1，2）模型对核电厂运行事件数量进行预测，并利用LSTM模型对非线性部分的残差进行预测，结合两模型预测值得到最终组合模型的预测值。通过与2019～2020年《国家核安全局年报》核电厂运行事件数量对比和验证，其预测结果能够准确地描述实测值，表明其预测结果具有良好的可信度。

ARIMA-LSTM组合模型比ARIMA单一模型的预测预测精度提高了3%，充分发挥了LSTM模型和ARIMA模型的优点，前者优化非线性部分，后者拟合线性部分，高效地处理了核电厂运行事件数量存在的非线性问题。混合模型对应用差分自然回归的ARIMA模型的残差进行优化，有效解决了单一模型存在的误差问题。研究结果表明，通过对原始数据的优化，组合模型能够提高核电厂运行事件数量预测的精度。

作者贡献声明侯秦脉：提出研究思路、设计研究方案、进行数据计算、数据分析及论文的撰写；朱伟：负责完善研究方案、数据收集、参与数据计算、稿件的审阅与修订；邹象：负责提供技术支持与指导；刘时贤：负责研究进度的监督与指导；吴彦农：负责收集参考文献。