吕雪峰 谢耀滨

工业控制系统(Industrial control system)广泛应用于国家基础设施,工控系统一旦遭到破坏可能会造成难以估量的经济损失甚至人员伤亡.一个典型的工业控制网路可分为企业网络层、控制网络层和现场网络层,如图1所示.

近年来,由于IT系统的软件和硬件技术不断集成到工控领域,IT领域的一些漏洞和后门也出现在工控系统上.以智能制造为核心的工业4.0,在推动工业转型的同时也使得工控系统面临更多来自互联网的威胁[1].2010年,世界首个网络超级破坏武器“震网”病毒[2]被检测出来,证明了工控系统是可被攻击并被利用的.此后又相继爆发了“毒区”、“火焰”等病毒和一连串的工控系统入侵事件.工控信息安全受到越来越多的重视,针对工控系统入侵检测的研究也成为目前的一个热点.

由于工控保护机制的存在,大的破坏性攻击容易被检测出来,因此对入侵检测的研究越来越集中于隐蔽攻击[3−5].当前,工控系统入侵检测的研究方法主要包含三类:1)基于概率统计的方法[6−7];2)基于机器学习的方法[8−12];3)基于状态的方法[13−16].

基于状态的方法因为更好地考虑了工控系统的物理特性,检测准确度高,目前大部分研究都是基于此类方法[13].基于状态的方法通过历史数据和专家经验定义系统的临界状态,实时监控系统当前状态与临界状态的距离以判断系统是否处于危险之中.Fovino等[15]通过监控系统状态的变化来检测复杂攻击,创立了系统的虚拟镜像作为系统的内部表示,并用规则语言描述系统的临界状态,但是仅考虑了离散输入输出数据,没有考虑连续数据.Carcano等[17]进一步考虑了连续输入和输出,但在输入输出数据量庞大的情况下,很难定义系统的临界状态.针对这一问题,Khalili等[13]提出了一个系统化的解决方案SysDetect,该方案基于Apriori算法,通过专家经验的判定可显著减少算法下一次迭代产生的候选临界状态数.但是该方法依然不能很好地处理高维数据,且在一定程度上依赖人工判定.

图1 典型的工业控制系统网络架构Fig.1 Typical network architecture for industrial control system

针对当前工控入侵检测系统存在的问题,本文提出一种新的基于状态迁移图的异常检测方法.利用工控系统与物理世界交互的特性,采用数据驱动的方法进行建模,即利用系统运行时系统变量的数据来建立检测模型.但考虑到数据维度可能过高,不直接以系统变量数据来描述系统运行状态,而是以相邻数据向量间的余弦相似度和欧氏距离来表征,因而可以处理高维数据.状态迁移图刻画系统运行过程中的正常模型,根据正常历史数据样本训练得出,所以不需要依赖专家预先定义系统的临界状态.

1 恶意数据攻击建模

一个工业控制过程可简单地用图2来表示.其中传感器和执行器易成为攻击者的攻击目标,因为它们直接与物理过程相连,一旦遭到攻击将会产生不可估量的后果,例如“震网”病毒恶意篡改了伊朗核工厂离心机的转速,使核工厂被迫停工.

图2 工业控制过程:传感器和执行器易成为攻击目标Fig.2 Industrial control process:sensors and actuators are vulnerable targets

由于工控系统的物理特性,恶意数据注入成为一种简单且收益高的方式.攻击者通过攻击执行器和传感器,注入恶意数据,修改变量的值,对系统造成破坏.

系统变量分为传感变量和操纵变量,传感变量的值由传感器测量取得,而操纵变量则由控制器通过系统参数和传感变量值计算得出.记系统含l个测量变量,用集合表示,含p个操纵变量,用集合表示.则系统变量集合可用表示.根据攻击者试图破坏系统速度的快慢,将恶意数据注入攻击分为快速注入和慢速注入.

快速注入使得系统变量在短时间内发生较大改变,以最快速度造成破坏,表现为最大/最小值攻击.

其中,Ta为攻击时段.

慢速注入缓慢改变系统变量的值,以达到潜伏的目的,可以表现为偏置注入和几何注入.偏置注入连续注入多个较小的常量.

其中,ci表示常量,通常取值较小.几何注入逐渐增大系统变量的改变幅度,具有一定的潜伏性,同时不失破坏性,可表现为指数形式的注入.

其中,a,b可以是常数,也可以是变量,t0为常数.

以上只是列举了三种可能的攻击形式,实际上恶意数据注入攻击可以表现为更多的形式.

2 基于状态迁移图的异常检测方法

由于工控系统的诸多限制,例如有限的内存、有限的计算能力、较高的实时性要求等,工控入侵检测系统必须是轻量级的,检测规则或检测模型应设置得相对简单.本文采用一个状态迁移图的检测模型,状态迁移图的状态数与数据向量的维度无关,因而可以处理高维的数据向量.

2.1 状态表示

实际的工控系统中,状态变量可能有很多个,因此数据向量的维度可能很高,而现有的基于状态的入侵检测算法无法很好地处理数据维度较高的情况.考虑控制系统的运行模式,控制器根据前一时刻传感变量的值,计算出操纵变量的值,操纵变量的值再作用于物理系统使之发生变化,传感变量随着系统的运动而变化,在下一个采样周期将传感变量的值传给控制器,如此循环往复.通过分析这个过程可以发现,系统变量在当前时刻的取值很大程度上取决于上一时刻的取值,因此正常情况下,相邻两个数据向量间的变化不会太大.

相邻数据向量之间的变化反应了系统变量运行的时间特性,本文以余弦相似度和欧氏距离来衡量这种变化.余弦相似度反应数据向量方向的变化,常用于文档聚类和信息提取.欧氏距离衡量两个数据向量绝对距离的远近.

选择作为系统状态表征的特征应该是相互关联的,即相邻数据向量间的余弦相似度和欧氏距离应该是相互关联的.这样当余弦相似度在某个区间取值时,欧氏距离只位于某些特定区间,而不会在所有区间取值.如图3所示,假设余弦相似度取值区间为A,B,C,欧氏距离取值区间为D,E,F,二者在二维空间所有可能的区域为1～9,实际取值区域为1,4,5,9.当余弦相似度取值区间为B时,欧氏距离只在区间D,E取值,而不会在F区间取值,若欧氏距离在F区间取值就可以认为发生了异常.这样可以更加细致地刻画数据的轮廓,而不是将所有可能的区域为1～9都视为可接受的区域.余弦相似度和欧氏距离的关联性如图4所示,其中余弦相似度和欧氏距离采用Normal数据集计算得出(经过归一化处理),图4展现了前30组数据.从图4可以看出,余弦相似度和欧氏距离的变化基本相反,当余弦相似度在高位取值时,欧氏距离总在低位取值,证明了系统变量数据向量间的余弦相似度和欧氏距离是相互关联的.

图3 余弦相似度和欧氏距离取值示意Fig.3 Possible taking-value for cosine similarity and Euclidian distance

图4 余弦相似度和欧氏距离关联性示意Fig.4 Relevance between cosine similarity and Euclidian distance

本文以数据向量间的余弦相似度和欧氏距离组成的二元组来表征系统运行状态,避免了直接对数据进行建模,将高维的数据向量转化为二维数据,有效降低了计算复杂度.同时,考虑数据的时间特性,更加符合工控系统的动态特性.

2.2 状态矩阵求解

为了求解出状态迁移图,需要知道图的顶点和各顶点之间的转换关系.状态迁移图的顶点表示了系统可能的状态,本文用一个状态矩阵求解状态迁移图,需要知道图的顶点和各顶点之间的转换关系.状态迁移图的顶点表示了系统可能的状态,本文以状态矩阵STATE表示可能的状态,STATE的每一行称为一个状态点,一个状态点表示一个状态,对应状态迁移图中的一个顶点.

2.2.1 余弦相似度和欧氏距离的求解

令

cossim是所有训练样本点任意相邻两点间余弦相似度组成的列向量.相邻样本点间的欧氏距离由式(6)计算得出.

2.2.2 cossim和ed量化

eeeddd′的量化阈值集为

按照各量化区间点数相等的原则计算和.定义函数k:R→R,k(P)表示某实数集合P中元素的个数.定义函数

和的计算如下:

2.2.3 状态矩阵生成

2.3 状态迁移图生成

状态迁移图是一个有向图,作为异常检测模型,其顶点由STATE确定,边用邻接矩阵表示.

2.3.1 顶点确定

由于对相邻样本点间的余弦相似度和欧氏距离进行了量化,状态矩阵STATE中含有相同的状态点.相同的状态点对应状态迁移图中的同一个顶点,应将其合并,具体做法为:

步骤1.用自然数对STATE中的每个状态点进行标号,称为状态号,并用向量来表示,例如表示STATE中的第k个状态点对应的状态号.状态号唯一标明了状态迁移图中的状态(即顶点).

步骤2.令STATE第1个状态点的状态号为1,即,令状态数sn=1.

步骤3.从第2个状态点开始,对第k个状态点遍历STATE中该状态点之前的状态点,将该状态点与之前的状态点进行比较,若找到相同的状态点j则停止遍历,将该点的状态号设为第j个状态点的状态号,即;若没有找到相同的状态点,则令sn=sn+1,将该点的状态号设为sn,即.

经过步骤1～3后,STATE中相同的状态点具有相同的状态号,最终的状态数为sn.状态迁移图的顶点对应STATE中状态号为1～sn的状态点.

2.3.2 邻接矩阵确定

用NM表示邻接矩阵,邻接矩阵反映了状态迁移图中的状态转移关系,若NMi,j=1,则状态迁移图中的第j个顶点到第j个顶点可达;若NMi,j=0,则第i个顶点到第i个顶点不可达.由于样本点是随时间增长逐渐采样而来的,因此前一样本点到当前样本点总是可达的.相应地,STATE中前一状态点到当前状态点也总是可达的.根据这一原则可以求解出邻接矩阵NM,具体做法为:

步骤1.将NM置为全零.

步骤2.根据前一样本点(状态点)到当前样本点(状态点)可达的原则,令1))=1,由此计算出NM.

确定了顶点和邻接矩阵后,就能生成状态迁移图.

2.4 基于状态迁移图的在线检测

状态迁移图生成后,可以据此进行在线异常检测.基于状态迁移图的在线检测流程如图所示,具体步骤如下:

步骤1.获取当前样本点.

步骤2.计算当前样本点与上一样本的余弦相似度t_cossimi和欧氏距离t_edi.

步骤3.按照训练样本的标准对t_cossimi和t-edi进行归一化处理,得到qt_cossimi和qt_edi.

步骤4.遍历STATE,判断STATE中是否有和当前点相同的点.若存在,根据确定当前状态点对应的状态号.

步骤5.若当前状态点对应的状态不在状态迁移图内,产生告警,记录异常类别为“异常1”;若当前点对应的状态位于状态图内,即STATE中有和当前点相同的点,则根据状态迁移图判断上一节点对应状态到当前节点对应状态是否可达.若不可达,则产生告警,记录异常类别为“异常2”;若可达,继续下一个样本点的检测.

3 实验与测试

目前工控领域没有一个标准的测试集,而实际工控系统在攻击下的数据样本很难获取.针对这一情况,本文采用田纳西–伊斯曼(Tennesseeeastman,TE)过程[18]的MATLAB模型作为仿真平台.根据第1节建立的恶意数据攻击模型,选定斜坡注入和偏置注入两种攻击形式,在特定时刻对TE过程实施攻击,生成系统正常运行数据和攻击情形下的数据,并用第2.4节描述的检测方法进行检测.

图5 基于状态迁移图的在线检测流程Fig.5 Online detection process based on state transition graph

3.1 TE仿真模型

TE过程的原型由Downs和Vogel在1993年根据一个真实的化工过程创建.TE过程主要部件有冷凝器、反应器、循环压缩机、汽提塔和气液分离器等.该过程包含两个放热反应和两个副反应,生成产品G,H和副产品F.四种化学反应如下所示:

TE过程包含41个测量变量和12个操纵变量,是一个多变量,多数据,复杂的非线性控制系统,常用于多变量控制、非线性控制和故障诊断等领域.TE过程的MATALB模型由Ricker[19]给出,该模型采用文献[20]中的控制策略,仿真时间为48h,采样时间为3min,每次仿真共产生960组数据,每组数据中都包含高斯白噪声.TE过程的工艺流程如图6所示.

3.2 测试数据

反应器是TE过程的核心部件,是发生化学反应的地方,对温度的要求很高.因此本文选定攻击对象为反应器的温度传感器,并根据第1节所述攻击模型,选定偏置注入和斜坡注入两种攻击方式.

其中,偏置注入中c为常数,而斜坡注入在t0时注入斜率为k的斜坡信号.反应器温度控制回路如图7所示,两个输入分别代表反应器温度初始设定值(用输入点9表示)和反应器温度(用xmeas9表示),输出为反应器冷却水流量(用xmv10表示).该回路通过反应器温度测量值和反应器温度初始设定值计算冷却水流量,以此来控制反应器温度.为达到攻击反应器温度的目的,在反应器温度输入端口增加Add模块,将其与攻击信号叠加,攻击信号通过延时模块在特定时刻发挥作用.图8为加入攻击信号后的反应器温度控制回路,延时模块时间设置为20h,攻击信号为常数,代表偏置注入信号.攻击信号可以替换为其他信号,以此设计更多的攻击形式.

图7 反应器温度控制回路Fig.7 Control loop for reactor temperature

图8 加入攻击信号的反应器温度控制回路Fig.8 Control loop for reactor temperature added with attack signal

采集系统正常运行时各变量的数据,生成数据集合Normal.选定式(17)中c为0.01,0.1和1,在系统仿真20h加入注入攻击,生成偏置注入下的数据集Dataset1,Dataset2,Dataset3.选定式(18)中k为0.01,0.1和1,在仿真20h实施斜坡攻击,生成斜坡注入下的数据集Dataset4,Dataset5,Dataset6.各数据集及其对应参数见表1.这些数据集将用来测试本文所提的方法.

表1 测试数据集及其参数Table 1 The test data set and the corresponding parameters

为直观感受加入攻击信号后对反应器温度的影响,用MATLAB采集生成了反应器温度在各种情况下的温度–时间变化曲线,分别如图9～11所示.图9为正常工况条件下的温度–时间变化曲线,可以看出反应器温度稳定在123℃左右.图10为斜坡注入k=0.01情况下的温度变化情况,可以看到在20h后反应器温度呈线性下降趋势.图11为偏置注入c=0.1情况下的温度–时间变化曲线,反应器温度在20h突然下降,之后则稳定下来.

图9 正常工况下反应器温度随时间变化情况Fig.9 Reactor temperature varies with time undernormal condition

3.3 状态迁移图生成测试

由于状态迁移图通过正常数据训练得出,所以用表1中的Normal数据集来进行训练.根据第2.2节的分析,生成的状态迁移图与余弦相似度和欧氏距离的量化区间数intervals相关,因此本文不断改变intervals的取值,进行了多组状态图测试.首先按照第2.2节和第2.3节的方式对Normal数据集进行处理,然后设定量化区间数intervals,计算状态矩阵,最后生成状态迁移图.状态迁移图的顶点数和边数随intervals的变化情况如表2所示.

图10 斜坡注入工况(k=0.01)下,反应器温度随时间变化情况Fig.10 Reactor temperature varies with time under ramp signal injection withkset at 0.01

图11 偏置注入工况(c=0.1)下,反应器温度随时间变化情况Fig.11 Reactor temperature varies with time under bias signal injection withcset at 0.1

从表2可以看出,随着intervals的增大,状态迁移图的顶点数和边数也不断增大.而intervals越大,意味着将余弦相似度和欧氏距离的量化区间划分得越小,由于状态迁移图的检测由顶点和边来决定,所以顶点数和边数越多,意味着检测规则的粒度越细,相应地对异常会更加敏感,此外还会额外增加资源开销.图12和图13分别显示了intervals为5和8时用MATLAB生成的状态迁移图.

3.4 异常检测性能测试

按照第2.4节的在线检测步骤,用表1的测试数据对状态迁移图检测模型的检测效果进行测试,选定intervals为5,8,10,15,检测结果分别见表3～6.本文选用从攻击到正确检测到异常的时间以及误报率作为评价的标准.

表2 状态迁移图的顶点数和边数随intervals的变化Table 2 The nodes and triangles number of state transition graph varies withintervals

图12 intervals为5时的状态迁移图Fig.12 The state transition graph whenintervalsis equal to 5

图13 intervals为8时的状态迁移图Fig.13 The state transition graph whenintervalsis equal to 8

从检测时间来看,在各种攻击情况下检测模型检测到异常的时间随着intervals的变化而变化.从检测最坏的结果对应的数据集来看(对应的数据集为Dataset4,此时c=0,k=0.01),当intervals为5和8时候,分别在第436和第411个样本点检测到异常(在第20小时进行注入攻击,对应的第一个异常样本点为第401个).而当intervals为10和15时,检测模型在第401个样本点即检测到异常,说明当intervals增大时,检测模型对异常更加敏感,检测异常的速度有所提升.

表3 intervals=5时的检测结果Table 3 Detection results whenintervalsis equal to 5

表4 intervals=8时的检测结果Table 4 Detection results whenintervalsis equal to 8

表5 intervals=10时的检测结果Table 5 Detection results whenintervalsis equal to 10

表6 intervals=15时的检测结果Table 6 Detection results whenintervalsis equal to 15

从误报率来看,当intervals逐渐增大时,相应地误报率也随之增大.当intervals为5时误报率均不超过1%.当intervals增大到8时,误报率迅速增加到5% 左右.而当intervals取15时,误报率超过20%,这时可认为状态迁移图不能正确检测异常.

3.5 主元分析法测试

由于工控领域缺乏标准的数据集,研究人员都用各自的模型和数据进行工控系统异常的研究,很难复现,难以对比各种方法的好坏.为了进一步评估状态迁移图检测的性能,将其与控制系统常用的异常检测方法—主元分析法(Principal component analysis,PCA)进行比较.PCA采用T2统计量和SPE统计量作为检测异常的指标,当样本的T2统计量和SPE统计量超过各自的控制限时,检测到异常.选定状态迁移图模型检测结果最差的数据集Dataset4作为PCA方法的测试集,其测试结果如图14和图15所示.从图14和图15可以看出,PCA虽然最终能检测到异常,但是在第500个样本点才检测出来.相比于状态迁移图检测模型,检测到异常的速度较慢.

图14 PCA方法T2统计量Fig.14 T2statistic of PCA method

3.6 讨论

本节从时空资源消耗以及检测性能两方面对状态迁移图检测模型进行讨论,由于训练过程是通过历史数据离线进行的,因此时空消耗只考虑检测时的情况.

图15 PCA方法SPE统计量Fig.15 SPE statistic of PCA method

1)从时间消耗上来看,对新样本点的检测需要判断新样本点是否位于状态迁移图内,即需要遍历一次状态图,若状态迁移图含n个顶点,则时间复杂度为O(n).实际检测过程中,顶点数与设置的量化区间数intervals有关.通过第3.4节的实验结果不难发现,当intervals取5或8时已有较好的检测结果,继续增大intervals只会增大资源消耗,同时增加误报率.且intervals越大,需要越多的训练样本进行训练.以intervals取10为例,顶点数n最多为100,而实际会小于100,因此这一步骤需要的计算量很小.第二个步骤判断上一状态到当前状态是否可达,只需查询邻接矩阵中对应的数值即可,时间复杂度为O(1).

2)从空间消耗上来看,本文的异常检测方法需要存储量化阈值集合,状态矩阵STATE,邻接矩阵NM,其中邻接矩阵NM占用的存储空间为n的二阶次,,状态矩阵STATE存储消耗为n的一阶次.但由于n取值很小(不超过100),且邻接矩阵存储单位为比特,当intervals取10时,估计整体消耗的存储空间不超过1MB.

3)从检测性能来看,即便对于微小的攻击(对应数据集Dataset1和Dataset4),检测模型也能较快检测出异常,相比于常规的异常检测算法,例如第3.4节的PCA算法,状态迁移图检测模型能够更好地刻画系统运行的动态变化,更快地检测出异常.但是当intervals较大(取10以上)时,检测误报率较高.而当intervals较小(例如,取5时)时,不能快速地对微小攻击做出反应.因此需要谨慎选择intervals以平衡误报率和检测速度之间的关系.

4 总结

本文提出了一种新的基于状态迁移图的异常检测方法,详细描述了其状态表示、状态图生成和在线检测过程,建立了工控系统的恶意数据攻击模型,并基于TE仿真模型进行了各项实验,最后从时空资源消耗和检测性能两方面对方法进行了讨论.该方法具有较小的时空消耗,且考虑了系统运行的动态特性,能够快速检测出异常,因而适用于资源受限和对实时性要求较高的工业控制系统.