CN112187528B - Industrial control system communication flow online monitoring method based on SARIMA - Google Patents

Abstract

The invention discloses an Industrial Control System (ICS) communication flow online monitoring method based on SARIMA. The method is implemented by carrying out SARIMA (P, D, Q) x (P, D, Q) with small period on the flow sequence of the industrial control network_sModeling analysis; and generating a flow threshold model based on different confidence degrees by using the training and prediction step length and the confidence interval defined by the small period. The monitoring host collects communication flow data from the industrial switch in real time, runs a plurality of SARIMA models in a distributed mode in a small period, generates a real-time threshold interval, and performs anomaly detection analysis on the current communication flow collected in real time. The invention carries out experimental analysis by using a target range test board combining industrial control safety, deficiency and reality in Zhejiang province, and gives detailed algorithm description for test data; finally, the method is deployed and applied on the spot in a certain chemical group in Zhejiang province, and the reliability and accuracy of the algorithm are verified.

Description

Industrial control system communication flow online monitoring method based on SARIMA

Technical Field

The invention relates to network flow prediction of an industrial control system, in particular to an online monitoring method for communication flow of the industrial control system based on SARIMA, and belongs to the field of industrial information safety detection.

Background

The key infrastructure of energy, refining, transportation and the like is the national stable operation nerve center and is the central importance of network safety in China. With the promotion of automation, interconnection and intelligent construction of national large-scale basic equipment (intelligent substations, intelligent chemical process industrial systems and industrial distributed control systems), the network space safety problem is increasingly highlighted. In recent years, a series of network attacks against the national key infrastructure have caused great national economic losses and irreversible damage to society. The top hackers frequently invade communication networks of hub substations, flow industrial systems and even nuclear power stations in a more concealed, more efficient and more powerful invasion mode. At present, defense and reinforcement aiming at national key infrastructure network systems have risen to the national strategic level, and communication traffic analysis is the most promising solution for the safety problem of industrial systems. The intelligent analysis of communication flow is a cross-disciplinary solution which organically combines the security solution in the traditional internet field with the characteristics of a modern power communication network and an industrial control system.

According to the relevant reports and literature, all attacks against industrial systems are reflected on the communication network. Most industrial control network attacks can cause damage to related communication networks, and the degree and the position of the damage to the networks can be different due to different attack types. Combined punch type attacks and a series of malicious code injections which are mainly 'Blackenergy' can cause the phenomena that a communication network is paralyzed, a key channel is blocked, a data acquisition and monitoring (SCADA) system is operated, and a control system delays recovery and state blindness occur. Since data traffic in ICS exhibits different traffic patterns and characteristics similar to internet traffic, the characteristics of ICS data traffic can be analyzed, developed and understood by generating mathematical models. For the ICS communication flow such a complex time series, a regression algorithm is generally adopted for modeling and statistical analysis. The existing ICS communication flow analysis model has the defects of high algorithm complexity, low interpretability on flow and low accuracy and prediction precision of regression analysis of an ICS communication flow model. And the flow detection mode is generally off-line detection, real-time analysis and modeling cannot be performed on ICS communication flow collected in real time, and meanwhile, the abnormal degree cannot be evaluated on line, so that operation and maintenance personnel cannot realize situation perception of the ICS network condition.

Disclosure of Invention

The SARIMA model is a popular algorithm for supervising and learning communication network flow in recent years, is matched with a distributed online detection algorithm, and can be used for quickly and accurately modeling ICS data flow acquired in real time, so that dynamic change of the ICS flow is analyzed and detected, and real-time situation perception of ICS network conditions is finally realized.

The invention aims to solve the problem of dynamic modeling of ICS communication flow acquired in real time on the premise of no prior knowledge, and provides a perfect analysis method aiming at the defects that the existing ICS communication flow abnormity detection algorithm depends too much on the prior knowledge and cannot be practically deployed due to high algorithm complexity; the generated dynamic ICS communication flow threshold model has guiding significance for network security protection and anomaly detection of national major industrial infrastructure.

The purpose of the invention can be realized by the following technical scheme:

an industrial control system communication flow online monitoring method based on SARIMA comprises the following steps:

1) deploying an industrial switch, a monitoring host and a testing host in an ICS communication network, wherein the monitoring host collects communication flow data from the industrial switch in real time;

2) establishing a small-period SARIMA (P, D, Q) x (P, D, Q) s model according to the selected flow aggregation scale;

3) generating an ICS dynamic flow threshold model according to different confidence degrees;

4) and dynamically analyzing the flow sequence acquired in real time.

The step 3) is specifically as follows:

3.1) collecting real-time flow data on an ICS industrial exchanger according to a set sampling frequency gamma_sampGenerating a time sequence by a polymerization scale, and taking a small period as an iteration period;

3.2) training the collected real-time flow data, and carrying out SARIMA (P, D, Q) x (P, D, Q) in a small period_sAfter the model is trained and adapted, outputting an optimal model and parameters for model adaptation; the model for the ith small period is defined as:

wherein f is_SARIMA() Is SARIMA (P, D, Q) x (P, D, Q)_sThe functional expression of the model is shown as,

small period training set, T, for the ith iteration_foreThe number of sequences predicted for a small period, s is a periodicity parameter,'BIC' is a criterion for choosing optimal (P, D, Q, P, D, Q) parameters for the economics of metrology,

is a time sequence predicted by the ith iteration of the SARIMA model;

calculating the predicted mean value of the ith iteration

3.3) distributed operation of short-period SARIMA (P, D, Q) x (P, D, Q)_sThe model is used for carrying out real-time dynamic rolling modeling on ICS flow data of the ith small period collected in real time, wherein the real-time flow data at the moment is equivalent to a verification set; comparing and analyzing the verification set with the upper and lower flow threshold limits based on the confidence interval generated in the prediction process; the upper and lower threshold values of the defined ith small period are as follows:

wherein

Is at 1-alpha_P.IFor the value of the z-distribution of significance,

is the upper bound of the dynamic threshold interval for the ith small period,

is the lower bound of the dynamic threshold interval,

are all of length T_foreTime series of (a)_P.IIs the confidence level;

the communication flow collected in real time and the extremely low algorithm complexity ensure

And

the temporal coincidence, the normal ICS communication traffic (each time-sequential sample within the time series) for the ith iteration is defined as:

the algorithm is essentially a time series comparison of dynamically acquired flows

Two threshold time series with same length

And

the magnitude relation of each corresponding time sequence sample; wherein

ICS communication traffic, T, collected in real time for the ith small cycle_foreIs the amount of sample collected;

ICS prediction sequence of i-th minor cycle

And training sequence

With real-time correlation, i-th minor period of predictionSequencing

Can be estimated as:

wherein the function # is the intersection of the two time series sets.

3.4) after the flow judgment is finished, continuing training iteration of the next small period, outputting a new optimal model and model adaptive parameters again, and judging newly input real-time flow data again;

3.5) the whole process is circulated until the set iteration number is reached.

Wherein the abnormal sequence of the ith iteration

The time series judgment calculation method of (1) is as follows (the calculation mode of the default time series is to calculate each time series sample corresponding to one in the time series):

wherein each iteration of the small-period sequence to be detected needs to pass through distributed operation SARIMA (P, D, Q) x (P, D, Q)_sA model for generating an anomaly matrix when the small-period sequence to be detected triggers an anomaly detection condition

When the value of the flow rate is in the normal interval range

Wherein

Assume that an ICS anomaly has occurredReal time to time as a sequence

Where n is the total number of occurrences of the anomaly. The sequence of the number of samples of the small period sequence number when the abnormality occurs

Comprises the following steps:

wherein t is_debugFor real-time debugging time before program run, gamma_sampIs the sampling frequency of the ICS traffic.

The nth ICS abnormality occurs at

The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs

(

Can be calculated retrospectively from the following system of equations:

wherein K_nIs an intermediate variable.

Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm

The following algorithm can be used:

where k is 1,2, …, n

Can be aligned with

And the integral threshold deviation of the minor-cycle iteration is measured, and the method has guiding significance for analyzing the abnormal degree of the ICS network.

Degree of ICS network traffic anomaly of ith iteration

The calculation method comprises the following steps:

wherein when an element is present

A slight ICS traffic anomaly in time,

for significant ICS traffic anomalies, when

The ICS flow is abnormal to the highest degree, and the ICS network is affected destructively at the moment.

The invention has the advantages that the problems of dynamic modeling and abnormal detection of ICS communication flow are solved; the generated ICS dynamic communication flow threshold model has guiding significance on network security protection of an ICS communication network. The ICS communication flow threshold model can monitor flow trend in real time, quickly scan the whole network, provide real-time and accurate analysis of network flow direction and flow components for daily network maintenance, and provide data basis for decision support for future network optimization, network adjustment and network construction. The invention carries out accurate, reliable and real-time dynamic modeling on the actually acquired communication network flow under different industrial control scenes, carries out intelligent anomaly detection aiming at different network attacks on the flow level, and can meet the intrusion detection and the security situation perception of typical ICS, thereby efficiently and accurately resisting the typical network attack aiming at the ICS.

Aiming at the characteristics of ICS communication flow, the distributed SARIMA model modeling is carried out on the ICS communication flow collected in real time by combining the conclusion of the existing SARIMA model modeling, and a real-time threshold model of the ICS communication network flow is dynamically generated. The method provided by the invention is used for collecting the ICS communication flow of the industrial control target range combining virtuality and reality in Zhejiang province on the spot, carrying out modeling analysis, establishing a flow model with appropriate parameters and verifying the network flow characteristics of the ICS. Using the idea of metrological economics, using optimized SARIMA (P, D, Q) x (P, D, Q)_sThe model establishes an ICS normal flow threshold model and analyzes the model fitting degree. And based on a self-defined and adaptable online detection algorithm, dynamic modeling analysis is carried out on the industrial control flow collected in real time, and a dynamic threshold interval is generated to distinguish abnormity. And analyzing the flow threshold models under different confidence degrees and combining the flow threshold models with ICS communication characteristics. The invention is finally applied to actual deployment in a process industrial control system of a chemical group in Zhejiang province, and has higher detection rate and lower false detection rate.

Drawings

FIG. 1 is a flow chart of a method of the present invention;

FIG. 2 is a simplified schematic of a single small cycle of the algorithm of the present invention;

FIG. 3 is a simplified schematic of a plurality of small-cycle iterations of the online algorithm of the present invention;

FIG. 4 is the original ICS communication traffic and the aggregate traffic at different aggregation scales;

FIG. 5 is a diagram of the original ICS network traffic differential and autocorrelation analysis;

FIG. 6 is a graph of ICS network traffic 1 order difference and autocorrelation with an aggregation scale of 10 s;

FIG. 7 is a graph of ICS network traffic 1 order difference and autocorrelation with an aggregation scale of 60 s;

FIG. 8 is a diagram of seasonal analysis of ICS network traffic at an aggregation scale of 60 s;

FIG. 9(a) is SARIMA (2,1,3) x (2,0,0)₁₀Fitting situation graph of the model;

FIG. 9(b) is SARIMA (2,1,3) x (2,0,0)₁₀A prediction effect graph of the model;

FIG. 10 is a graph of the effect of different confidence intervals on the threshold interval of ICS traffic;

FIG. 11(a) shows SARIMA (1,0,1) x (0,0,1)₆The model fitting situation is as shown in the figure;

FIG. 11(b) is SARIMA (1,0,1) x (0,0,1)₆A prediction effect and abnormality detection map of the model;

FIG. 12 is a variation of T_trai/T_foreThe influence of the ratio on the algorithm delay time and the optimal BIC value.

Detailed Description

The objects and effects of the present invention will become more apparent from the following detailed description of the present invention with reference to the accompanying drawings. FIG. 1 is a flow chart of a method of the present invention; FIG. 2 is a simplified schematic of a single small cycle of the algorithm of the present invention; FIG. 3 is a simplified schematic of a plurality of small-cycle iterations of the online algorithm of the present invention.

The communication network flow of an ICS target field combined with the deficiency and the excess of Zhejiang university is collected in the early stage of the experiment. An industrial PLC controller, an industrial Ethernet switch and an industrial control upper host are configured in a laboratory. And a TCP/IP communication protocol is adopted between the upper computer and the PLC. An industrial Modbus protocol is adopted between the PLC and the field device layer. And collecting and storing the actual ICS communication network flow, and performing offline analysis on the flow characteristics. The arrangement position of the flow probe is arranged on an industrial exchanger between an industrial control upper computer and a controller, and the analysis flow type is the local flow of a single exchanger port and the full flow of an exchanger mirror image port. The experiment carries out modeling and data analysis on the full flow data of the mirror port, and simultaneously carries out related experiments on the local flow. Under the condition of no manual operation and no interference, the experiment captures 111815 data packets in the ICS communication network which normally runs, the running time is about 16 hours, and no packet loss phenomenon exists; and several scenarios of typical industrial control network communication are selected for analysis. The original flow rate and the polymerization flow rate at different polymerization scales are shown in FIG. 4.

The actually acquired original ICS network traffic is subjected to difference and autocorrelation analysis, as shown in fig. 5. As can be seen from fig. 5, the original sequence is not stable in time sequence, the sequence after 1 st order of difference of the original sequence is already close to stable, and the sequence after 2 nd order of difference needs to be stable under strong assumption. The autocorrelation function (ACF) of the original sequence, the 1 st order differential sequence and the 2 nd order differential sequence are all approximate; it can be seen that the ACF value exceeds the predetermined significance level threshold within the smaller lag term (lag), and the subsequent lag terms also exceed the threshold. The original sequence obtained by the analysis needs to be subjected to 1-order difference to achieve better ARIMA (P, D, Q) modeling requirements, but due to the fact that the lag term exceeds a critical value, the original sequence needs to be periodically decomposed, namely, the sequence needs to be subjected to SARIMA (P, D, Q) x (P, D, Q)_sAnd (6) modeling.

And performing 1-order difference and autocorrelation analysis on the ICS network traffic with the aggregation scale of 10s, which is obtained by actual acquisition, as shown in FIG. 6. The fluctuation degree of the ICS network flow after aggregation (taking 10s as a flow aggregation scale) is larger than the original flow, and the lag term (lag) of the autocorrelation function of the ICS network flow exceeds a significant level critical value on the whole time scale. The hysteresis term of the autocorrelation function is less likely to significantly exceed the threshold over the entire time scale (except for the 1-2 order hysteresis term) as compared to the stationary aggregation sequence after the 1 st order difference.

And performing 1-order difference and autocorrelation analysis on the ICS network traffic with the aggregation scale of 60s, which is obtained by actual acquisition, as shown in FIG. 7. The ICS network flow after aggregation (taking 60s as a flow aggregation scale) fluctuates to a greater extent than the original flow, and the lag term (lag) of the autocorrelation function is smaller on the whole time scale. Compared with the 1 st order difference, the polymerization sequence is more stable, but the local part has a sudden increase and a sudden decrease; the lag term of the autocorrelation function hardly appears to significantly exceed the critical value over the entire time scale (except for the 1-2 order lag term). The traffic data of the ICS network after aggregation (60 s is a traffic aggregation scale) is subjected to seasonal analysis, as shown in fig. 8. It can be seen that its aggregate sequence can be split seasonally into a trend component as well as a seasonal component. Therefore, the flow of the ICS communication network obtained through analysis has certain periodicity.

Referring to fig. 1, fig. 2 and fig. 3, the online monitoring method for communication flow of an industrial control system based on SARIMA of the present invention comprises the following steps:

1) deploying an industrial switch, a monitoring host and a testing host in an ICS communication network, wherein the monitoring host collects communication flow data from the industrial switch in real time;

2) establishing a small-period SARIMA (P, D, Q) x (P, D, Q) s model according to the selected flow aggregation scale;

2.1) defining selected flow aggregation scales and small period analysis scales using SARIMA (P, D, Q) x (P, D, Q)_sSequence definition to generate SARIMA (P, D, Q) x (P, D, Q)_sTime series of (a):

the SARIMA (P, D, Q) x (P, D, Q) s model is obtained by respectively carrying out D-order difference calculation and D-order seasonal difference calculation on an ARMA (P, Q) model, and the ARMA (P, Q) model is formed by combining an AR (P) model and an MA (Q) model;

the autoregressive moving average model ARMA (p, q) is defined as follows:

X_t＝φ₁X_t-1+φ₂X_t-2+…+φ_pX_t-p+ε_t-θ₁ε_t-1-…-θ_qε_t-q

the formula is as follows: x_tThe time sequence is a stable time sequence of a small period after the equalization processing, and the length of the time sequence is shorter; phi is a_pIs the coefficient of the autoregressive term AR; theta_qCoefficient of the moving average term MA; epsilon_tIs random disturbance; p is the order of AR; q is the order of MA;

defining a delay operator B, BX_t＝X_t-1If the polynomial of the AR coefficient phi (B) is 1-phi₁B-…-φ_p(B)^pThe MA coefficient polynomial Θ (B) is 1- θ₁B-…-θ_q(B)^q；

Introduction of a difference operator Δ^d＝(1-B)^dThen the ARIMA (p, d, q) model is expressed as:

Φ(B)Δ^dX_t＝Θ(B)ε_t

the SARIMA model is obtained by carrying out seasonal difference operation on the ARIMA model, and is defined as follows:

Φ_p(B)Φ_P(B^s)Δ^dΔ_s ^DX_t＝Θ_q(B)Θ_Q(B^s)ε_t

wherein epsilon_tIs a white noise sequence, D is the order of trend difference, D is the seasonal difference order with period s as compensation, B^sFor delay operators of order s, Δ_s ^DIs a seasonal difference operator; b is^sX_t＝X_t-s，Δ_s ^D＝1-B^s，Φ_P(B^s) Is B^sPolynomial of order Q,. phi_P(B^s) Is B^sA polynomial of order P;

2.2) pairs of SARIMA (P, D, Q) x (P, D, Q) using Bayesian information criterion BIC_sAnd (5) carrying out supervision analysis and order determination on the P, D, Q, P, D and Q orders of the model. The BIC information standard is an information standard for selecting an optimal parameter model according to observation data, and the BIC is defined as follows:

BIC＝kIn(n)-2In(L)

the model complexity penalty term (precision advantage term of model)

On the right of the above expression, the first term represents the complexity of the model, and the second term reflects the goodness of the fit;

2.3) applying least squares method to SARIMA (P, D, Q) x (P, D, Q)_sP-order coefficient of (phi)_k( k 1,2, …, p), q-order coefficient θ_k(k-1, 2, …, q), and seasonal P-th order coefficient

Seasonal Q-order coefficient

Carrying out estimation;

2.4) SARIMA (P, D, Q) x (P, D, Q) using the optimal BIC criterion_sThe model performs fitting analysis on the original sequence and performs residual error detection; if the residual error is white noise, performing inverse filtering processing on the fitting sequence to obtain a fitting value or a predicted value of the original sequence; if the residual error is not white noise, the BIC information criterion is adopted again to carry out order fixing on the ARMA (p, q) model;

2.5) obtaining a small period of SARIMA (P, D, Q) x (P, D, Q)_sThe mathematical expression of (1).

3) Generating an ICS dynamic flow threshold model according to different confidence degrees;

3.1) collecting real-time flow data on an ICS industrial exchanger according to a set sampling frequency gamma_sampGenerating a time sequence by a polymerization scale, and taking a small period as an iteration period;

3.2) training the collected real-time flow data, as shown in FIG. 2, SARIMA (P, D, Q) x (P, D, Q) is performed in a small period_sAfter the model is trained and adapted, outputting an optimal model and parameters for model adaptation; the model for the ith small period is defined as:

wherein f is_SARIMA() Is SARIMA (P, D, Q) x (P, D, Q)_sThe functional expression of the model is shown as,

small period training set, T, for the ith iteration_foreThe number of sequences predicted for a small period, s is a periodic parameter, 'BIC' is a criterion for selecting optimal (P, D, Q, P, D, Q) parameters for the economics of metrology,

is a time sequence predicted by the ith iteration of the SARIMA model;

calculating the predicted mean value of the ith iteration

3.3) distributed operation of SARIMA (P, D, Q) x (P, D, Q) for a short period as shown in FIG. 3_sThe model is used for carrying out real-time dynamic rolling modeling on ICS flow data of the ith small period collected in real time, wherein the real-time flow data at the moment is equivalent to a verification set; comparing and analyzing the verification set with the upper and lower flow threshold limits based on the confidence interval generated in the prediction process; the upper and lower threshold values of the defined ith small period are as follows:

wherein

Is at 1-alpha_P.IFor the value of the z-distribution of significance,

is the upper bound of the dynamic threshold interval for the ith small period,

is the lower bound of the dynamic threshold interval,

are all of length T_foreTime series of (a)_P.IIs the confidence level;

the normal ICS traffic for the ith iteration is defined as:

wherein

ICS communication traffic, T, collected in real time for the ith small cycle_foreIs the amount of sample collected;

ICS prediction sequence of i-th minor cycle

And training sequence

There is a real-time correlation, the predicted sequence of the ith minor period

Can be estimated as:

wherein the function # is the intersection of the two time series sets.

3.4) after the flow judgment is finished, continuing training iteration of the next small period, outputting a new optimal model and model adaptive parameters again, and judging newly input real-time flow data again;

3.5) the whole process is circulated until the set iteration number is reached.

In step 3), the exception sequence of the ith iteration

The time series judgment calculation method of (2) is as follows:

wherein each iteration of the small-period sequence to be detected needs to pass through the distributionRunning SARIMA (P, D, Q) x (P, D, Q)_sA model for generating an anomaly matrix when the small-period sequence to be detected triggers an anomaly detection condition

When the value of the flow rate is in the normal interval range

Wherein

Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence

Where n is the total number of occurrences of the anomaly. The sequence of the number of samples of the small period sequence number when the abnormality occurs

Comprises the following steps:

wherein t is_debugFor real-time debugging time before program run, gamma_sampIs the sampling frequency of the ICS traffic.

The nth ICS abnormality occurs at

The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs

(

Can be calculated retrospectively from the following system of equations:

wherein K_nIs an intermediate variable.

Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm

The following algorithm can be used:

where k is 1,2, …, n

Can be aligned with

And the integral threshold deviation of the minor-cycle iteration is measured, and the method has guiding significance for analyzing the abnormal degree of the ICS network.

Degree of ICS network traffic anomaly of ith iteration

The calculation method comprises the following steps:

wherein when an element is present

A slight ICS traffic anomaly in time,

for significant ICS traffic anomalies, when

The ICS flow is abnormal to the highest degree, and the ICS network is affected destructively at the moment.

4) And dynamically analyzing the flow sequence acquired in real time.

On the basis of carrying out early-stage off-line analysis on communication network flow actually acquired on an ICS target site combined with virtuality and reality of Zhejiang university, aiming at the characteristics of the ICS communication network flow, an industrial switch, a monitoring host and a testing host which are deployed on the target site are further subjected to on-line test analysis. An in-situ testing machine (Core i5, Macbook pro) was added to the in-situ target field environment for on-line modeling and anomaly analysis. The field test machine is connected to the switch through a network cable, runs the compiled Python script to acquire the full flow in the industrial Ethernet switch in real time, dynamically models and analyzes, generates a flow threshold interval and carries out real-time identification and early warning on the acquired ICS network flow. Sampling frequency gamma of defined ICS flow_sampIs 1(s).

SARIMA (P, D, Q) x (P, D, Q) for ICS communication network traffic_sAnd (6) modeling. Defining the small period training item number as T_trai＝300，T_fore30 and for a single small period of SARIMA (P, D, Q) x (P, D, Q)_sThe model performs a metric-economic analysis, assuming this small-period training as the first iteration, i.e., i is 1.

Modeling the ICS flow data of single small period collected on line to obtain a model SARIMA (2,1,3) x (2,0,0) under the assumption that the original sequence is not aggregated₁₀. The parameters of the economics of the measurement are shown in the following table:

TABLE 1

Its SARIMA (2,1,3) x (2,0,0)₁₀The table of the parameter fit of the model is as follows:

TABLE 2

SARIMA(2,1,3)x(2,0,0)₁₀The fitting of the model is shown in fig. 9(a), and the prediction of the model is shown in fig. 9 (b). Fitting parameters of its model

R-Square ═ 0.89. SARIMA (2,1,3) x (2,0,0) is thus obtainable₁₀For the small period T_traiThe ICS traffic of 300 has strong explanatory property, and the threshold interval represents T_fore30 ICS communication traffic. As can be seen from the threshold interval and the verification sequence in fig. 9(b), the ICS communication traffic sequence in the time sequence thereof conforms to the traffic characteristics and rules of the training set, and therefore does not exceed the upper and lower bounds of the threshold interval. The upper bound time series of the threshold intervals is:

the time stamps corresponding to the time sequence are as follows:

{13:52:19,13:52:20,13:52:21, …,13:52:50,13:52:51}, the time stamps indicating the ICS communication traffic and the actual environment of the industrial site coincide with each other, and the physical information spaces correspond to each other.

Thus in the period T _trai300 ICS traffic for SARIMA (2,1,3) x (2,0,0)₁₀After modeling, it verifies the period T_foreThe ICS communication flow of 30 has no abnormal condition and belongs to the ICS communication flow in normal operation. The realistic correspondence time of occurrence of ICS anomaly is sequence

As an empty set, the sequence of the number of iterations of the sequence in a small period when an abnormal event occurs

Is an empty set. Variance at this time

Is 0. Sequence of degree of abnormality thereof

Different confidence intervals have different effects on the threshold interval of the flow, the effect of which is shown in fig. 10. The larger the confidence coefficient is, the larger the corresponding threshold interval is, the lower the sensitivity of the system to abnormal flow is, so that the false alarm rate of the abnormal detection system is lower; the smaller the confidence coefficient is, the smaller the corresponding threshold interval is, the higher the sensitivity of the system to abnormal flow is, the easier the small-scale flow abnormality is found, but the false alarm rate of the system is higher.

Performing real-time anomaly detection on the next ICS flow sequence with different sampling frequencies, wherein the defined sampling frequency gamma of the ICS flow_sampIs 60(s). At this point the aggregation scale is modified to 60s and the time span of the original sequence is increased. In this case, SARIMA (1,0,1) x (0,0,1) is obtained₆Is an optimal model. The operation starting time of the model is 04:14:26, and the operation state of the flow online monitoring model is maintained. And intercepting a certain small-period operation model parameter.

SARIMA(1,0,1)x(0,0,1)₆The economic parameters of the model are shown in table 3 below:

TABLE 3

SARIMA(1,0,1)x(0,0,1)₆The fitting parameters of the model are shown in table 4 below:

TABLE 4

SARIMA(1,0,1)x(0,0,1)₆Model (model)The fitting situation is shown in fig. 11(a), and the prediction situation of the model is shown in fig. 11 (b). At the moment, the selected ICS flow verification sequence is an abnormal flow sequence and is obtained by the test host machine through TCP-flooding attack on the industrial switch; the abnormal flow sequence is reflected in a time sequence that the flow in unit time is suddenly increased, and the flow is recovered to be normal after the attack. The TCP-flooding attack injects at 14:36:50 with a traffic bump, and then stops injecting abnormal traffic at 14:38: 00.

The upper bound time series of the threshold interval at this time is:

the lower bound time sequence of the threshold interval at this time is:

the time stamps corresponding to the time sequence are as follows:

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

through the model, the normal flow threshold value of a certain practical moment in a short time interval can be obtained approximately. For example, at time 14:36:20, the normal flow rate is in the interval 350.3,1650.3 when the confidence interval is 95%, and in the interval 398.3,1621.6 when the confidence interval is 90%. Therefore, a threshold model of the communication flow of the ICS at a certain time under the normal condition can be obtained. The real-time of the abnormal event is as follows:

where n is 3, indicating the presence of three anomalies.

t_debug＝4，γ_samp60, available as

By the formula

It can be known that

Small-cycle number of iterations for ICS exception event occurrence

Are all 17.

By computational analysis, now present

The training sequence for the 16 th epoch remains approximately continuous in time with the 17 th epoch. Therefore, the real-time performance of the ICS communication flow online detection method is well guaranteed.

Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm

The following algorithm can be used:

(wherein k is 1,2, 3) can be represented by

The integral threshold deviation of the minor cycle iteration is measured, and the evaluation of the abnormal degree of the ICS network has the indexLeading significance. The analysis of variance shows that the normal fluctuation range of the normal flow is 0-1475.079.

From the above analysis, the number of small-cycle iterations when the ICS abnormal event occurs is 17, and the timing when the flow abnormality occurs is 14: 37. Since the flow rate is aggregated on the aggregation scale of 1 minute, the anomaly detection is judged to be anomalous, and the flow rate exceeds the upper threshold limit thereof, in the time period of 14:37 to 14: 38.

There are cases where both of the first two sequences of the ICS network traffic verification set exceed the upper threshold. The abnormal parameters of the ICS communication network in the iteration are as follows:

therefore, a condition of significant ICS flow abnormality occurs in the period of 14:36:20-14:38:20, and a condition of slight ICS flow abnormality occurs in the period of 14:38:20-14:39: 20. ICS flow abnormal degree sequence

And subsequent flow has no abnormal condition.

By the end of this 17 th iteration mini-cycle, the ICS communication network traffic analysis for the next iteration mini-cycle can be performed, by definition.

FIG. 12 is a graph of T's that differ among 15 iterations of small period analysis_trai/T_foreThe proportion has the lowest delay time ratio of 100/30 and the lowest delay time ratio of 500/30 is larger than 100/30, but the accuracy of the algorithm is higher than 100/30. The difference of the BIC values under different training scales is small, and the algorithm is proved to have better convergence.

The ICS flow modes reflected by different network anomalies are different, abnormal conditions of the flow model can be analyzed according to different side points in subsequent anomaly detection, and a machine learning algorithm is introduced to classify abnormal flows.

Claims (3)

1. An industrial control system communication flow online monitoring method based on SARIMA is characterized by comprising the following steps:

1) deploying an industrial switch, a monitoring host and a testing host in an ICS communication network, wherein the monitoring host collects communication flow data from the industrial switch in real time;

2) establishing a short period SARIMA (P, D, Q) x (P, D, Q) according to the selected flow aggregation scale_sA model; the method specifically comprises the following steps:

2.1) defining selected flow aggregation scales and small period analysis scales using SARIMA (P, D, Q) x (P, D, Q)_sSequence definition to generate SARIMA (P, D, Q) x (P, D, Q)_sTime series of (a):

SARIMA(p,d,q)x(P,D,Q)_sthe model is obtained by respectively carrying out D-order difference calculation and D-order seasonal difference calculation on an ARMA (p, q) model, wherein the ARMA (p, q) model is formed by combining an AR (p) model and an MA (q) model;

the autoregressive moving average model ARMA (p, q) is defined as follows:

X_t＝φ₁X_t-1+φ₂X_t-2+…+φ_pX_t-p+ε_t-θ₁ε_t-1-…-θ_qε_t-q

the formula is as follows: x_tA stationary time sequence of small periods after the averaging process; phi is a_pIs the coefficient of the autoregressive term AR; theta_qCoefficient of the moving average term MA; epsilon_tIs random disturbance; p is ARThe order; q is the order of MA;

defining a delay operator B, BX_t＝X_t-1Then the polynomial of the AR coefficient phi_p(B)＝1-φ₁B-…-φ_p(B)^pCoefficient of MA polynomial theta_q(B)＝1-θ₁B-…-θ_q(B)^q；

Introduction of a difference operator Δ^d＝(1-B)^dThen the ARIMA (p, d, q) model is expressed as:

Φ_p(B)Δ^dX_t＝Θ_q(B)ε_t

the SARIMA model is obtained by carrying out seasonal difference operation on the ARIMA model, and is defined as follows:

Φ_p(B)Φ_P(B^s)Δ^dΔ_s ^DX_t＝Θ_q(B)Θ_Q(B^s)ε_t

wherein epsilon_tIs a white noise sequence, D is the order of trend difference, D is the seasonal difference order with period s as compensation, B^sFor delay operators of order s, Δ_s ^DIs a seasonal difference operator; b is^sX_t＝X_t-s，Δ_s ^D＝1-B^s，Θ_Q(B^s) Is B^sPolynomial of order Q,. phi_P(B^s) Is B^sA polynomial of order P;

2.2) pairs of SARIMA (P, D, Q) x (P, D, Q) using Bayesian information criterion BIC_sCarrying out supervision analysis and order determination on the P, D, Q, P, D and Q orders of the model;

2.3) applying least squares method to SARIMA (P, D, Q) x (P, D, Q)_sP-order coefficient of (phi)_k(k 1,2, …, p), q-order coefficient θ_k(k-1, 2, …, q), and seasonal P-th order coefficient

Seasonal Q-order coefficient

Carrying out estimation;

2.4) SARIMA (P, D, Q) x (P, D, Q) using the optimal BIC criterion_sThe model performs fitting analysis on the original sequence and performs residual error detection; if the residual error is white noise, performing inverse filtering processing on the fitting sequence to obtain a fitting value or a predicted value of the original sequence; if the residual error is not white noise, the BIC information criterion is adopted again to carry out order fixing on the ARMA (p, q) model;

2.5) obtaining a small period of SARIMA (P, D, Q) x (P, D, Q)_sThe mathematical expression of (a);

3) generating an ICS dynamic flow threshold model according to different confidence degrees;

4) and dynamically analyzing the flow sequence acquired in real time.

2. The online monitoring method for communication flow of SARIMA-based industrial control system as claimed in claim 1, wherein the step 3) is specifically as follows:

3.1) collecting real-time flow data on an ICS industrial exchanger according to a set sampling frequency gamma_sampGenerating a time sequence by a polymerization scale, and taking a small period as an iteration period;

3.2) training the collected real-time flow data, and carrying out SARIMA (P, D, Q) x (P, D, Q) in a small period_sAfter the model is trained and adapted, outputting an optimal model and parameters for model adaptation; the model for the ith small period is defined as:

wherein f is_SARIMA() Is SARIMA (P, D, Q) x (P, D, Q)_sThe functional expression of the model is shown as,

small period training set, T, for the ith iteration_traiFor a small period of training, T_foreThe number of sequences predicted for a small period, s is a periodic parameter, 'BIC' is a criterion for selecting optimal (P, D, Q, P, D, Q) parameters for the economics of metrology,

is a time sequence predicted by the ith iteration of the SARIMA model;

calculating the predicted mean value of the ith iteration

3.3) distributed operation of short-period SARIMA (P, D, Q) x (P, D, Q)_sThe model is used for carrying out real-time dynamic rolling modeling on ICS flow data of the ith small period collected in real time, wherein the real-time flow data at the moment is equivalent to a verification set; comparing and analyzing the verification set with the upper and lower flow threshold limits based on the confidence interval generated in the prediction process; the upper and lower threshold values of the defined ith small period are as follows:

wherein

Is at 1-alpha_P.IFor the value of the z-distribution of significance,

is the upper bound of the dynamic threshold interval for the ith small period,

is the lower bound of the dynamic threshold interval,

are all of length T_foreTime series of (a)_P.IIs the confidence level;

the normal ICS traffic for the ith iteration is defined as:

wherein

ICS communication traffic, T, collected in real time for the ith small cycle_foreIs the amount of sample collected;

ICS prediction sequence of i-th minor cycle

And training sequence

There is a real-time correlation, the predicted sequence of the ith minor period

Can be estimated as:

wherein the function n is the intersection of the two time series sets;

3.4) after the flow judgment is finished, continuing training iteration of the next small period, outputting a new optimal model and model adaptive parameters again, and judging newly input real-time flow data again;

3.5) the whole process is circulated until the set iteration number is reached.

3. The SARIMA-based industrial control system communication flow online monitoring method as claimed in claim 2, wherein the abnormal sequence of the ith iteration

The time series judgment calculation method of (2) is as follows:

wherein each iteration of the small-period sequence to be detected needs to pass through distributed operation SARIMA (P, D, Q) x (P, D, Q)_sA model for generating an anomaly matrix when the small-period sequence to be detected triggers an anomaly detection condition

When the value of the flow rate is in the normal interval range

Wherein

Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence

Wherein n is the total number of occurrences of the anomaly; the sequence of the number of samples of the small period sequence number when the abnormality occurs

Comprises the following steps:

wherein t is_debugFor real-time debugging before program runningM, y_sampIs the sampling frequency of the ICS flow;

the nth ICS abnormality occurs at

The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs

The following equation system can be used to obtain the calculation:

wherein K_nIs an intermediate variable;

variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm

The following algorithm can be used:

wherein

Can be aligned with

The integral threshold deviation of the minor-cycle iteration is measured, and the method has guiding significance for analyzing the abnormal degree of the ICS network;

integrally analyzing abnormal degree of ICS network flow of ith iteration

The calculation method comprises the following steps:

wherein when an element is present

A slight ICS traffic anomaly in time,

for significant ICS traffic anomalies, when

The ICS flow is abnormal to the highest degree, and the ICS network is affected destructively at the moment.

Images