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

Industrial control system communication flow online monitoring method based on SARIMA Download PDF

Info

Publication number
CN112187528A
CN112187528A CN202010969136.9A CN202010969136A CN112187528A CN 112187528 A CN112187528 A CN 112187528A CN 202010969136 A CN202010969136 A CN 202010969136A CN 112187528 A CN112187528 A CN 112187528A
Authority
CN
China
Prior art keywords
ics
model
sarima
flow
sequence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010969136.9A
Other languages
Chinese (zh)
Other versions
CN112187528B (en
Inventor
杨强
郝唯杰
杨涛
阮伟
王文海
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202010969136.9A priority Critical patent/CN112187528B/en
Publication of CN112187528A publication Critical patent/CN112187528A/en
Application granted granted Critical
Publication of CN112187528B publication Critical patent/CN112187528B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L41/00—Arrangements for maintenance, administration or management of data switching networks, e.g. of packet switching networks
    • H04L41/08—Configuration management of networks or network elements
    • H04L41/0803—Configuration setting
    • H04L41/0823—Configuration setting characterised by the purposes of a change of settings, e.g. optimising configuration for enhancing reliability
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L41/00—Arrangements for maintenance, administration or management of data switching networks, e.g. of packet switching networks
    • H04L41/14—Network analysis or design
    • H04L41/145—Network analysis or design involving simulating, designing, planning or modelling of a network
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L43/00—Arrangements for monitoring or testing data switching networks
    • H04L43/08—Monitoring or testing based on specific metrics, e.g. QoS, energy consumption or environmental parameters
    • H04L43/0876—Network utilisation, e.g. volume of load or congestion level
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L63/00—Network architectures or network communication protocols for network security
    • H04L63/14—Network architectures or network communication protocols for network security for detecting or protecting against malicious traffic
    • H04L63/1408—Network architectures or network communication protocols for network security for detecting or protecting against malicious traffic by monitoring network traffic

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 networksModeling 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 verifiedThe reliability and accuracy of the algorithm.

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 gammasampGenerating 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 periodsAfter 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:
Figure RE-GDA0002764332400000021
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure RE-GDA0002764332400000022
small period training set, T, for the ith iterationforeThe 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,
Figure RE-GDA0002764332400000023
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure RE-GDA0002764332400000031
Figure RE-GDA0002764332400000032
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; it defines the i-th sub-periodThe bound threshold is:
Figure RE-GDA0002764332400000033
Figure RE-GDA0002764332400000034
wherein
Figure RE-GDA0002764332400000035
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure RE-GDA0002764332400000036
is the upper bound of the dynamic threshold interval for the ith small period,
Figure RE-GDA0002764332400000037
is the lower bound of the dynamic threshold interval,
Figure RE-GDA0002764332400000038
are all of length TforeTime series of (a)P.IIs the confidence level;
the communication flow collected in real time and the extremely low algorithm complexity ensure
Figure RE-GDA0002764332400000039
And
Figure RE-GDA00027643324000000310
the temporal coincidence, the normal ICS communication traffic (each time-sequential sample within the time series) for the ith iteration is defined as:
Figure RE-GDA00027643324000000311
the algorithm is essentially a time series comparison of dynamically acquired flows
Figure RE-GDA00027643324000000312
Two threshold time series with same length
Figure RE-GDA00027643324000000313
And
Figure RE-GDA00027643324000000314
the magnitude relation of each corresponding time sequence sample; wherein
Figure RE-GDA00027643324000000315
ICS communication traffic, T, collected in real time for the ith small cycleforeIs the amount of sample collected;
ICS prediction sequence of i-th minor cycle
Figure RE-GDA00027643324000000316
And training sequence
Figure RE-GDA00027643324000000317
There is a real-time correlation, the predicted sequence of the ith minor period
Figure RE-GDA00027643324000000318
Can be estimated as:
Figure RE-GDA00027643324000000319
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
Figure RE-GDA0002764332400000041
Time series ofThe judgment calculation method is as follows (the calculation mode of the default time sequence is to calculate each time sequence sample corresponding to one in the time sequence):
Figure RE-GDA0002764332400000042
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
Figure RE-GDA0002764332400000043
When the value of the flow rate is in the normal interval range
Figure RE-GDA0002764332400000044
Wherein
Figure RE-GDA0002764332400000045
Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence
Figure RE-GDA0002764332400000046
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
Figure RE-GDA0002764332400000047
Comprises the following steps:
Figure RE-GDA0002764332400000048
wherein t isdebugFor real-time debugging time before program run, gammasampIs the sampling frequency of the ICS traffic.
The nth ICS abnormality occurs at
Figure RE-GDA0002764332400000049
In a sub-small cycle iterative calculation, therefore ICNumber of small-cycle iterations of S-anomaly occurrence
Figure RE-GDA00027643324000000410
(
Figure RE-GDA00027643324000000411
Can be calculated retrospectively from the following system of equations:
Figure RE-GDA00027643324000000412
wherein KnIs an intermediate variable.
Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure RE-GDA00027643324000000413
The following algorithm can be used:
Figure RE-GDA0002764332400000051
where k is 1,2, …, n
Figure RE-GDA0002764332400000052
Figure RE-GDA0002764332400000053
Can be aligned with
Figure RE-GDA0002764332400000054
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
Figure RE-GDA0002764332400000055
The calculation method comprises the following steps:
Figure RE-GDA0002764332400000056
wherein when an element is present
Figure RE-GDA0002764332400000057
A slight ICS traffic anomaly in time,
Figure RE-GDA0002764332400000058
for significant ICS traffic anomalies, when
Figure RE-GDA0002764332400000059
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)10Fitting situation graph of the model;
FIG. 9(b) is SARIMA (2,1,3) x (2,0,0)10A 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)6The model fitting situation is as shown in the figure;
FIG. 11(b) is SARIMA (1,0,1) x (0,0,1)6A prediction effect and abnormality detection map of the model;
FIG. 12 is a variation of Ttrai/TforeThe 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 the original sequence needs to be subjected to hysteresis lag term exceeding a critical valuePeriodic decomposition, i.e. requiring SARIMA (P, D, Q) x (P, D, Q) of the sequencesAnd (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)sDefinition of sequence toProduction of 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:
Xt＝φ1Xt-1+φ2Xt-2+…+φpXt-p+t-θ1 t-1-…-θq t-q
the formula is as follows: xtThe 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 apIs the coefficient of the autoregressive term AR; thetaqCoefficient of the moving average term MA;tis random disturbance; p is the order of AR; q is the order of MA;
defining a delay operator B, BXt＝Xt-1If the polynomial of the AR coefficient phi (B) is 1-phi1B-…-φp(B)pThe MA coefficient polynomial Θ (B) is 1- θ1B-…-θq(B)q；
Introduction of a difference operator Δd＝(1-B)dThen the ARIMA (p, d, q) model is expressed as:
Φ(B)ΔdXt＝Θ(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(Bs)ΔdΔs DXt＝Θq(B)ΘQ(Bs)t
whereintIs a white noise sequence, D is the order of trend difference, D is the seasonal difference order with period s as compensation, BsFor delay operators of order s, Δs DIs a seasonal difference operator; b issXt＝Xt-s，Δs D＝1-Bs，ΦP(Bs) Is composed ofBsPolynomial of order Q,. phiP(Bs) Is BsA polynomial of order P;
2.2) pairs of SARIMA (P, D, Q) x (P, D, Q) using Bayesian information criterion BICsAnd (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
Figure RE-GDA0002764332400000096
Seasonal Q-order coefficient
Figure RE-GDA0002764332400000097
Carrying out estimation;
2.4) SARIMA (P, D, Q) x (P, D, Q) using the optimal BIC criterionsThe 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 gammasampGenerating a time sequence by a polymerization scale, and taking a small period as an iteration period;
3.2) training the collected real-time traffic data, as shown in FIG. 2, atSARIMA (P, D, Q) x (P, D, Q) in a small cyclesAfter 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:
Figure RE-GDA0002764332400000091
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure RE-GDA0002764332400000092
small period training set, T, for the ith iterationforeThe 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,
Figure RE-GDA0002764332400000093
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure RE-GDA0002764332400000094
Figure RE-GDA0002764332400000095
3.3) distributed operation of SARIMA (P, D, Q) x (P, D, Q) for a short period as shown in FIG. 3sThe 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:
Figure RE-GDA0002764332400000101
Figure RE-GDA0002764332400000102
wherein
Figure RE-GDA0002764332400000103
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure RE-GDA0002764332400000104
is the upper bound of the dynamic threshold interval for the ith small period,
Figure RE-GDA0002764332400000105
is the lower bound of the dynamic threshold interval,
Figure RE-GDA0002764332400000106
are all of length TforeTime series of (a)P.IIs the confidence level;
the normal ICS traffic for the ith iteration is defined as:
Figure RE-GDA0002764332400000107
wherein
Figure RE-GDA0002764332400000108
ICS communication traffic, T, collected in real time for the ith small cycleforeIs the amount of sample collected;
ICS prediction sequence of i-th minor cycle
Figure RE-GDA0002764332400000109
And training sequence
Figure RE-GDA00027643324000001010
There is a real-time correlation, the predicted sequence of the ith minor period
Figure RE-GDA00027643324000001011
Can be estimated as:
Figure RE-GDA00027643324000001012
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
Figure RE-GDA00027643324000001013
The time series judgment calculation method of (2) is as follows:
Figure RE-GDA00027643324000001014
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
Figure RE-GDA0002764332400000111
When the value of the flow rate is in the normal interval range
Figure RE-GDA0002764332400000112
Wherein
Figure RE-GDA0002764332400000113
Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence
Figure RE-GDA0002764332400000114
Where n is the total number of occurrences of the anomaly. Samples of the small cycle sequence number when an anomaly occursNumber sequence
Figure RE-GDA0002764332400000115
Comprises the following steps:
Figure RE-GDA0002764332400000116
wherein t isdebugFor real-time debugging time before program run, gammasampIs the sampling frequency of the ICS traffic.
The nth ICS abnormality occurs at
Figure RE-GDA0002764332400000117
The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs
Figure RE-GDA0002764332400000118
(
Figure RE-GDA0002764332400000119
Can be calculated retrospectively from the following system of equations:
Figure RE-GDA00027643324000001110
wherein KnIs an intermediate variable.
Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure RE-GDA00027643324000001111
The following algorithm can be used:
Figure RE-GDA00027643324000001112
where k is 1,2, …, n
Figure RE-GDA00027643324000001113
Figure RE-GDA00027643324000001114
Can be aligned with
Figure RE-GDA00027643324000001115
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
Figure RE-GDA00027643324000001116
The calculation method comprises the following steps:
Figure RE-GDA00027643324000001117
wherein when an element is present
Figure RE-GDA00027643324000001118
A slight ICS traffic anomaly in time,
Figure RE-GDA00027643324000001119
for significant ICS traffic anomalies, when
Figure RE-GDA00027643324000001120
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, and dynamically models and analyzes to generateAnd carrying out real-time identification and early warning on the acquired ICS network flow in the flow threshold interval. Sampling frequency gamma of defined ICS flowsampIs 1(s).
SARIMA (P, D, Q) x (P, D, Q) for ICS communication network trafficsAnd (6) modeling. Defining the small period training item number as Ttrai＝300，Tfore30 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 aggregated10. The parameters of the economics of the measurement are shown in the following table:
TABLE 1
Figure RE-GDA0002764332400000121
Its SARIMA (2,1,3) x (2,0,0)10The table of the parameter fit of the model is as follows:
TABLE 2
Figure RE-GDA0002764332400000122
Figure RE-GDA0002764332400000131
SARIMA(2,1,3)x(2,0,0)10The 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
Figure RE-GDA0002764332400000132
R-Square ═ 0.89. SARIMA (2,1,3) x (2,0,0) is thus obtainable10For the small period TtraiThe ICS traffic of 300 has strong explanatory property, and the threshold interval represents Tfore30 ICS communication traffic. From the threshold interval and experiment in FIG. 9(b)The verification sequence can know that the ICS communication flow sequence on the time sequence accords with the flow characteristics and rules of the training set, so that the upper and lower boundaries of the threshold interval are not exceeded. The upper bound time series of the threshold intervals is:
Figure RE-GDA0002764332400000133
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)10After modeling, it verifies the period TforeThe 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
Figure RE-GDA0002764332400000134
As an empty set, the sequence of the number of iterations of the sequence in a small period when an abnormal event occurs
Figure RE-GDA0002764332400000135
Is an empty set. Variance at this time
Figure RE-GDA0002764332400000136
Is 0. Sequence of degree of abnormality thereof
Figure RE-GDA0002764332400000137
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 flowsampIs 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 obtained6Is 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)6The economic parameters of the model are shown in table 3 below:
TABLE 3
Figure RE-GDA0002764332400000141
SARIMA(1,0,1)x(0,0,1)6The fitting parameters of the model are shown in table 4 below:
TABLE 4
Figure RE-GDA0002764332400000142
SARIMA(1,0,1)x(0,0,1)6The model fitting situation is shown in fig. 11(a), and the model prediction situation 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:
Figure RE-GDA0002764332400000151
the lower bound time sequence of the threshold interval at this time is:
Figure RE-GDA0002764332400000152
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:
Figure RE-GDA0002764332400000153
where n is 3, indicating the presence of three anomalies.
tdebug＝4，γsamp60, available as
Figure RE-GDA0002764332400000154
Figure RE-GDA0002764332400000155
By the formula
Figure RE-GDA0002764332400000156
It can be known that
Figure RE-GDA0002764332400000157
Small-cycle number of iterations for ICS exception event occurrence
Figure RE-GDA0002764332400000158
Are all 17.
By computational analysis, now present
Figure RE-GDA0002764332400000159
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
Figure RE-GDA0002764332400000161
The following algorithm can be used:
Figure RE-GDA0002764332400000162
Figure RE-GDA0002764332400000163
(wherein k is 1,2, 3) can be represented by
Figure RE-GDA0002764332400000164
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. 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.
Figure RE-GDA0002764332400000165
Figure RE-GDA0002764332400000166
Figure RE-GDA0002764332400000167
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:
Figure RE-GDA0002764332400000168
Figure RE-GDA0002764332400000169
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
Figure RE-GDA00027643324000001610
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 analysistrai/TforeThe 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 (4)

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 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.
2. The online monitoring method for communication flow of SARIMA-based industrial control system as claimed in claim 1, wherein the step 2) is specifically as follows:
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:
Xt＝φ1Xt-1+φ2Xt-2+…+φpXt-p+t-θ1 t-1-…-θq t-q
the formula is as follows: xtThe 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 apIs the coefficient of the autoregressive term AR; thetaqCoefficient of the moving average term MA;tis random disturbance; p is the order of AR; q is the order of MA;
defining a delay operator B, BXt＝Xt-1If the polynomial of the AR coefficient phi (B) is 1-phi1B-…-φp(B)pThe MA coefficient polynomial Θ (B) is 1- θ1B-…-θq(B)q；
Introduction of a difference operator Δd＝(1-B)dThen the ARIMA (p, d, q) model is expressed as:
Φ(B)ΔdXt＝Θ(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(Bs)ΔdΔs DXt＝Θq(B)ΘQ(Bs)t
whereintIs a white noise sequence, D is the order of trend difference, D is the seasonal difference order with period s as compensation, BsFor delay operators of order s, Δs DIs a seasonal difference operator; b issXt＝Xt-s，Δs D＝1-Bs，ΦP(Bs) Is BsPolynomial of order Q,. phiP(Bs) Is BsA polynomial of order P;
2.2) pairs of SARIMA (P, D, Q) x (P, D, Q) using Bayesian information criterion BICsCarrying 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
Figure FDA0002683417000000021
Seasonal Q-order coefficient
Figure FDA0002683417000000027
Carrying out estimation;
2.4) SARIMA (P, D, Q) x (P, D, Q) using the optimal BIC criterionsThe 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. 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 gammasampGenerating 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 periodsAfter 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:
Figure FDA0002683417000000022
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure FDA0002683417000000023
small period training set, T, for the ith iterationforeThe 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,
Figure FDA0002683417000000024
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure FDA0002683417000000025
Figure FDA0002683417000000026
3.3) distributed fortuneSARIMA (P, D, Q) x (P, D, Q) with small periodssThe 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:
Figure FDA0002683417000000031
Figure FDA0002683417000000032
wherein
Figure FDA0002683417000000033
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure FDA0002683417000000034
is the upper bound of the dynamic threshold interval for the ith small period,
Figure FDA0002683417000000035
is the lower bound of the dynamic threshold interval,
Figure FDA0002683417000000036
are all of length TforeTime series of (a)P.IIs the confidence level;
the normal ICS traffic for the ith iteration is defined as:
Figure FDA0002683417000000037
wherein
Figure FDA0002683417000000038
ICS communication traffic, T, collected in real time for the ith small cycleforeIs the amount of sample collected;
ICS prediction sequence of i-th minor cycle
Figure FDA0002683417000000039
And training sequence
Figure FDA00026834170000000310
There is a real-time correlation, the predicted sequence of the ith minor period
Figure FDA00026834170000000311
Can be estimated as:
Figure FDA00026834170000000312
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.
4. The SARIMA-based industrial control system communication flow online monitoring method as claimed in claim 3, wherein the abnormal sequence of the ith iteration
Figure FDA00026834170000000313
The time series judgment calculation method of (2) is as follows:
Figure FDA0002683417000000041
in which the small-period sequence to be detected for each iteration needs to be run in a distributed mannerSARIMA(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
Figure FDA0002683417000000042
When the value of the flow rate is in the normal interval range
Figure FDA0002683417000000043
Wherein
Figure FDA0002683417000000044
Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence
Figure FDA0002683417000000045
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
Figure FDA0002683417000000046
Comprises the following steps:
Figure FDA0002683417000000047
wherein t isdebugFor real-time debugging time before program run, gammasampIs the sampling frequency of the ICS traffic.
The nth ICS abnormality occurs at
Figure FDA0002683417000000048
The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs
Figure FDA0002683417000000049
(
Figure FDA00026834170000000410
Of) can be traced back by the following system of equationsAnd calculating to obtain:
Figure FDA00026834170000000411
wherein KnIs an intermediate variable.
Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure FDA00026834170000000417
The following algorithm can be used:
Figure FDA00026834170000000412
wherein
Figure FDA00026834170000000413
Figure FDA00026834170000000414
Can be aligned with
Figure FDA00026834170000000415
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.
Integrally analyzing abnormal degree of ICS network flow of ith iteration
Figure FDA00026834170000000416
The calculation method comprises the following steps:
Figure FDA0002683417000000051
wherein when an element is present
Figure FDA0002683417000000052
Slight ICS flow at the timeIn the event of an abnormality,
Figure FDA0002683417000000053
for significant ICS traffic anomalies, when
Figure FDA0002683417000000054
The ICS flow is abnormal to the highest degree, and the ICS network is affected destructively at the moment.
CN202010969136.9A 2020-09-15 2020-09-15 Industrial control system communication flow online monitoring method based on SARIMA Active CN112187528B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010969136.9A CN112187528B (en) 2020-09-15 2020-09-15 Industrial control system communication flow online monitoring method based on SARIMA

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010969136.9A CN112187528B (en) 2020-09-15 2020-09-15 Industrial control system communication flow online monitoring method based on SARIMA

Publications (2)

Publication Number Publication Date
CN112187528A true CN112187528A (en) 2021-01-05
CN112187528B CN112187528B (en) 2021-10-08

Family

ID=73921186

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010969136.9A Active CN112187528B (en) 2020-09-15 2020-09-15 Industrial control system communication flow online monitoring method based on SARIMA

Country Status (1)

Country Link
CN (1) CN112187528B (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113411216A (en) * 2021-06-21 2021-09-17 国网宁夏电力有限公司信息通信公司 Network flow prediction method based on discrete wavelet transform and FA-ELM
CN113534731A (en) * 2021-07-16 2021-10-22 珠海市鸿瑞信息技术股份有限公司 Download data security analysis system and method based on industrial control
CN114154333A (en) * 2021-12-06 2022-03-08 西安电子科技大学 Atmospheric temperature profile prediction method
CN114499934A (en) * 2021-12-16 2022-05-13 西安交通大学 Intrusion detection method and system based on fusion learning in industrial Internet of things
CN115931055A (en) * 2023-01-06 2023-04-07 长江信达软件技术（武汉）有限责任公司 Rural water supply operation diagnosis method and system based on big data analysis

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160219067A1 (en) * 2015-01-28 2016-07-28 Korea Internet & Security Agency Method of detecting anomalies suspected of attack, based on time series statistics
US20170329314A1 (en) * 2014-11-26 2017-11-16 Shenyang Institute Of Automation, Chinese Academy Of Sciences Modbus tcp communication behaviour anomaly detection method based on ocsvm dual-outline model
CN107547269A (en) * 2017-08-14 2018-01-05 浙江大学 The construction method of intelligent substation communication flow threshold model based on FARIMA
CN107733905A (en) * 2017-10-24 2018-02-23 北京威努特技术有限公司 A kind of detection method of industry control network unit exception flow

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170329314A1 (en) * 2014-11-26 2017-11-16 Shenyang Institute Of Automation, Chinese Academy Of Sciences Modbus tcp communication behaviour anomaly detection method based on ocsvm dual-outline model
US20160219067A1 (en) * 2015-01-28 2016-07-28 Korea Internet & Security Agency Method of detecting anomalies suspected of attack, based on time series statistics
CN107547269A (en) * 2017-08-14 2018-01-05 浙江大学 The construction method of intelligent substation communication flow threshold model based on FARIMA
CN107733905A (en) * 2017-10-24 2018-02-23 北京威努特技术有限公司 A kind of detection method of industry control network unit exception flow

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
MACIEJ SZMIT等: "sage of Pseudo-estimator LAD and SARIMA Models for Network Traffic Prediction", 《SPRINGER BERLIN HEIDELBERG》 *
TAO YANG等: "Fuzzy information granulation and ED-LSTM based traffic prediction of industrial control systems", 《 2020 CHINESE CONTROL AND DECISION CONFERENCE (CCDC》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113411216A (en) * 2021-06-21 2021-09-17 国网宁夏电力有限公司信息通信公司 Network flow prediction method based on discrete wavelet transform and FA-ELM
CN113411216B (en) * 2021-06-21 2022-11-04 国网宁夏电力有限公司信息通信公司 Network flow prediction method based on discrete wavelet transform and FA-ELM
CN113534731A (en) * 2021-07-16 2021-10-22 珠海市鸿瑞信息技术股份有限公司 Download data security analysis system and method based on industrial control
CN114154333A (en) * 2021-12-06 2022-03-08 西安电子科技大学 Atmospheric temperature profile prediction method
CN114499934A (en) * 2021-12-16 2022-05-13 西安交通大学 Intrusion detection method and system based on fusion learning in industrial Internet of things
CN114499934B (en) * 2021-12-16 2022-12-09 西安交通大学 Intrusion detection method and system based on fusion learning in industrial Internet of things
CN115931055A (en) * 2023-01-06 2023-04-07 长江信达软件技术（武汉）有限责任公司 Rural water supply operation diagnosis method and system based on big data analysis
CN115931055B (en) * 2023-01-06 2023-06-16 长江信达软件技术（武汉）有限责任公司 Rural water supply operation diagnosis method and system based on big data analysis

Also Published As

Publication number Publication date
CN112187528B (en) 2021-10-08

Similar Documents

Publication Publication Date Title
CN112202736B (en) Communication network anomaly classification method based on statistical learning and deep learning
CN112187528B (en) Industrial control system communication flow online monitoring method based on SARIMA
Coble et al. Applying the general path model to estimation of remaining useful life
CN109739214B (en) Method for detecting intermittent faults in industrial process
CN111259717B (en) Method and system for judging abnormal state of rotating equipment
CN112799898B (en) Interconnection system fault node positioning method and system based on distributed fault detection
CN113239132B (en) Online out-of-tolerance identification method for voltage transformer
CN112149877B (en) Multi-source data driven fault prediction method and system for multi-element complex urban power grid
Wang et al. An accurate false data detection in smart grid based on residual recurrent neural network and adaptive threshold
Li et al. Deep learning based covert attack identification for industrial control systems
Fujita et al. An approach for intelligent evaluation of the state of complex autonomous objects based on the wavelet analysis
Horvath et al. Sensor fault diagnosis of inland navigation system using physical model and pattern recognition approach
CN113269327A (en) Flow anomaly prediction method based on machine learning
CN116796261B (en) Closed switch equipment mechanical characteristic prediction method based on artificial intelligence
Li et al. Real-time detecting false data injection attacks based on spatial and temporal correlations
CN111934903A (en) Docker container fault intelligent prediction method based on time sequence evolution genes
Akbarian et al. Attack resilient cloud-based control systems for industry 4.0
Zhang et al. Set-membership filtering approach for fault detection of systems with unknown-but-bounded noises
CN115685045A (en) Online evaluation method for voltage transformer
CN111966966B (en) Method and system for analyzing feasible domain of sensor measurement error model parameters
Cao et al. An adaptive UKF algorithm for process fault prognostics
CN112085043A (en) Intelligent monitoring method and system for network security of transformer substation
Cross et al. Model-based condition monitoring for wind turbines
Aftabi et al. A variational autoencoder framework for robust, physics-informed cyberattack recognition in industrial cyber-physical systems
Athalye et al. Model-based cps attack detection techniques: Strengths and limitations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant