Industrial control system anomaly detection method combining time-frequency analysis and pure consistency measurement

Abstract

The invention discloses an industrial control system anomaly detection method combining time-frequency analysis and pure consistency measurement, and belongs to the technical field of industrial control system monitoring and intelligent diagnosis. The problems that in the prior art, when industrial control system data are processed, precision is not high, robustness is poor and the like are solved, discrete fractional Fourier transform (DFRFT) is adopted to extract time-frequency features, important features are screened in combination with an energy extraction strategy, and complex nonlinear and non-stationary data are effectively processed. On the basis, inverse fractional Fourier transform (iDFRFT) is used for recovering time domain features, and the inherent rule of the data is further disclosed. Meanwhile, through K-S inspection and random consistency analysis, the consistency and stability of feature distribution are quantified, and the accuracy and reliability of feature extraction are ensured. According to the method, a normal behavior mode is modeled by using a generative adversarial auto-encoding network (AGAE), and potential abnormal behaviors are identified according to reconstruction errors, so that the accuracy and robustness of anomaly detection are improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of industrial control system monitoring and intelligent diagnosis, and particularly relates to an abnormal detection method for industrial control systems combining time-frequency analysis and pure consistency measurement. Background Art

[0002] With the rapid development of industrial automation and informatization, industrial control systems (ICS) play a crucial role in modern production processes. It is widely used in multiple fields such as power, chemical industry, metallurgy, transportation, and intelligent manufacturing, supporting the normal operation of various important facilities. However, with the continuous expansion and complexity of industrial control systems, the operational safety and stability of the systems are facing more and more challenges. Abnormalities and faults in the system may lead to production interruptions, equipment damage, or even serious safety accidents. Therefore, real-time monitoring and fault diagnosis of industrial control systems have important practical significance.

[0003] Existing abnormal detection methods for industrial control systems usually rely on traditional statistical analysis, model identification, or machine learning techniques. Although these methods have improved the accuracy of abnormal detection to a certain extent, they have certain limitations in the following aspects:

[0004] (1) Complexity of data features: The data in industrial control systems usually has characteristics such as high dimension, multi-scale, non-linearity, and noise, making it difficult for traditional methods to effectively extract the potential laws of the data. For example, there may be singular data or abnormal patterns in the data, and these abnormal patterns often have large differences from normal data, making it difficult for traditional methods to accurately capture.

[0005] (2) Insufficient time-frequency information: Most traditional methods mainly rely on time-domain or frequency-domain features for analysis, but ignore the rich information that can be provided by the combination of time-frequency domains. In industrial control systems, abnormal behaviors often show various change characteristics of signals in the time domain and frequency domain. Therefore, traditional time-domain or frequency-domain analysis methods are difficult to comprehensively reflect the system state.

[0006] (3) Lack of consistency measurement: Existing abnormal detection methods usually lack in-depth quantitative evaluation of the consistency of data feature distributions, especially for non-linear and heterogeneous data that may exist in complex systems, whose distribution characteristics are often more complex. The detection of such data is easily affected by data deviation and noise, thus affecting the reliability and accuracy of the detection results. Summary of the Invention

[0007] Aiming at the problems of low precision and poor robustness existing in the prior art when dealing with industrial control system data, the present invention provides an abnormal detection method for industrial control systems combining time-frequency analysis and pure consistency measurement.

[0008] To achieve the above object, the present invention adopts the following technical solutions:

[0009] An abnormal detection method for industrial control systems that combines time-frequency analysis and pure consistency measurement, the method comprising the following steps:

[0010] Step 1: Time-frequency feature extraction: For the singular data in the industrial control system, first use the discrete fractional Fourier transform (DFRFT) to perform frequency-domain analysis on each dimension of features and extract frequency-domain features. The DFRFT can flexibly adjust the transformation order, and thus extract multi-level frequency-domain features, which is especially suitable for processing non-stationary, non-linear or noisy data. During the process of extracting frequency-domain features, by setting different orders, the information of the system data at multiple scales is obtained, revealing the potential laws and abnormal patterns of the singular data;

[0011] Further, the specific operation of the said Step 1 is:

[0012] For each dimension of feature, calculate the fractional Fourier transform of different orders respectively, and the formula is as follows:

[0013]

[0014] where, F α (u) represents the result of the fractional Fourier transform, where α is the order of the fractional order and u is the transformed variable; f(x) represents the input time-domain signal; represents the quadratic phase factor in the fractional Fourier transform, where j is the imaginary unit and x is the time-domain variable; e 2πjux represents the complex exponential factor in the Fourier transform. These symbols together constitute the mathematical expression of the fractional Fourier transform for analyzing and processing signals.

[0015] For the discrete signal x[n], assume that x[n] is a signal with a length of N. The calculation formula for the discrete fractional Fourier transform is:

[0016]

[0017] where, X α [k] represents the result of the discrete fractional Fourier transform, and k is the frequency-domain index after transformation; K α (n,k) is the discrete fractional kernel function, which is the discretized version of the continuous kernel function and is defined as:

[0018]

[0019] where, α is the fractional order, representing the rotation angle between the time domain and the frequency domain, and the range is 0 ≤ α ≤ 2π. Usually, the order α is expressed in the form of α = p * π / 2, where p is any real number, representing the fraction of the Fourier transform;

[0020] When α = 0, the DFRFT degenerates to the original signal, i.e.:

[0021] X0[k] = x[k]

[0022] When α = π / 2, the DFRFT degenerates to the standard discrete Fourier transform, i.e.:

[0023]

[0024] When α = π, the DFRFT degenerates to a time-reversed signal, i.e.:

[0025] X π [k] = x[N - k]

[0026] When α = 2π, the DFRFT returns to the original signal again;

[0027] Through the transformation, time-frequency characteristics at different orders are obtained, and these characteristics can help us capture potential anomalies in the industrial control system traffic.

[0028] Step 2: Energy extraction and data feature screening: Based on the frequency-domain feature extraction, use the energy extraction strategy to screen the extracted frequency-domain features, and select the part with the cumulative energy accounting for more than 80% of the total energy. This strategy helps to eliminate redundant information, extract the most representative features for anomaly detection, thereby reducing the computational complexity while retaining the core features of the data.

[0029] Furthermore, the specific operation of the said Step 2 is as follows:

[0030] Extract the part containing more than 80% energy from the fractional Fourier transform result, filter out the low-energy noise, and extract the energy part with a proportion of more than 80% by calculating the energy distribution of the transformed signal. The formula is:

[0031] E(f) = ∫|F α (u)| 2 du

[0032]

[0033] Among them, E(f) represents the total energy of the signal, which is obtained by integrating the squared modulus of the fractional Fourier transform result; F α (u) represents the result of the fractional Fourier transform; E p represents the cumulative energy; the cumulative energy reaches the required energy percentage p, which satisfies where p = 0.8.

[0034] Step 3: Time-domain Restoration: The inverse fractional Fourier transform (iDFRFT) is used to restore the extracted frequency-domain features back to the time domain. This process reveals the inherent time-domain features of the data, compensates for the deficiencies of time-frequency domain analysis, and provides a richer data basis for subsequent anomaly detection. Through the joint analysis of time-frequency features, abnormal changes in the system can be more comprehensively captured, improving the detection accuracy.

[0035] Further, the specific operation of Step 3 is as follows:

[0036] After extracting the energy, perform the inverse fractional Fourier transform on it to restore the time-frequency signal to the time domain, obtaining the industrial control system data after noise elimination:

[0037]

[0038] Step 4: Consistency Measurement and Feature Stability Analysis: Kolmogorov-Smirnov (K-S) test and random consistency analysis are used to ensure the stability and consistency of features. The K-S test is used to quantify the difference in feature distributions between singular and non-singular data, evaluating the stability of data features; random consistency analysis deeply analyzes the distribution consistency of features to ensure that the selected features have strong representativeness and consistency. This method effectively avoids the influence of data distribution deviation on the detection results and improves the reliability of feature extraction.

[0039] Further, the specific operation of Step 4 is as follows:

[0040] Use the K-S test to compare the distribution similarity between singular and non-singular data in the industrial control system. The K-S test calculates the maximum difference between the two distributions by comparing their empirical cumulative distribution functions (ECDFs):

[0041] D = max|F train (x) - F test (x)|

[0042] where D represents the maximum difference between the two empirical cumulative distribution functions (ECDFs), used to measure the similarity between the two distributions; F train (x) represents the empirical cumulative distribution function of the training data, where x is the data point; F test (x) represents the empirical cumulative distribution function of the test data;

[0043] Perform random consistency adjustment on the K-S statistic of each feature. The adjusted K-S statistic is defined as:

[0044] KS′ f = KS f - γ f (t)

[0045] Among them, KS f represents the K-S statistic of the i-th feature after random consistency adjustment, and γ f (t) is the random consistency statistic obtained through fractional Fourier transform, energy extraction, and inverse transformation, representing the impact of the random consistency noise between singular data and non-singular data on the i-th feature;

[0046] Perform multiple random samplings on each feature, generate the corresponding K-S test statistic, and construct its reference distribution. Compare the K-S statistic of the actual data with the random reference distribution and calculate the difference:

[0047] D adj = D actual - D random

[0048] Among them, D adj represents the adjusted difference, that is, the difference between the K-S statistic of the actual data and the random reference distribution; D actual represents the K-S statistic of the actual data, that is, the maximum difference between the empirical cumulative distribution function of the actual data and the reference distribution; D random represents the K-S statistic of the random reference distribution, that is, the maximum difference between the reference distribution generated by random sampling and the actual data distribution.

[0049] Measure the stability of the feature through the difference value and sort according to the difference value, and retain the feature with the smallest distribution difference;

[0050] After the test, obtain the K-S statistic of each feature, which reflects the distribution difference of the feature between singular data and non-singular data.

[0051] Step 5: Anomaly detection model training and evaluation: Based on the results of random consistency analysis, select a stable feature set, use this feature set as the input of the generative adversarial autoencoder network AGAE, and train the generative adversarial autoencoder network; AGAE can effectively learn and model normal behavior patterns, so as to identify abnormal behaviors with large differences from the normal pattern. Through the reconstruction error analysis of the features of the test data, accurately judge potential abnormal behaviors; when the reconstruction error exceeds a certain threshold, it is determined as abnormal.

[0052] In the present invention, the Generative Adversarial Autoencoder (GAAE) is used for anomaly detection in industrial control systems, combined with time-frequency analysis and pure consistency metrics to improve the accuracy and robustness of detection. The present invention utilizes the Generative Adversarial Autoencoder in the following manner: The objective of the autoencoder stage is to perform high-quality reconstruction of normal data through an encoder and a decoder. The encoder maps the input data into the latent space, while the decoder reconstructs the data from the latent space. The training objective is to minimize the reconstruction error to ensure that the model can accurately capture the feature distribution of normal data.

[0053] In the adversarial training stage, the present invention introduces an adversarial mechanism to enhance the robustness of the model. Specifically, one autoencoder attempts to mislead another autoencoder through the generated latent variables, making it difficult to distinguish between real data and reconstructed data. At the same time, the objective of the other autoencoder is to distinguish whether the input data is the original data or the reconstructed data. This adversarial process forces the model to generate reconstructed data closer to the real data, thereby improving the generalization ability of the model. The present invention integrates the objectives of the autoencoder stage and the adversarial training stage into a unified optimization objective. Through this unified optimization objective, the present invention can simultaneously optimize the reconstruction error and the adversarial loss to ensure that the model performs excellently in the anomaly detection task.

[0054] The training process of the Generative Adversarial Autoencoder is divided into two stages:

[0055] (1) Autoencoder stage: The main objective of training the encoder E and the decoder D is to minimize the reconstruction error so that the model can perform high-quality reconstruction of normal data; the training objective is:

[0056] L AE1 = ||X - D1(E(X))|| 2

[0057] L AE2 = ||X - D2(E(X))|| 2

[0058] where D1 and D2 are two decoders; E is the shared encoder; X is the input data;

[0059] (2) Adversarial training stage: The model is further trained through the adversarial mechanism; AE1 attempts to mislead AE2 through the generated latent variable Z so that it cannot distinguish between real data and the reconstructed data of AE1; AE2 attempts to distinguish whether the input data is the original data or the reconstructed data of AE1; the objective function is:

[0060] min AE1 max AE2 = ||X - D2(E(X))||2

[0061] In the two-stage training, the objectives of AE1 and AE2 are integrated into a unified optimization objective. The training objective of AE1 is:

[0062]

[0063] The training objective of AE2 is:

[0064]

[0065] where n is the number of training iterations, which is used to dynamically adjust the weights of the reconstruction error and the adversarial error;

[0066] The trained generative adversarial autoencoder network is used for anomaly detection. The final result of anomaly detection is determined by calculating the anomaly score. The formula for the anomaly score is:

[0067] A(Y) = α||Y - D2(E(Y))|| 2 + β||Y - D2(D1(E(Y)))|| 2

[0068] where Y represents the test set data, E represents the encoder, and D1 and D2 represent two different decoders; the anomaly score combines the errors of two parts. α||Y - D2(E(Y))|| 2 comes from the direct reconstruction error to detect whether the data deviates from the normal distribution; β||Y - D2(D1(E(Y)))|| 2 further enhances the discrimination ability of anomaly detection by re-reconstructing the discriminator's reconstruction result. α and β are weighting coefficients, satisfying α + β = 1, which are used to balance the contributions of the two parts of the error to the anomaly score; when α > β, it focuses on the direct reconstruction error and reduces false positives, which is suitable for scenarios with high requirements for false positives; when α < β, it focuses on the anomaly detection ability of the adversarial detection and improves the detection rate of subtle anomalies, which is suitable for scenarios with high requirements for the anomaly detection rate;

[0069] Calculate the reconstruction error of each sample in the validation set, and select an appropriate threshold τ such that the samples with errors higher than this threshold are determined to be anomalies; the threshold is set as the mean of the reconstruction error plus a constant multiple of the standard deviation:

[0070]

[0071] where k is a tuning parameter, and the optimal value is selected through cross-validation, represents the reconstruction error of the samples in the validation set.

[0072] Anomaly detection decision: For each test sample X″ i , calculate its reconstruction error:

[0073]

[0074] When then X′ i is determined to be abnormal.

[0075] Step 6: Comprehensive evaluation and output of anomaly detection results: Based on training and testing, combined with the detection ability of the model, comprehensively evaluate the anomaly detection results. Finally, output the detection results, and evaluate the adaptability and effectiveness of the model according to the detection accuracy and robustness.

[0076] Compared with the prior art, the present invention has the following advantages:

[0077] (1) Application of fractional Fourier transform: By performing fractional Fourier transforms of different orders on the data of the industrial control system, fully utilize the time-frequency domain information, extract important features in the traffic, and enhance the capture ability of the detection system for complex anomaly patterns.

[0078] (2) Energy extraction and inverse transform: After the fractional Fourier transform, extract signals in a specific energy range and perform inverse transform, retain the important information of the original data, and at the same time filter out noise and irrelevant features.

[0079] (3) K-S test and random consistency analysis: Use the Kolmogorov-Smirnov (K-S) test to compare the distribution similarities of training and test data, combined with random consistency analysis, screen out features with stability and robustness, and avoid the overfitting problem in feature selection.

[0080] (4) Combination of generative adversarial autoencoder network: Input the screened stable features into the autoencoder model for training and testing, improve the accuracy, recall rate, and F1 score of anomaly detection in the industrial control system, and ensure that the model can accurately identify and respond to attacks and abnormal situations. Description of the Drawings

[0081] Figure 1 is the time-frequency energy distribution and cumulative energy threshold graph of the embodiment;

[0082] Figure 2 is the comparison graph of the ECDF curves of Feature 1 of the embodiment;

[0083] Figure 3 is the comparison graph of the ECDF curves of Feature 2 of the embodiment;

[0084] Figure 4 is the comparison graph of the ECDF curves of Feature 3 in the embodiment;

[0085] Figure 5Bar chart of the K-S statistic for the features in the embodiment

[0086] Figure 6 Schematic diagram of the reconstruction of the generative adversarial autoencoder network in the embodiment;

[0087] Figure 7 Overall framework diagram of the model in the embodiment;

[0088] Figure 8 Training and validation loss diagram in the embodiment;

[0089] Figure 9 Diagram of the index results obtained using the traditional autoencoder in the embodiment;

[0090] Figure 10 Diagram of the index results obtained using the K-S statistic as feature selection in the embodiment;

[0091] Figure 11 Diagram of the index results obtained using the feature selection method of the present invention in the embodiment. Detailed implementation manners

[0092] To gain a deep understanding of the present invention, we will describe it comprehensively and meticulously. However, the present invention has various implementation manners and is not limited to the specific examples listed herein. The presentation of these examples aims to deepen the comprehensive understanding of the disclosed content of the present invention.

[0093] The present invention provides an industrial control system anomaly detection method combining time-frequency analysis and pure consistency measurement. Taking the SWaT dataset as an example, the SWaT (Secure Water Treatment) dataset is a publicly available dataset specifically for researching and testing the security and anomaly detection of industrial control systems. It was created by the Singapore University of Technology and Design (SUTD) and is based on their real water treatment system experimental platform. The SWaT dataset is very valuable in researching anomaly detection, intrusion detection, and fault diagnosis in industrial control systems. We use the first 80% as the training set and the last 20% as the test set. The main steps are as follows:

[0094] Step 1: Collect flow data and perform preprocessing. Remove outliers and perform data cleaning to ensure the consistency and comparability of the input flow data. In the data preprocessing stage, the original data was manually processed, and some unimportant feature columns were removed. For the case of missing numerical data, the average value of the nearest two data was used for filling; for the missing non-numerical features, the entire data was deleted.

[0095] Step 2: Perform fractional Fourier transform (FRFT) of order 0.1 to 0.9 on the industrial control system data and extract time-frequency features. In this step, we first obtain the system flow data and perform fractional Fourier transform (FRFT) on each feature. Through the fractional Fourier transform, both the time-domain and frequency-domain information of the data can be retained, and then the time-frequency features of the features can be extracted.

[0096] For each feature, calculate the fractional Fourier transform with orders ranging from 0.1 to 0.9 respectively. The formula is as follows:

[0097]

[0098] where F α (u) represents the result of the fractional Fourier transform, where α is the order of the fraction, and u is the transformed variable; f(x) represents the input time-domain signal; represents the quadratic phase factor in the fractional Fourier transform, where j is the imaginary unit and x is the time-domain variable; e 2πjux represents the complex exponential factor in the Fourier transform. These symbols together constitute the mathematical expression of the fractional Fourier transform for signal analysis and processing.

[0099] For the discrete signal x[n], assume x[n] is a signal of length N. The calculation formula for the discrete fractional Fourier transform is:

[0100]

[0101] where X α [k] represents the result of the discrete fractional Fourier transform, and k is the frequency-domain index after transformation; K α (n,k) is the discrete fractional kernel function, which is the discretized version of the continuous kernel function and is defined as:

[0102]

[0103] where α is the order of the fraction, representing the rotation angle between the time domain and the frequency domain, with the range 0 ≤ α ≤ 2π. Usually, the order α is expressed in the form of α = p * π / 2, where p is any real number representing the fraction of the Fourier transform.

[0104] When α = 0, the DFRFT degenerates to the original signal, that is:

[0105] X0[k] = x[k]

[0106] When α = π / 2, the DFRFT degenerates to the standard discrete Fourier transform (DFT), that is:

[0107]

[0108] When α = π, the DFRFT degenerates into a signal with time reversal, that is:

[0109] X π [k] = x[N - k]

[0110] When α = 2π, the DFRFT returns to the original signal again.

[0111] Through the transformation, time-frequency characteristics at different orders are obtained, and these characteristics can capture potential anomalies in the industrial control system traffic.

[0112] Step 3, from the fractional Fourier transform results, extract the part containing more than 80% of the energy, filter out the low-energy noise, and by calculating the energy distribution of the transformed signal, extract the part with more than 80% of the energy. The formula is:

[0113] E(f) = ∫|F α (u)| 2 du

[0114]

[0115] Among them, E(f) represents the total energy of the signal, which is obtained by integrating the square of the modulus of the fractional Fourier transform result; F α (u) represents the result of the fractional Fourier transform; E p represents the cumulative energy; the cumulative energy reaches the required energy percentage p, which satisfies where p = 0.8.

[0116] Then select the frequency range that makes the energy accumulation reach 80%. After extracting the energy, perform the inverse fractional Fourier transform (IFRFT) on it to restore the time-frequency signal to the time domain, and obtain the industrial control system data after noise elimination:

[0117]

[0118] Step 4, use the Kolmogorov - Smirnov (K - S) test to compare the distribution similarity between the singular data and non - singular data in the industrial control system. The K - S test calculates the maximum difference between the two distributions by comparing their empirical cumulative distribution functions (ECDFs):

[0119] D = max|F train (x) - F test (x)|

[0120] Among them, D represents the maximum difference between the two empirical cumulative distribution functions (ECDFs), which is used to measure the similarity between the two distributions; F train (x) represents the empirical cumulative distribution function of the training data, where x is the data point; Ftest (x) represents the empirical cumulative distribution function of the test data;

[0121] To consider the influence of randomness in the data on the difference in feature distributions, the K-S statistic for each feature is adjusted for random consistency. The adjusted K-S statistic is defined as:

[0122] KS′ f = KS f -γ f (t)

[0123] where KS f represents the K-S statistic of the i-th feature after random consistency adjustment, and γ f (t) is the random consistency statistic obtained through fractional Fourier transform, energy extraction, and inverse transformation, representing the influence of the random consistency noise between singular and non-singular data on the i-th feature;

[0124] Step 5: Obtain the distribution of the K-S test results through random consistency analysis, calculate the stability of each feature, and rank them. Random consistency analysis is used to evaluate the stability of the K-S test results. A reference distribution is generated through multiple random samplings and compared with the K-S statistic of the actual data.

[0125] Perform multiple random samplings for each feature, generate the corresponding K-S test statistic, and construct its reference distribution. Compare the K-S statistic of the actual data with the random reference distribution and calculate the difference:

[0126] D adj = D actual -D random

[0127] where D adj represents the adjusted difference, that is, the difference between the K-S statistic of the actual data and the random reference distribution; D actual represents the K-S statistic of the actual data, that is, the maximum difference between the empirical cumulative distribution function of the actual data and the reference distribution; D random represents the K-S statistic of the random reference distribution, that is, the maximum difference between the reference distribution generated through random sampling and the actual data distribution.

[0128] Measure the stability of the feature by the difference value and rank according to the difference value, retaining the feature with the smallest distribution difference;

[0129] After the test, the K-S statistic of each feature is obtained, reflecting the distribution difference of the feature between singular and non-singular data.

[0130] The random consistency RC i represents the feature Xi Statistical noise generated under random factors. The results of the following theorem are cited here, indicating that the impact of random noise can be reduced through random consistency adjustment, thereby improving the stability of the feature selection process. The specific formula is as follows:

[0131]

[0132] where represents the feature distribution of the f-th time series feature in the k-th detection in the industrial control system anomaly detection, represents the probability loss function of the f-th time series feature conforming to the anomaly detection in the k-th detection, L * then represents the optimal value of the loss function for the industrial control system anomaly detection.

[0133] Let q(T f ) represent the adjustment function of the f-th time series feature. Then for any q(T f ), we have:

[0134] 1 - RA = p + (1 - 2p)q(T f ) ≤ 1 - p

[0135] Therefore,

[0136] L = (1 - RA)PL ≤ (1 - p)PL

[0137] Further magnifying the upper bound therein will yield:

[0138] L * ≤ L(hT f ) ≤ (1 - p)PL(hT f )

[0139] hT f represents the construction or model state under the practice feature T f . For all distributions, assume then the pure loss of hT f satisfies:

[0140]

[0141] Explanation of the adjusted K-S test result. The K-S statistic KS′ f (t) represents the true distribution difference between the training set and the test set on feature i after eliminating the random noise. According to the formula:

[0142]

[0143] KS′ f(t) represents the true distribution difference of the f-th time series feature at time t after eliminating random noise. In addition, the following explanations can be made:

[0144] If KS′ f (t) has a large value, it indicates that there is a large difference in the distribution of feature f between the training set and the test set, and it is an unstable feature.

[0145] If KS′ f (t) has a small value, it indicates that the distribution of feature f is relatively consistent between the training set and the test set, and it is a stable feature.

[0146] The random consistency adjustment can significantly reduce the generalization error, and its formula is as follows:

[0147]

[0148] In the time series data scenario, T f represents the set of the f-th time series features in anomaly detection, and γ f (t) represents the influence of random noise of feature f at time t. The above formula shows that as the sample size N increases, the adjusted KS′ f (t) will be closer to the true KS f(t) , that is, it has better feature distribution stability, so as to effectively reduce the interference brought by random noise and ensure the robustness of feature selection. The proof is as follows:

[0149] First, according to the union bound to handle the upper limit of probability:

[0150]

[0151] PL N (T f ) represents the probability loss function of the f-th time series feature in a certain detection. ∈ represents the threshold within the error range, measuring the fluctuation of the probability loss, and |H| represents the size of the feature set, that is, the number of features in the time series data. Then for any we can get:

[0152]

[0153] According to the bound derivation, the probability is at most 1 - δ / 4, and the following estimate can be obtained:

[0154]

[0155] Let Then,

[0156]

[0157] By introducing Rc(H), the noise impact on features is adjusted to reduce the interference of noise on the detection performance and ensure the stability of the feature distribution. The adjustment of the coefficient αε can make the probability loss of the feature gradually converge to the true distribution difference in the case of a large sample N.

[0158] For the second term of the above formula, by the triangle inequality, we can get:

[0159]

[0160] For each function with probability at least 1 - δ / 4, we can get:

[0161]

[0162] Also, because there is Then for each function with probability at least 1 - δ / 4, we can get:

[0163]

[0164] The above represents the true and estimated probability correlation of the f-th time series feature, which is used to evaluate the similarity between features. L(T f ) represents the loss function of the f-th time series feature. Substituting δ, we can get that for any there is:

[0165]

[0166] Combining multiple formulas gives the final result.

[0167] According to the adjusted K-S statistic KS′ i , we can sort the features from small to large and select the features with smaller KS′ i as stable features:

[0168] sorted(KS′1, KS′2, …, KS′ n )

[0169] where n is the total number of features, and the sorted features can be used for anomaly detection in subsequent models. According to the following formula, we can get how the random consistency adjustment reduces the noise in the feature distribution:

[0170]

[0171] This formula indicates that if the difference between the estimated values and the expected values of features Z1 and Z2 exceeds a certain threshold ∈, then through random consistency adjustment, this difference can be decomposed into two parts to control the noise of each feature separately, thereby reducing the overall impact.

[0172] For the estimated values Z1, Z2 ∈ [0, 1] of two features in anomaly detection, there is a threshold ε ∈ (0, 1], which represents the allowable error range and measures whether the difference between the ratio and the expected value exceeds this threshold. So that We can obtain:

[0173]

[0174] β and γ represent the coefficients in the decomposition, which are used to further refine the control of random variables and the application of inequalities. The first inequality is obtained through the following conditions:

[0175]

[0176] It is used to decompose the joint probability distribution into two parts to control the noise effects of a and b separately.

[0177] The second inequality uses conditional probability for decomposition and is obtained through the following conditions:

[0178]

[0179] And the noise effects are further restricted by controlling events B1 and B2.

[0180] For any events B1, B2, where is the naivety of B2, we will take it as B2 to divide the two terms of the first inequality, so that We can obtain:

[0181]

[0182] Finally, by assuming ε ≤ 1, the final result γ ≥ α is obtained. In the anomaly detection of industrial control systems, this formula and the derivation process show how to effectively reduce the impact of noise on features through conditional probability and noise adjustment. By adjusting β and γ, the estimation error between features can be ensured to be controlled within a reasonable range, thereby improving the detection accuracy and robustness.

[0183] According to the following derivation, prove that the feature selection method has good generalization ability under specific conditions:

[0184]

[0185] As the sample set increases, the influence of random noise in the feature selection process will be further reduced, ensuring the generalization ability of the model on different data sets. Here, τ = 6.0 is selected. The training set and test set of these features have similar distributions, and after random consistency adjustment, they can provide more stable inputs for the anomaly detection model.

[0186] Step 6: Design a deep learning model, including a generator and a discriminator. Based on the results of the random consistency analysis in Step 5, select the features with the best stability, that is, the features with the highest distribution similarity, as the input of the generative adversarial autoencoder network.

[0187] Select the features with the highest stability from the sorting results. These features have the smallest distribution difference between the training set and the test set, and have high stability and representativeness. Input these features into the generative adversarial autoencoder network for training and testing.

[0188] Encoder: Consists of multiple convolutional layers, used to extract high-level features. After each convolutional layer, a non-linear activation function (ReLU) and a pooling layer are connected to reduce the dimension.

[0189] Hyperparameter settings:

[0190] Convolution kernel size: 3x3; Number of convolutional layers: 3 layers; Pooling window: 2x2; Activation function: ReLU.

[0191] Decoder: Consists of multiple deconvolutional layers, used to reconstruct the low-dimensional representation into the original feature dimension. After each deconvolutional layer, a non-linear activation function (such as ReLU) is connected.

[0192] Hyperparameter settings:

[0193] Deconvolution kernel size: 3x3; Number of deconvolutional layers: 3 layers; Activation function: ReLU.

[0194] Step 7: Use the standardized feature matrix to train the generative adversarial autoencoder network. The goal is to minimize the reconstruction error between the input feature matrix and the reconstructed feature matrix. The mean squared error (MSE) is used as the loss function, and the model parameters are updated through the backpropagation algorithm, using the Adam optimizer.

[0195] Hyperparameter settings:

[0196] Learning rate: 0.01; Batch size: 32; Number of training epochs: 50.

[0197] Step 8: Use the trained generative adversarial autoencoder network to predict the singular data and detect anomalies by comparing the reconstruction errors.

[0198] Threshold determination: By calculating the distribution of reconstruction errors on the validation set, a reasonable threshold τ is set to distinguish normal data from abnormal data. The specific method is to calculate the reconstruction error of each sample in the validation set and select an appropriate threshold such that samples with errors higher than this threshold are determined to be abnormal. The threshold is set as the mean of the reconstruction errors plus a constant multiple of the standard deviation:

[0199]

[0200] where k is a tuning parameter and the optimal value is selected through cross-validation.

[0201] Anomaly detection decision: For each test sample X″ i , calculate its reconstruction error

[0202]

[0203] If then X′ i is determined to be abnormal.

[0204] Step 9: And the corresponding bar chart is drawn to visualize the results. Then, according to the evaluation results, the model parameters are optimized to ensure robust anomaly detection in various situations, while considering the comprehensive performance of calculating precision, recall, accuracy, and F1-score. According to the evaluation results, the model parameters are optimized to ensure robust anomaly detection in various situations.

[0205] Step 10: Deploy the optimized model to the actual industrial control system environment, monitor the traffic data in real time, and perform anomaly detection. For the detected abnormal behaviors, corresponding security response measures are taken in a timely manner, such as blocking the traffic, issuing an alarm, etc., to ensure the security of the system.

[0206] It is verified through experimental data and trained and tested on the SWaT dataset. According to Figure 8 It can be seen that due to the large training data, when the number of training rounds reaches 10, the training loss and validation loss have dropped to a very small value and tend to be flat. The algorithm of the present invention has the advantages that with the increase of the number of training rounds, the accuracy steadily increases, and the main patterns and features of the data can be learned within a relatively small number of rounds.

[0207] Figure 9The results of ablation experiments using traditional autoencoders (without adversarial training, i.e., β=0β=0) are shown. From the results, the accuracy of the model is 0.91, indicating that it performs well in the overall classification task and can correctly classify most samples. However, the recall rate is 0.67, indicating that the model has limited ability to identify positive (abnormal) samples, resulting in a considerable number of abnormal samples being misclassified as normal. The precision rate is 0.60, indicating that only 60% of the samples predicted as positive are truly positive, and there are many false positives. The F1 score is 0.63, reflecting the poor performance of the model in balancing precision and recall. Overall, although the traditional autoencoder performs well in overall classification accuracy without adversarial training, it performs poorly in detecting abnormal samples. This shows that adversarial training plays an important role in improving the model's ability to identify abnormal samples. In order to further improve the recall and precision of the model, it may be necessary to introduce adversarial training or other enhancement techniques to improve the model's ability to detect abnormal samples.

[0208] Figure 10 The results of using only the KS test as the feature selection criterion are shown, ignoring the mitigation of random consistency in feature selection. As can be seen from the figure, the accuracy, recall, precision, and F1 score of the model are relatively low, which indicates that relying solely on the KS test for feature selection may not be sufficient to effectively improve model performance. Although the KS test can help identify differences between feature distributions, it ignores the random consistency in the feature selection process, which may lead to the selection of unstable features, thereby affecting the stability and generalization ability of the model. In order to improve the effect of feature selection, it may be necessary to combine other methods or techniques to alleviate the random consistency problem, such as by integrating multiple feature selection methods or introducing regularization techniques to improve the overall performance of the model.

[0209] Figure 11 The results of the patented method are shown, with an accuracy of 0.95, a recall of 0.59, a precision of 1.0, and an F1 score of 0.74. These results show that the patented method performs very well in terms of overall classification accuracy and precision, especially when the precision reaches 1.0, indicating that all samples predicted by the model as positive are truly positive and there are no false positives. However, the recall rate is 0.59, indicating that the model still has room for improvement in identifying all positive samples, and some abnormal samples are missed. Nevertheless, the F1 score is 0.74, showing that the model performs well in balancing precision and recall. Overall, the patented method works best in anomaly detection tasks, especially in reducing false positives and improving classification accuracy. In order to further improve the recall rate, it may be necessary to further optimize the model or introduce additional techniques to enhance the ability to identify abnormal samples.

[0210] The content not detailedly described in the specification of the present invention belongs to the prior art well-known to those skilled in the art. Although the illustrative specific embodiments of the present invention are described above for the understanding of those skilled in the art of the present technology, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

Claims

1. An industrial control system anomaly detection method combining time-frequency analysis with pure consistency measurement, characterized in that: The method comprises the following steps: Step 1: Time-frequency feature extraction: For the singular data in the industrial control system, the discrete fractional Fourier transform (DFRFT) is first used to perform frequency domain analysis on each dimensional feature to extract the frequency domain features. In the process of extracting the frequency domain features, different orders are set to obtain the information of the system data at multiple scales, revealing the potential laws and abnormal patterns of the singular data; Step 2: Energy extraction and data feature screening: Based on the frequency domain feature extraction, the energy extraction strategy is used to screen the extracted frequency domain features to screen out the part whose cumulative energy accounts for more than 80% of the total energy; Step 3: Time domain recovery: Use the inverse fractional Fourier transform iDFRFT to restore the extracted frequency domain features back to the time domain; Step 4: Consistency measurement and feature stability analysis: KS test and random consistency analysis are used to ensure the stability and consistency of the features; Step 5: Anomaly detection model training and evaluation: Based on the results of random consistency analysis, a stable feature set is selected, the feature set is input into the generative adversarial autoencoder network AGAE, and the generative adversarial autoencoder network is trained; the potential abnormal behavior is accurately judged through reconstruction error analysis of test data features; when the reconstruction error exceeds a certain threshold, it is judged as abnormal; Step 6: Comprehensive evaluation and output of anomaly detection results: Based on training and testing, the anomaly detection results are comprehensively evaluated in combination with the detection capabilities of the model. Finally, the detection results are output and the adaptability and effect of the model are evaluated based on the detection accuracy and robustness.

2. The method for detecting anomalies in an industrial control system by combining time-frequency analysis with pure consistency measurement according to claim 1 is characterized in that: The specific operations of step 1 are: For each dimension of features, fractional Fourier transforms of different orders are calculated respectively. The formula is as follows: Among them, F α (u) represents the result of fractional Fourier transform, where α is the order of the fractional order and u is the transformed variable; f(x) represents the input time domain signal; represents the quadratic phase factor in the fractional Fourier transform, where j is the imaginary unit and x is the time domain variable; e 2πjux Represents the complex exponential factor in Fourier transform, which is used to analyze and process signals; For a discrete signal x[n], assuming x[n] is a signal of length N, the calculation formula for the discrete fractional Fourier transform is: Among them, X α [k] represents the result of discrete fractional Fourier transform, k is the frequency domain index after the transform; K α (n,k) is a discrete fractional kernel function, which is a discretized version of the continuous kernel function and is defined as: Where α is the fractional order, representing the rotation angle between the time domain and the frequency domain, ranging from 0≤α≤2π; usually, the order α is expressed as α=p*π / 2, where p is an arbitrary real number representing the fraction of the Fourier transform; When α = 0, DFRFT degenerates into the original signal, that is: X0[k]=x[k] When α = π / 2, DFRFT degenerates into the standard discrete Fourier transform, that is: When α=π, DFRFT degenerates into a time-reversed signal, that is: X π [k]=x[N-k] When α=2π, DFRFT returns to the original signal; through transformation, time-frequency characteristics of different orders are obtained.

3. The method for detecting anomalies in an industrial control system by combining time-frequency analysis with pure consistency measurement according to claim 2 is characterized in that: The specific operations of step 2 are: From the fractional Fourier transform result, extract the part containing more than 80% of the energy, filter out the low-energy noise, and extract the energy part accounting for more than 80% by calculating the energy distribution of the transformed signal. The formula is: E(f)=∫|F α (un)| 2 you Where E(f) represents the total energy of the signal, which is obtained by integrating the square of the modulus of the fractional Fourier transform result; F α (u) represents the result of fractional Fourier transform; E p Indicates the cumulative energy; the cumulative energy reaches the required energy percentage p to satisfy Where p = 0.

8.

4. The method for detecting anomalies in an industrial control system by combining time-frequency analysis with pure consistency measurement according to claim 3 is characterized in that: The specific operations of step 3 are: After extracting the energy, perform an inverse fractional Fourier transform on it to restore the time-frequency signal to the time domain, and obtain the industrial control system data after noise elimination:

5. The method for detecting anomalies in an industrial control system by combining time-frequency analysis with pure consistency measurement according to claim 4 is characterized in that: The specific operations of step 4 are: The KS test is used to compare the distribution similarity of singular data and non-singular data in industrial control systems. The KS test calculates the maximum difference between the two distributions by comparing their empirical cumulative distribution functions ECDF: D=max|F train (x-F test (x| Where D represents the maximum difference between two empirical cumulative distribution functions (ECDFs), which is used to measure the similarity of two distributions; F train (x) represents the empirical cumulative distribution function of the training data, where x is a data point; F test (x) represents the empirical cumulative distribution function of the test data; The KS statistic of each feature is adjusted for random consistency, and the adjusted KS statistic is defined as: KS′ f =KS f -γ f (t) Among them, KS f represents the KS statistic of the i-th feature after random consistency adjustment, γ f (t) is the random consistency statistic obtained by fractional Fourier transform, energy extraction and inverse transform, which represents the influence of random consistency noise between singular data and non-singular data on the i-th feature; Perform multiple random sampling on each feature, generate the corresponding KS test statistic, and construct its reference distribution. Compare the KS statistic of the actual data with the random reference distribution and calculate the difference: D adj =D actual -D random Among them, D adj represents the adjusted difference, that is, the difference between the KS statistic of the actual data and the random reference distribution; D actual The KS statistic represents the actual data, that is, the maximum difference between the empirical cumulative distribution function of the actual data and the reference distribution; D random The KS statistic represents the random reference distribution, that is, the maximum difference between the reference distribution generated by random sampling and the actual data distribution; The stability of the feature is measured by the difference value, and the features are sorted according to the difference value, and the features with the smallest distribution difference are retained; After the test, the KS statistic of each feature is obtained, which reflects the distribution difference of the feature on singular data and non-singular data.

6. The method for detecting anomalies in an industrial control system by combining time-frequency analysis with pure consistency measurement according to claim 5, characterized in that: The specific operations of step 5 are: The training process of the generative adversarial autoencoder network is divided into two stages: (1) Autoencoder stage: The main goal of training the encoder E and decoder D is to minimize the reconstruction error so that the model can reconstruct normal data with high quality. The training objectives are: IT AE1 =||X-D1(E(X))|| 2 L AE2 =||X-D2(E(X))|| 2 Where D1 and D2 are two decoders; E is the shared encoder; X is the input data; (2) Adversarial training stage: The model is further trained through an adversarial mechanism; AE1 attempts to mislead AE2 through the generated latent variable Z, making it unable to distinguish between the real data and the reconstructed data of AE1; AE2 attempts to distinguish whether the input data is the original data or the reconstructed data of AE1; the objective function is: min AE1 max AE2 =||X-D2(E(X))|| 2 In the two-stage training, the objectives of AE1 and AE2 are integrated into a unified optimization objective. The training objective of AE1 is: The training objectives of AE2 are: Where n is the number of training iterations, which is used to dynamically adjust the weights of the reconstruction error and the adversarial error; The trained generative adversarial autoencoder network is used for anomaly detection. The final result of anomaly detection is determined by calculating the anomaly score. The formula for the anomaly score is: A(Y)=α||Y-D2(E(Y))|| 2 +β||Y-D2(D1(E(Y)))|| 2 Among them, Y represents the test set data, E represents the encoder, D2 and D2 represent two different decoders; the anomaly score combines the errors of the two parts, α||Y-D2(E(Y))|| 2 From direct reconstruction error, detect whether the data deviates from the normal distribution; β||Y-D2(D1(E(Y)))|| 2 Through the reconstruction of the discriminator reconstruction results, the adversarial mechanism is further used to enhance the ability of anomaly detection to distinguish. α and β are weighted coefficients, satisfying α+β=1, which are used to balance the contribution of the two parts of the error to the anomaly score. When α>β, the focus is on the direct reconstruction error; when α<β, the focus is on the anomaly detection ability of adversarial detection. Calculate the reconstruction error of each sample in the validation set, and select an appropriate threshold τ so that samples with errors higher than the threshold are judged as abnormal; the threshold is set to the mean of the reconstruction error plus a constant times the standard deviation: Among them, k is a tuning parameter, and the best value is selected through cross-validation. Represents the reconstruction error of the samples in the validation set; Anomaly detection judgment: For each test sample X″ i , calculate its reconstruction error: when Then X′ i Deemed abnormal.

Images