CN112187528B - 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
CN112187528B
CN112187528B CN202010969136.9A CN202010969136A CN112187528B CN 112187528 B CN112187528 B CN 112187528B CN 202010969136 A CN202010969136 A CN 202010969136A CN 112187528 B CN112187528 B CN 112187528B
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.)
Active
Application number
CN202010969136.9A
Other languages
Chinese (zh)
Other versions
CN112187528A (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

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Computer Security & Cryptography (AREA)
  • Environmental & Geological Engineering (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Exchanges In Wide-Area Networks (AREA)

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 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 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 GDA0002764332400000021
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure GDA0002764332400000022
small period training set, T, for the ith iterationforeThe 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,
Figure GDA0002764332400000023
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure GDA0002764332400000031
Figure 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; the upper and lower threshold values of the defined ith small period are as follows:
Figure GDA0002764332400000033
Figure GDA0002764332400000034
wherein
Figure GDA0002764332400000035
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure GDA0002764332400000036
is the upper bound of the dynamic threshold interval for the ith small period,
Figure GDA0002764332400000037
is the lower bound of the dynamic threshold interval,
Figure 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 GDA0002764332400000039
And
Figure 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 GDA00027643324000000311
the algorithm is essentially a time series comparison of dynamically acquired flows
Figure GDA00027643324000000312
Two threshold time series with same length
Figure GDA00027643324000000313
And
Figure GDA00027643324000000314
the magnitude relation of each corresponding time sequence sample; wherein
Figure 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 GDA00027643324000000316
And training sequence
Figure GDA00027643324000000317
With real-time correlation, i-th minor period of predictionSequencing
Figure GDA00027643324000000318
Can be estimated as:
Figure 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 GDA0002764332400000041
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):
Figure 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 GDA0002764332400000043
When the value of the flow rate is in the normal interval range
Figure GDA0002764332400000044
Wherein
Figure GDA0002764332400000045
Assume that an ICS anomaly has occurredReal time to time as a sequence
Figure 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 GDA0002764332400000047
Comprises the following steps:
Figure 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 GDA0002764332400000049
The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs
Figure GDA00027643324000000410
(
Figure GDA00027643324000000411
Can be calculated retrospectively from the following system of equations:
Figure GDA00027643324000000412
wherein KnIs an intermediate variable.
Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure GDA00027643324000000413
The following algorithm can be used:
Figure GDA0002764332400000051
where k is 1,2, …, n
Figure GDA0002764332400000052
Figure GDA0002764332400000053
Can be aligned with
Figure 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 GDA0002764332400000055
The calculation method comprises the following steps:
Figure GDA0002764332400000056
wherein when an element is present
Figure GDA0002764332400000057
A slight ICS traffic anomaly in time,
Figure GDA0002764332400000058
for significant ICS traffic anomalies, when
Figure 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 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:
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; epsilontIs 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
wherein epsilontIs 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 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 GDA0002764332400000096
Seasonal Q-order coefficient
Figure 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 flow data, as shown in FIG. 2, SARIMA (P, D, Q) x (P, D, Q) is performed 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 GDA0002764332400000091
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure 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 GDA0002764332400000093
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure GDA0002764332400000094
Figure 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 GDA0002764332400000101
Figure GDA0002764332400000102
wherein
Figure GDA0002764332400000103
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure GDA0002764332400000104
is the upper bound of the dynamic threshold interval for the ith small period,
Figure GDA0002764332400000105
is the lower bound of the dynamic threshold interval,
Figure 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 GDA0002764332400000107
wherein
Figure 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 GDA0002764332400000109
And training sequence
Figure GDA00027643324000001010
There is a real-time correlation, the predicted sequence of the ith minor period
Figure GDA00027643324000001011
Can be estimated as:
Figure 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 GDA00027643324000001013
The time series judgment calculation method of (2) is as follows:
Figure GDA00027643324000001014
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
Figure GDA0002764332400000111
When the value of the flow rate is in the normal interval range
Figure GDA0002764332400000112
Wherein
Figure GDA0002764332400000113
Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence
Figure GDA0002764332400000114
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 GDA0002764332400000115
Comprises the following steps:
Figure 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 GDA0002764332400000117
The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs
Figure GDA0002764332400000118
(
Figure GDA0002764332400000119
Can be calculated retrospectively from the following system of equations:
Figure GDA00027643324000001110
wherein KnIs an intermediate variable.
Variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure GDA00027643324000001111
The following algorithm can be used:
Figure GDA00027643324000001112
where k is 1,2, …, n
Figure GDA00027643324000001113
Figure GDA00027643324000001114
Can be aligned with
Figure 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 GDA00027643324000001116
The calculation method comprises the following steps:
Figure GDA00027643324000001117
wherein when an element is present
Figure GDA00027643324000001118
A slight ICS traffic anomaly in time,
Figure GDA00027643324000001119
for significant ICS traffic anomalies, when
Figure 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, 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 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 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 GDA0002764332400000122
Figure 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 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. 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:
Figure 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 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 GDA0002764332400000135
Is an empty set. Variance at this time
Figure GDA0002764332400000136
Is 0. Sequence of degree of abnormality thereof
Figure 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 GDA0002764332400000141
SARIMA(1,0,1)x(0,0,1)6The fitting parameters of the model are shown in table 4 below:
TABLE 4
Figure GDA0002764332400000142
SARIMA(1,0,1)x(0,0,1)6Model (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:
Figure GDA0002764332400000151
the lower bound time sequence of the threshold interval at this time is:
Figure 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 GDA0002764332400000153
where n is 3, indicating the presence of three anomalies.
tdebug＝4，γsamp60, available as
Figure GDA0002764332400000154
Figure GDA0002764332400000155
By the formula
Figure GDA0002764332400000156
It can be known that
Figure GDA0002764332400000157
Small-cycle number of iterations for ICS exception event occurrence
Figure GDA0002764332400000158
Are all 17.
By computational analysis, now present
Figure 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 GDA0002764332400000161
The following algorithm can be used:
Figure GDA0002764332400000162
Figure GDA0002764332400000163
(wherein k is 1,2, 3) can be represented by
Figure GDA0002764332400000164
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.
Figure GDA0002764332400000165
Figure GDA0002764332400000166
Figure 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 GDA0002764332400000168
Figure 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 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 (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 scalesA 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:
Xt＝φ1Xt-1+φ2Xt-2+…+φpXt-p+εt-θ1εt-1-…-θqεt-q
the formula is as follows: xtA stationary time sequence of small periods after the averaging process; phi is apIs the coefficient of the autoregressive term AR; thetaqCoefficient of the moving average term MA; epsilontIs random disturbance; p is ARThe order; q is the order of MA;
defining a delay operator B, BXt＝Xt-1Then the polynomial of the AR coefficient phip(B)＝1-φ1B-…-φp(B)pCoefficient of MA polynomial thetaq(B)＝1-θ1B-…-θq(B)q；
Introduction of a difference operator Δd＝(1-B)dThen the ARIMA (p, d, q) model is expressed as:
Φp(B)ΔdXt＝Θ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(Bs)ΔdΔs DXt＝Θq(B)ΘQ(Bs)εt
wherein epsilontIs 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，ΘQ(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 FDA0003216564430000021
Seasonal Q-order coefficient
Figure FDA0003216564430000022
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 (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 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 FDA0003216564430000023
wherein f isSARIMA() Is SARIMA (P, D, Q) x (P, D, Q)sThe functional expression of the model is shown as,
Figure FDA0003216564430000024
small period training set, T, for the ith iterationtraiFor a small period of training, TforeThe 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 FDA0003216564430000025
is a time sequence predicted by the ith iteration of the SARIMA model;
calculating the predicted mean value of the ith iteration
Figure FDA0003216564430000026
Figure FDA0003216564430000027
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:
Figure FDA0003216564430000031
Figure FDA0003216564430000032
wherein
Figure FDA0003216564430000033
Is at 1-alphaP.IFor the value of the z-distribution of significance,
Figure FDA00032165644300000314
is the upper bound of the dynamic threshold interval for the ith small period,
Figure FDA0003216564430000034
is the lower bound of the dynamic threshold interval,
Figure FDA0003216564430000035
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 FDA0003216564430000036
wherein
Figure FDA0003216564430000037
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 FDA0003216564430000038
And training sequence
Figure FDA0003216564430000039
There is a real-time correlation, the predicted sequence of the ith minor period
Figure FDA00032165644300000310
Can be estimated as:
Figure FDA00032165644300000311
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
Figure FDA00032165644300000312
The time series judgment calculation method of (2) is as follows:
Figure FDA00032165644300000313
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 FDA0003216564430000041
When the value of the flow rate is in the normal interval range
Figure FDA0003216564430000042
Wherein
Figure FDA0003216564430000043
Suppose that the actual corresponding time of occurrence of an ICS anomaly is a sequence
Figure FDA0003216564430000044
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
Figure FDA0003216564430000045
Comprises the following steps:
Figure FDA0003216564430000046
wherein t isdebugFor real-time debugging before program runningM, ysampIs the sampling frequency of the ICS flow;
the nth ICS abnormality occurs at
Figure FDA0003216564430000047
The number of small-period iterations in the sub-small-period iteration calculation, for which the ICS abnormal event occurs
Figure FDA0003216564430000048
The following equation system can be used to obtain the calculation:
Figure FDA0003216564430000049
wherein KnIs an intermediate variable;
variance of dynamic ICS flow threshold interval generated by SARIMA online detection algorithm
Figure FDA00032165644300000410
The following algorithm can be used:
Figure FDA00032165644300000411
wherein
Figure FDA00032165644300000412
Can be aligned with
Figure FDA00032165644300000413
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 FDA00032165644300000414
The calculation method comprises the following steps:
Figure FDA00032165644300000415
wherein when an element is present
Figure FDA00032165644300000416
A slight ICS traffic anomaly in time,
Figure FDA00032165644300000417
for significant ICS traffic anomalies, when
Figure FDA00032165644300000418
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 CN112187528A (en) 2021-01-05
CN112187528B true 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)

Families Citing this family (5)

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

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107733905A (en) * 2017-10-24 2018-02-23 北京威努特技术有限公司 A kind of detection method of industry control network unit exception flow

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105704103B (en) * 2014-11-26 2017-05-10 中国科学院沈阳自动化研究所 Modbus TCP communication behavior abnormity detection method based on OCSVM double-contour model
KR101621019B1 (en) * 2015-01-28 2016-05-13 한국인터넷진흥원 Method for detecting attack suspected anomal event
CN107547269B (en) * 2017-08-14 2020-06-30 浙江大学 Method for constructing intelligent substation communication flow threshold model based on FARIMA

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107733905A (en) * 2017-10-24 2018-02-23 北京威努特技术有限公司 A kind of detection method of industry control network unit exception flow

Also Published As

Publication number Publication date
CN112187528A (en) 2021-01-05

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
CN103824137B (en) A kind of complex mechanical equipment multi-state failure prediction method
CN112799898B (en) Interconnection system fault node positioning method and system based on distributed fault detection
CN107016236A (en) Power network false data detection method for injection attack based on non-linear measurement equation
CN113191485B (en) Power information network security detection system and method based on NARX neural network
CN107517205B (en) Intelligent substation network abnormal flow detection model construction method based on probability
CN113239132A (en) Online out-of-tolerance identification method for voltage transformer
CN118157132B (en) Data mining method and device for voltage monitoring system based on neural network
Li et al. Deep learning based covert attack identification for industrial control systems
CN113269327A (en) Flow anomaly prediction method based on machine learning
CN118130984A (en) Cable partial discharge fault real-time monitoring method based on data driving
Akbarian et al. Attack resilient cloud-based control systems for industry 4.0
Horvath et al. Sensor fault diagnosis of inland navigation system using physical model and pattern recognition approach
Li et al. Real-time detecting false data injection attacks based on spatial and temporal correlations
CN116796261A (en) Closed switch equipment mechanical characteristic prediction method based on artificial intelligence
CN112085043A (en) Intelligent monitoring method and system for network security of transformer substation
Aftabi et al. A Variational Autoencoder Framework for Robust, Physics-Informed Cyberattack Recognition in Industrial Cyber-Physical Systems
CN111966966B (en) Method and system for analyzing feasible domain of sensor measurement error model parameters
CN114995338A (en) Industrial process micro-fault detection method based on normative variable analysis and JS divergence fusion
Cao et al. An adaptive UKF algorithm for process fault prognostics
Wang et al. Remaining useful life prediction method by integrating two-phase accelerated degradation data and field information

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