基于SARIMA的工业控制系统通信流量在线监测方法与流程

本发明涉及工业控制系统网络流量预测，尤其是涉及一种基于sarima的工业控制系统通信流量在线监测方法，属于工业信息安全检测领域。

背景技术：

能源、炼化和交通等关键基础设施是国家稳定运行的神经中枢，是我国网络安全的重中之重。随着国家大型基础设备(智能变电站，智能化工流程工业系统，工业分布式控制系统)的自动化，互联化以及智能化建设的推进，其网络空间安全问题日益凸显。近年来，一系列针对国家关键基础设施的网络攻击造成了极大的国民经济损失和对社会不可逆的破坏。这些顶尖黑客通过更隐蔽，更高效，杀伤力更大的入侵方式频繁入侵枢纽变电站，流程化工业系统甚至核电站的通信网络。目前，针对国家关键基础设施网络系统的防御与加固已经上升到国家战略层面，而通信流量分析则是工业系统安全问题公认最有前景的解决方案。通信流量智能化分析是将传统互联网领域的安全解决方案与现代化的电力通信网络和工业控制系统特性有机结合的交叉学科解决方案。

根据相关报告和文献，所有针对工业系统的攻击都会在通信网络上有所体现。大多数的工控网络攻击都会导致相关的通信网络受损，不同的攻击类型导致网络受损的程度和位置会有所不同。以“blackenergy”为主导的组合拳式的攻击以及一系列恶意代码注入会导致通信网络瘫痪，关键信道被阻塞，数据采集与监控(scada)系统被操纵，控制系统迟滞恢复与状态致盲等现象。由于ics中的数据流量表现出不同的流量模式和类似于互联网流量的特征，因此可以通过生成数学模型加以分析，开发和理解ics数据流量的特征。对于ics通信流量这种复杂时间序列，一般采用回归算法进行建模及统计学分析。现有的ics通信流量分析模型具有算法复杂程度高，对流量的解释性不强以及对ics通信流量模式的回归分析的准确性与预测精度不高的缺陷。且流量的检测方式一般为离线检测，不能对实时采集的ics通信流量进行实时分析与建模，同时也无法在线评估异常的程度，从而导致运维人员无法实现对ics网络状况的态势感知。

技术实现要素：

sarima模型是近年对通信网络流量进行监督学习的较为热门的算法，搭配分布式的在线检测算法，sarima模型可对实时采集的ics数据流量进行快速，精准地建模，从而对ics流量的动态变化进行分析与检测，最终实现对ics网络状况的实时态势感知。

本发明的目的在于解决在没有先验知识的前提下，对实时采集的ics通信流量的动态建模问题，针对现有ics通信流量异常检测算法过于依赖先验知识及算法复杂度高无法实际部署的不足提出了完善的分析方法；生成的动态ics通信流量阈值模型对国家重大工业基础设施的网络安全防护及异常检测具有指导意义。

本发明的目的可以通过以下技术方案来实现：

一种基于sarima的工业控制系统通信流量在线监测方法，其包括如下步骤：

1)在ics通信网络之中部署工业交换机、监控主机以及测试主机，监控主机从工业交换机中实时采集通信流量数据；

2)根据选定的流量聚合尺度，建立小周期sarima(p,d,q)x(p,d,q)s模型；

3)根据不同的置信度生成ics动态流量阈值模型；

4)对实时采集的流量序列进行动态分析。

所述的步骤3)具体为：

3.1)采集ics工业交换机上的实时流量数据，根据设定的采样频率γsamp，聚合尺度生成时间序列，并以一个小周期为一个迭代周期；

3.2)对采集到的实时流量数据进行训练，在一个小周期内进行sarima(p,d,q)x(p,d,q)s模型训练与适配之后，输出最优的模型以及模型适配的参数；第i个小周期的模型定义为：

其中fsarima()为sarima(p,d,q)x(p,d,q)s模型的函数表达式，为第i次迭代的小周期训练集，tfore为小周期预测的序列个数，s为周期性参数，'bic'为计量经济学选取最优的(p,d,q,p,d,q)参数的准则，是sarima模型第i次迭代预测出的时间序列；

计算第i次迭代的预测均值

3.3)分布式运行小周期的sarima(p,d,q)x(p,d,q)s模型，对实时采集的第i个小周期的ics流量数据进行实时动态滚动建模，此时的实时流量数据相当于验证集；将验证集与预测过程生成的基于置信区间的流量阈值上下限进行对比分析；其定义的第i个小周期的上下界阈值为：

其中是以1-αp.i为显著性的z分布值，为第i个小周期的动态阈值区间的上界，为动态阈值区间的下界，均为长度为tfore的时间序列，αp.i为置信度；

实时采集的通信流量以及极低的算法复杂度保证了与在时间上的重合性，第i次迭代的正常ics通信流量(时间序列内的每一个时序样本)定义为：

算法实质上是比较动态采集流量的时间序列与两个长度相同的阈值时间序列与所对应的每一个时序样本的大小关系；其中为第i次小周期实时采集的ics通信流量，tfore是采集的样本量；

第i次小周期的ics预测序列与训练序列存在实时的相关性，第i次小周期的预测序列可估计为：

其中函数∩为对两个时间序列集合取交集。

3.4)流量判定结束之后，继续下一个小周期的训练迭代，并重新输出新的最优模型以及模型适配的参数，对新输入的实时流量数据进行再次判定；

3.5)循环整个过程，直到达到设定的迭代次数。

其中，第i次迭代的异常序列的时间序列判定计算方法如下(默认时间序列的计算方式为对时间序列中一一对应的每一个时序样本进行计算)：

其中每一次迭代的小周期待检测序列都需要通过分布式运行的sarima(p,d,q)x(p,d,q)s模型，当该小周期待检测序列触发异常检测条件时，生成异常矩阵当流量的值处于正常区间范围时其中

假设ics异常发生的现实对应时间为序列其中n为异常发生的总次数。则异常发生时小周期序列号的样本数目序列为：

其中tdebug为程序运行之前的实时调试时间，γsamp为ics流量的采样频率。

第n次ics异常发生在第次小周期迭代计算之中，因此ics异常事件发生的小周期迭代次数(的元素)可以由以下的方程组追溯计算得到：

其中kn为中间变量。

sarima在线检测算法所生成的动态ics流量阈值区间的方差可以由以下的算法得到：

其中k＝1,2，…，n能对第次小周期迭代的阈值整体偏差进行衡量，对分析ics网络异常程度具有指导意义。

第i次迭代的ics网络流量异常程度的计算方法为：

其中当元素时为轻微的ics流量异常，为显著的ics流量异常，当时ics流量异常的程度最高，此时对ics网络具有毁灭性的影响。

本发明的有益效果在于解决了ics通信流量的动态建模以及异常检测问题；生成的ics动态通信流量阈值模型对ics通信网络的网络安全防护具有指导意义。ics通信流量阈值模型能实时监测流量走向，快速扫描全网并为日常网络维护提供实时、精准的网络流量流向和流量成分的分析，为未来网络优化、网络调整、网络建设提供决策支持的数据依据。本发明对实际采集的不同工控场景下的通信网络流量进行精确可靠的实时动态建模，在流量层面上针对不同的网络攻击进行智能化的异常检测，能满足对典型ics的入侵检测与安全态势感知，从而高效准确地抵御针对ics的典型网络攻击。

本发明针对ics通信流量的特点，结合现有sarima模型建模的结论，对实时采集的ics通信流量进行分布式sarima模型建模，并动态生成实时的ics通信网络流量的阈值模型。本发明实地采集浙江省某虚实结合的工控靶场ics通信流量，并进行建模分析，建立具有适当参数的流量模型，验证ics的网络流量特点。利用计量经济学的思想，使用优化的sarima(p,d,q)x(p,d,q)s模型建立ics的正常流量阈值模型，并分析模型切合程度。基于自定义、可适配的在线检测算法，对实时采集的工控流量进行动态建模分析，生成动态阈值区间以区分异常。分析不同置信度情况下的流量阈值模型，并与ics通信特征加以结合。本发明最终在浙江省某化工集团的过程工业控制系统中进行实际部署应用，具有较高的检测率与较低的误检率。

附图说明

图1是本发明的方法流程图；

图2是本发明的算法的单个小周期简明示意图；

图3是本发明的在线算法的多个小周期迭代的简明示意图；

图4是原始ics通信流量以及不同聚合尺度下的聚合流量；

图5是原始ics网络流量差分及自相关性分析图；

图6是聚合尺度为10s的ics网络流量1阶差分及自相关性图；

图7是聚合尺度为60s的ics网络流量1阶差分及自相关性图；

图8是聚合尺度为60s的ics网络流量季节性分析图；

图9(a)是sarima(2,1,3)x(2,0,0)10模型的拟合情况图；

图9(b)是sarima(2,1,3)x(2,0,0)10模型的预测效果图；

图10是不同的置信区间对ics流量的阈值区间的影响图；

图11(a)是sarima(1,0,1)x(0,0,1)6模型拟合情况如图；

图11(b)是sarima(1,0,1)x(0,0,1)6模型的预测效果及异常检测图；

图12是不同的ttrai/tfore比例对算法延迟时间以及最优bic值的影响情况。

具体实施方式

下面根据附图详细说明本发明，本发明的目的和效果将变得更加明显。图1是本发明的方法流程图；图2是本发明的算法的单个小周期简明示意图；图3是本发明的在线算法的多个小周期迭代的简明示意图。

实验前期采集了浙江大学某虚实结合的ics靶场的通信网络流量。实验室配置了工业plc控制器，工业以太网交换机以及工控上位主机。其中上位机与plc之间采用tcp/ip的通信协议。plc到现场设备层之间采用工业类modbus协议。对实际ics通信网络流量进行采集并存储，对流量的特性进行离线分析。流量探针部署的位置在工控上位机与控制器之间的工业交换机上，分析流量的类型为单个交换机口的局部流量以及交换机镜像口的全流量。实验对镜像口全流量数据进行建模与数据分析，并同时对局部流量进行相关的实验。在无任何人为操作与干扰的情况下，本次实验一共抓取了正常运行的ics通信网络中的111815个数据包，大约运行时长为16小时，无丢包现象；并选择了几种典型的工业控制网络通信的场景进行分析。其原始流量以及不同聚合尺度下的聚合流量如图4所示。

对实际采集得到的原始ics网络流量进行差分以及自相关性分析，如图5所示。由图5可以看出，原始序列在时序上并不平稳，对原始序列进行1阶差分之后的序列已经接近平稳，2阶差分之后的序列需要在强假设的前提下才能近似平稳。原始序列、1阶差分序列、2阶差分序列的自相关函数(acf)都较为近似；可以看出其acf值在较小的滞后项(lag)内均超过划定的显著性水平临界值，后续的滞后项也有超过临界值的现象。由此分析得原始序列需要经过1阶差分才能达到较好的arima(p,d,q)建模要求，但由于存在滞后项超过临界值的现象，需要对原始序列进行周期性分解，即需要对序列进行sarima(p,d,q)x(p,d,q)s建模。

对实际采集得到的聚合尺度为10s的ics网络流量进行1阶差分以及自相关性分析，如图6所示。聚合后(以10s为一个流量聚合尺度)的ics网络流量的波动程度大于原始流量，且其自相关函数的滞后项(lag)在整个时间尺度上都有超过显著水平临界值的现象。相较而言1阶差分之后的聚合序列较为平稳，自相关函数的滞后项在整个时间尺度上较少出现显著超过临界值的现象(1-2阶滞后项除外)。

对实际采集得到的聚合尺度为60s的ics网络流量进行1阶差分以及自相关性分析，如图7所示。聚合后(以60s为一个流量聚合尺度)的ics网络流量的波动程度大于原始流量，其自相关函数的滞后项(lag)在整个时间尺度上较小。相较而言1阶差分之后的聚合序列较为平稳，但局部有存在突增与突减；自相关函数的滞后项在整个时间尺度上几乎不出现显著超过临界值的现象(1-2阶滞后项除外)。对聚合后(以60s为一个流量聚合尺度)的ics网络流量数据进行季节性分析，如图8所示。可以看出，其聚合序列可以季节性拆分为趋势分量以及季节性分量。因此分析得到ics通信网络流量具有一定的周期性。

参照图1，图2和图3，本发明基于sarima的工业控制系统通信流量在线监测方法包括如下步骤：

1)在ics通信网络之中部署工业交换机、监控主机以及测试主机，监控主机从工业交换机中实时采集通信流量数据；

2)根据选定的流量聚合尺度，建立小周期sarima(p,d,q)x(p,d,q)s模型；

2.1)定义选定的流量聚合尺度以及小周期分析尺度，用sarima(p,d,q)x(p,d,q)s序列的定义法来产生sarima(p,d,q)x(p,d,q)s的时间序列：

sarima(p,d,q)x(p,d,q)s模型是通过对arma(p,q)模型分别进行d阶差分以及d阶季节性差分计算得到的，arma(p,q)模型是由ar(p)以及ma(q)模型组合而成；

自回归滑动平均模型arma(p,q)定义如下：

xt＝φ1xt-1+φ2xt-2+…+φpxt-p+εt-θ1εt-1-…-θqεt-q

上式：xt为均值化处理之后的小周期的平稳时间序列，其时间序列的长度较短；φp为自回归项ar的系数；θq为滑动平均项ma的系数；εt为随机扰动；p为ar的阶数；q为ma的阶数；

定义一个延迟算子b，bxt＝xt-1，则ar系数多项式φ(b)＝1-φ1b-…-φp(b)^p，ma系数多项式θ(b)＝1-θ1b-…-θq(b)^q；

引入差分算子δ^d＝(1-b)^d，则arima(p,d,q)模型表示为：

φ(b)δ^dxt＝θ(b)εt

sarima模型通过对arima模型进行季节性差分运算得到，sarima模型的定义如下：

φp(b)φp(b^s)δ^dδs^dxt＝θq(b)θq(b^s)εt

其中εt为白噪声序列，d为趋势差分的阶数，d为周期s为补偿的季节差分阶数，b^s为s阶延迟算子，δs^d为季节性差分算子；b^sxt＝xt-s，δs^d＝1-b^s，φp(b^s)为b^s的q阶多项式，φp(b^s)为b^s的p阶多项式；

2.2)使用贝叶斯信息准则bic对sarima(p,d,q)x(p,d,q)s模型的p,d,q,p,d,q,阶数进行监督分析并定阶。bic信息标准是根据观测数据选择最佳参数模型的信息标准，bic定义如下：

bic＝kin(n)-2in(l)

＝(模型复杂度的惩罚项)-(模型的精度优势项)

在上述表达式的右边，第一项表示模型的复杂性，第二项反映拟合的优劣；

2.3)采用最小二乘法对sarima(p,d,q)x(p,d,q)s的p阶系数φk(k＝1,2,…,p)，q阶系数θk(k＝1,2,…,q)，以及季节性p阶系数季节性q阶系数进行估计；

2.4)采用最优bic准则下的sarima(p,d,q)x(p,d,q)s模型对原始序列进行拟合分析，并进行残差检验；若残差为白噪声，则对拟合序列进行反滤波处理，得到原序列的拟合值或预测值；若残差不为白噪声，则重新采用bic信息准则对arma(p,q)模型进行定阶；

2.5)得到小周期的sarima(p,d,q)x(p,d,q)s的数学表达式。

3)根据不同的置信度生成ics动态流量阈值模型；

3.1)采集ics工业交换机上的实时流量数据，根据设定的采样频率γsamp，聚合尺度生成时间序列，并以一个小周期为一个迭代周期；

3.2)对采集到的实时流量数据进行训练，如图2所示，在一个小周期内进行sarima(p,d,q)x(p,d,q)s模型训练与适配之后，输出最优的模型以及模型适配的参数；第i个小周期的模型定义为：

其中fsarima()为sarima(p,d,q)x(p,d,q)s模型的函数表达式，为第i次迭代的小周期训练集，tfore为小周期预测的序列个数，s为周期性参数，'bic'为计量经济学选取最优的(p,d,q,p,d,q)参数的准则，是sarima模型第i次迭代预测出的时间序列；

计算第i次迭代的预测均值

3.3)如图3所示，分布式运行小周期的sarima(p,d,q)x(p,d,q)s模型，对实时采集的第i个小周期的ics流量数据进行实时动态滚动建模，此时的实时流量数据相当于验证集；将验证集与预测过程生成的基于置信区间的流量阈值上下限进行对比分析；其定义的第i个小周期的上下界阈值为：

其中是以1-αp.i为显著性的z分布值，为第i个小周期的动态阈值区间的上界，为动态阈值区间的下界，均为长度为tfore的时间序列，αp.i为置信度；

第i次迭代的正常ics通信流量定义为：

其中为第i次小周期实时采集的ics通信流量，tfore是采集的样本量；

第i次小周期的ics预测序列与训练序列存在实时的相关性，第i次小周期的预测序列可估计为：

其中函数∩为对两个时间序列集合取交集。

3.4)流量判定结束之后，继续下一个小周期的训练迭代，并重新输出新的最优模型以及模型适配的参数，对新输入的实时流量数据进行再次判定；

3.5)循环整个过程，直到达到设定的迭代次数。

在步骤3)中，第i次迭代的异常序列的时间序列判定计算方法如下：

其中每一次迭代的小周期待检测序列都需要通过分布式运行的sarima(p,d,q)x(p,d,q)s模型，当该小周期待检测序列触发异常检测条件时，生成异常矩阵当流量的值处于正常区间范围时其中

假设ics异常发生的现实对应时间为序列其中n为异常发生的总次数。则异常发生时小周期序列号的样本数目序列为：

其中tdebug为程序运行之前的实时调试时间，γsamp为ics流量的采样频率。

第n次ics异常发生在第次小周期迭代计算之中，因此ics异常事件发生的小周期迭代次数(的元素)可以由以下的方程组追溯计算得到：

其中kn为中间变量。

sarima在线检测算法所生成的动态ics流量阈值区间的方差可以由以下的算法得到：

其中k＝1,2，…，n能对第次小周期迭代的阈值整体偏差进行衡量，对分析ics网络异常程度具有指导意义。

第i次迭代的ics网络流量异常程度的计算方法为：

其中当元素时为轻微的ics流量异常，为显著的ics流量异常，当时ics流量异常的程度最高，此时对ics网络具有毁灭性的影响。

4)对实时采集的流量序列进行动态分析。

在对浙江大学某虚实结合的ics靶场上实际采集的通信网络流量进行前期离线分析的基础上，针对ics通信网络流量的特征，进一步在靶场上部署的工业交换机，监控主机及测试主机上进行在线测试分析。在原靶场环境中添加现场测试机(corei5,macbookpro)进行在线建模与异常分析。现场测试机通过网线连接到交换机上，运行编译完毕的python脚本实时采集工业以太网交换机内的全流量，并动态建模分析，生成流量阈值区间对采集的ics网络流量进行实时识别与预警。定义的ics流量的采样频率γsamp为1(s)。

对ics通信网络流量进行sarima(p,d,q)x(p,d,q)s建模。定义小周期训练项数为ttrai＝300，tfore＝30，并对单个小周期的sarima(p,d,q)x(p,d,q)s模型进行计量经济学分析，假定本次小周期训练为第一次迭代，即i＝1。

假设对原始序列不进行聚合，对在线采集的单个小周期的ics流量数据进行建模，得到模型sarima(2,1,3)x(2,0,0)10。其计量经济学参数如下表所示：

表1

其sarima(2,1,3)x(2,0,0)10模型的参数拟合表如下所示：

表2

sarima(2,1,3)x(2,0,0)10模型的拟合情况如图9(a)所示，模型的预测情况如图9(b)所示。其模型的拟合参数r-square＝0.89。由此可得sarima(2,1,3)x(2,0,0)10对该小周期ttrai＝300的ics通信流量具有较强的解释性，其阈值区间代表了tfore＝30的ics通信流量的情况。由图9(b)中的阈值区间以及验证序列可知，其时序上的ics通信流量序列符合训练集的流量特征与规律，因此并未超过阈值区间的上下界。阈值区间的上界时间序列为：

其时间序列所对应的时间戳为：

{13:52:19,13:52:20,13:52:21,…,13:52:50,13:52:51}，该时间戳表明ics通信流量与工业现场实际环境相互吻合，物理信息空间上相互对应。

因此在对小周期ttrai＝300的ics通信流量进行sarima(2,1,3)x(2,0,0)10建模后，其验证周期tfore＝30的ics通信流量不存在异常情况，属于正常运行的ics通信流量。ics异常发生的现实对应时间为序列为空集，异常事件产生时的序列小周期迭代次数序列为空集。此时的方差为0。其异常程度序列

不同的置信区间对流量的阈值区间有不同的影响，其影响情况如图10所示。置信系数越大，所对应的阈值区间越大，系统对异常流量的敏感性越低，因此异常检测系统的误报率较低；置信系数越小，所对应的阈值区间越小，系统对异常流量的敏感性较高，更容易发现小规模的流量异常，但是系统的误报率较高。

对下一个不同采样频率的ics流量序列进行实时异常检测，定义的ics流量的采样频率γsamp为60(s)。此时修改聚合尺度为60s，原始序列的时间跨度增加。此时取得sarima(1,0,1)x(0,0,1)6为最优的模型。模型的运行起始时间为04:14:26，并维持流量在线监测模型的运行状态。截取其中某一次小周期运行模型参数。

sarima(1,0,1)x(0,0,1)6模型的经济学参数如下表3所示：

表3

sarima(1,0,1)x(0,0,1)6模型的拟合参数如下表4所示：

表4

sarima(1,0,1)x(0,0,1)6模型拟合情况如图11(a)所示，模型的预测情况如图11(b)所示。此时选择的ics流量验证序列为异常流量序列，由测试主机对工业交换机进行tcp-flooding攻击得到；该异常流量序列在时序上的体现为单位时间内流量会存在突增的现象，攻击过后流量会恢复正常。tcp-flooding攻击于14:36:50注入，流量突增，随后于14:38:00停止注入异常流量。

此时阈值区间的上界时间序列为：

此时阈值区间的下界时间序列为：

其时间序列所对应的时间戳为：

{14:36:20,14:37:20,14:38:20,…,14:54:20,14:55:20}

通过上述模型，我们可以近似得到较短时间间隔内某实际时刻的正常流量阈值。例如在时刻14:36:20，当置信区间为95％时其正常流量的区间为[350.3,1650.3]，当置信区间为90％时其正常流量的区间为[398.3,1621.6]。由此我们可以得到正常情况下ics某时刻通信流量的阈值模型。此时异常事件发生的实时时间为：

其中n为3，表示存在三处异常。

tdebug＝4，γsamp＝60，可得

由公式

可知ics异常事件发生的小周期迭代次数均为17。

通过计算分析，此时存在第16个小周期的训练序列在时间上与第17个小周期保持近似连续。因此该ics通信流量在线检测方法的实时性得到了很好的保证。

sarima在线检测算法所生成的动态ics流量阈值区间的方差可以由以下的算法得到：

(其中k＝1,2，3)能对第次小周期迭代的阈值整体偏差进行衡量，对分析ics网络异常程度具有指导意义。从上述方差分析可知其正常流量的正常波动范围区间为0～1475.079。

由上述分析可知，ics异常事件发生时的小周期迭代次数为17，另外此时的流量异常发生的时刻为14:37。由于流量是按照1分钟为聚合尺度进行聚合的，因此在14:37-14:38这个时段，异常检测判定为异常，其流量超过其阈值上限。

因此ics网络流量验证集的前两个序列均存在超过阈值上限的情况。此次迭代中的ics通信网络异常参数为：因此在14:36:20-14:38:20这个时段出现显著的ics流量异常的状况，14:38:20-14:39:20出现轻微的ics流量异常的状况。ics流量异常程度序列后续的流量不存在异常情况。

到此第17个迭代小周期结束，依照定义可以进行下一个迭代小周期的ics通信网络流量分析。

图12为15次迭代小周期分析之中不同的ttrai/tfore比例对算法延迟时间以及最优bic值的影响情况，其中100/30的延迟时间比率最低，500/30的延迟比率相较100/30较大一些，但是其算法的精度比100/30高。不同训练尺度下的bic值相差不大，证明算法有较好的收敛性。

不同的网络异常所反映的ics流量模式有所差别，在后续的异常检测中可以根据不同的侧重点对流量模型的非正常情况做分析，并引入机器学习算法对异常流量进行分类。