Non-stationary industrial process anomaly detection method based on slow feature decomposition and copula high-dimensional space prediction

Abstract

This invention relates to a method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction, belonging to the field of industrial process time-series anomaly detection technology. The method includes: using slow feature analysis to decompose industrial time-series data into fast and slow changing components; mapping the fast and slow features to a linearly observable Koopman high-dimensional space based on Koopman theory, constructing a deep prediction model to achieve stepwise prediction of the feature sequence; constructing a SPE statistic that approximately follows a weighted chi-square distribution by combining the deviation between the prediction results and the actual values, and setting SPE control limits for anomaly detection accordingly; introducing KL divergence to measure the distribution difference between the current operating condition and the normal operating condition, and adaptively updating the SPE control limits based on its magnitude; and during the online deployment phase of the model, calculating the SPE value in real time and comparing it with the control limits to achieve rapid identification and dynamic early warning of anomalies under non-stationary conditions.

Description

Technical Field

[0001] This invention belongs to the field of industrial process time-series anomaly detection technology, and relates to a non-stationary industrial process anomaly detection method based on slow feature decomposition and Koopman high-dimensional space prediction. Background Technology

[0002] With the increasing intelligence level in industrial manufacturing, the massive amounts of time-series data generated in real time during production bring new opportunities and challenges to anomaly detection. By deeply mining this data, abnormal patterns can be identified in a timely manner, ensuring the safety and stability of the production process. However, due to the complexity and variability of actual industrial environments, how to quickly and accurately detect anomalies in massive non-stationary time series has become a key research direction for industrial intelligence. The core objective of anomaly detection is to identify abnormal patterns that are significantly different from normal patterns. Traditional methods mainly rely on statistical analysis techniques, judging anomalies through indicators such as mean and variance. However, due to the high dimensionality, dynamism, and non-stationarity of industrial data, traditional methods have limitations in dealing with complex dynamic changes. Existing anomaly detection methods can be divided into three categories: the first category is statistical analysis-based methods (such as the 3σ method and box plots), which rely on the assumption of normal distribution and are difficult to adapt to nonlinear and non-stationary characteristics; the second category is clustering-based machine learning methods (such as DBSCAN and K-means), which cluster normal data and identify a small number of data points that deviate from the group as anomalies, but they are not effective for complex non-stationary data; the third category is deep learning-based methods, which learn the patterns of normal working conditions to make predictions or reconstructions and judge anomalies based on errors. Although these methods improve detection accuracy to a certain extent, they rely on a large amount of labeled data, while industrial anomaly samples are scarce, and the non-stationarity and variable working conditions of the data make model training difficult and prone to overfitting, thus affecting the accuracy and stability of detection. Currently, anomaly detection in non-stationary time series data in industry faces two major challenges: (1) Due to the limited number of collected data variables, industrial data is usually located in a low-dimensional space, making it difficult to capture highly nonlinear relationships between variables. Therefore, it is difficult to fully characterize the complexity of dynamic systems, resulting in limited model prediction capabilities; (2) Industrial data has multi-scale variation characteristics. Directly projecting all variables into a high-dimensional space will ignore the differences in the frequency of fast and slow changes, thereby reducing prediction accuracy and fault detection sensitivity. To address these issues, a new approach is proposed based on a combination of Slow Feature Decomposition (SFA) and Koopman high-dimensional space prediction representation. SFA can decompose slow and fast change features from non-stationary time series, achieving effective separation of multi-scale characteristics. Meanwhile, Koopman operator theory maps low-dimensional nonlinear dynamic systems to a high-dimensional linear space, enabling linear dynamic evolution and prediction in the high-dimensional space. Therefore, this method constructs a unified framework that can extract fast and slow change features and perform linear modeling in a high-dimensional space, effectively addressing the non-stationarity and dynamic complexity of data in industrial processes, and significantly improving the accuracy and real-time performance of anomaly detection. This method has broad application prospects in intelligent manufacturing and industrial intelligent operation and maintenance. Summary of the Invention

[0003] In view of this, the purpose of this invention is to provide a nonstationary industrial process anomaly detection method based on slow feature decomposition and Koopman high-dimensional space prediction, which addresses the problem of poor anomaly detection performance caused by the multi-scale variation characteristics of nonstationary industrial time series and the nonlinear and nonstationary characteristics in low-dimensional space.

[0004] To achieve the above objectives, the present invention provides the following technical solution:

[0005] A method for detecting anomalies in nonstationary industrial processes based on slow eigenvalue decomposition and Koopman high-dimensional space prediction includes the following steps:

[0006] Step 1: Collect industrial time-series data from the industrial process, and divide it into normal operating condition data and various variable operating condition data. Each variable operating condition data is first operated normally for a period of time, then an anomaly occurs, and the abnormal operation is maintained for a period of time before returning to normal operation. Divide the normal operating condition data into training set, test set and validation set. Preprocess all collected industrial time-series data.

[0007] Step 2: Feature engineering based on slow feature analysis, decompose the data in the normal operating condition training set into two parts of features: fast-changing and slow-changing features;

[0008] Step 3: Train a deep prediction neural network based on SFA-Koopa using normal operating condition training data. The SFA-Koopa-based deep prediction neural network contains multiple layers of SFA-Koopman prediction blocks. The SFA-Koopa-based deep prediction neural network uses a stackable prediction structure with multi-level residual correction to correct and compensate for the reconstruction residuals of the previous layer. The SFA-Koopman prediction blocks are used to map the fast-changing and slow-changing features to an observable high-dimensional linear space using Koopman theory, and make predictions using local and global Koopman operators respectively. The prediction results of each layer are used as the input features of the next layer. Finally, the prediction values ​​of each layer are added together to obtain the final prediction result.

[0009] Step 4: Test the trained SFA-Koopa-based deep prediction neural network using the test set data from the normal operating condition data, calculate the prediction error between the prediction result and the true value, and calculate the SPE control limit for anomaly detection based on the fact that the SPE control limit approximately follows a weighted chi-square distribution.

[0010] Step 5: Use the KL divergence between the training data and test data under normal operating conditions as a threshold; calculate the KL divergence between the normal operating data and the test data under normal operating conditions in various variable operating conditions, and compare it with the threshold to determine whether to update the SPE control limit.

[0011] Step 6: In the online detection phase, the data of the new working condition to be detected is input into the trained SFA-Koopa deep prediction neural network to predict the working condition data in real time, calculate the prediction error between the prediction result and the actual value, and then obtain the real-time SPE value. This value is compared with the SPE control limit to determine in real time whether an anomaly has occurred.

[0012] Furthermore, in step 1, the data from normal operating conditions is divided into a training set, a test set, and a validation set in a 7:2:1 ratio; assuming the industrial time-series data has J measurement variables and K sampling points, a K×J matrix is ​​formed; using X... n (N n ×J) and X a (N a ×J) respectively describe the normal operating condition dataset and the variable operating condition dataset; N n and N a These represent the number of samples in the normal dataset and the variable operating condition dataset, respectively; for dataset X... n and X a The following maximum-minimum normalization method is used:

[0013]

[0014] Where, x k,j This represents the data at the k-th time step of variable j in the dataset, where k∈[1,K] is the sampling time and j∈[1,J] is the variable index. x represents k,j The normalized value, Let $\mathbf{j}$ and $\mathbf{j}$ represent the maximum and minimum values ​​of the $j$-th variable in the training set of the normal operating condition dataset, respectively.

[0015] Then, an average sliding window filter is used to remove the effects of noise:

[0016]

[0017] Where N represents the time step of the sliding window. This represents the data of the j-th variable at time ti after normalization. This represents the data at time t after the j-th variable has undergone average moving filter.

[0018] Furthermore, in step 2, the feature engineering based on slow feature analysis specifically includes the following steps:

[0019] Training set data from preprocessed normal operating condition data Perform nonlinear mapping

[0020]

[0021] Calculate the rate of change of the features over time after nonlinear mapping.

[0022]

[0023] Calculate the covariance matrix C0 of the characteristic signal and the covariance matrix C of the time derivative:

[0024] C0=<Φ(X)Φ(X) T >

[0025]

[0026] The optimization objective is to minimize the rate of change of the output features over time.

[0027]

[0028] The constraints are:

[0029] W T C0W=I

[0030] Where W is the projection matrix, which is obtained by solving the generalized eigenvalue problem:

[0031] CW = C0WΛ

[0032] The smaller the eigenvalue Λ, the slower the corresponding feature changes;

[0033] Based on the scree plot, the eigenvectors corresponding to the smaller eigenvalues ​​before the mutation of eigenvalue Λ are selected to form the slow feature projection matrix W1, and the eigenvectors corresponding to the larger eigenvalues ​​after the mutation are selected to form the fast feature projection matrix W2. The slow feature projection matrix W1 is used to project the nonlinear mapping matrix Φ(X) to output the slow feature Y1(t); the fast feature projection matrix W2 is used to project the nonlinear mapping matrix Φ(X) to output the fast feature Y2(t), thus realizing the decomposition of the fast and slow changing parts.

[0034]

[0035] Furthermore, step 3 specifically includes the following steps:

[0036] Step 3.1: Using the observation function of the MLP encoder structure, map the slowly changing part Y1 obtained from the decomposition to the Koopman high-dimensional linear space. The observation function of the MLP encoder structure is:

[0037] Z1 = MLP Encoder(Y1)

[0038] Utilizing the globally learnable Koopman operator K slow Predicting the slowly changing components obtained from the decomposition:

[0039]

[0040] This represents the predicted value of slow features in high-dimensional space.

[0041] Then, the predicted values ​​in the high-dimensional space are mapped back to the original data space using the inverse observation function of the decoder structure to obtain the predicted values ​​for the slowly changing part:

[0042]

[0043] Step 3.2: Divide the rapidly changing part Y2 obtained from the decomposition into segments of length S, with an input time step of T. The segmented data is y. j :

[0044]

[0045] The segmented, rapidly changing data is mapped to the Koopman high-dimensional linear space using the observation function of the MLP encoder structure:

[0046] z i =MLP Encoder(y i )

[0047] Using the extended dynamic mode decomposition (eDMD) method, the local Koopman operator K, which adapts to rapid dynamic changes, is calculated by the least squares method. fast Using K fast Reconstructing and predicting known input data yields high-dimensional reconstructed values ​​and high-dimensional predicted values:

[0048]

[0049] Where T represents the number of input time steps, This represents the high-dimensional reconstruction value of the known input data. Represents the predicted value in a high dimension

[0050] Then, the predicted values ​​in the high-dimensional space are mapped back to the original data space through the inverse observation function of the decoder structure:

[0051]

[0052] Step 3.3: Use the stackable prediction structure with multi-level residual correction to correct and compensate the reconstruction residual of the previous layer, and add the prediction values ​​of the slow-changing part and the fast-changing part of the Koopman predictor to obtain the final prediction output.

[0053] Furthermore, in step 3.2, the local Koopman operator K, adapted to rapid dynamic changes, is calculated using the extended dynamic mode decomposition (eDMD) method via the least squares method. fastSpecifically, it includes:

[0054]

[0055] Z back This represents the high-dimensional representation of the first T / S-1 segment, Z. fore This represents the high-dimensional representation from the second paragraph to the last paragraph. It is Z back The Moore-Penrose pseudoinverse matrix is ​​obtained by solving for K using the least squares method. fast The value of .

[0056] Furthermore, in the training process of a deep predictive neural network based on SFA-Koopa, the Encoder, Decoder, and K... slow These are learnable parameters that share the Encoder and Decoder in both the slowly changing and rapidly changing parts, denoted as φ. slow and φ fast The model training uses the mean squared error loss function (MSE) to calculate the predicted value Y and the true value Y. gt The mean square error between the two values ​​is used to optimize the learnable parameter K by minimizing this error. slow φ slow and φ fast :

[0057]

[0058] Furthermore, in step 4, the trained SFA-Koopa-based deep prediction neural network is tested using test set data from the normal operating condition data. The prediction results obtained from the normal operating condition data are used to calculate the error between the prediction results and the true values. Based on the fact that the SPE control limit approximately follows a weighted chi-square distribution, the SPE control limit used for anomaly detection is calculated.

[0059]

[0060] g = v / 2m

[0061] h = 2m 2 / v

[0062]

[0063] Where, x i,k y represents the true value of the k-th variable in the i-th row of the original input data. i,k Let represent the predicted value of the k-th variable in the i-th row, m represent the mean of the SPE of all samples, v represent the variance of the SPE, and h represent the degrees of freedom parameter of the weighted chi-square distribution. It is the critical value of the chi-square distribution, with a confidence level of α.

[0064] Further, in step 5, the KL divergence, representing the difference in distribution between the two data sets, is calculated using the training and test data under normal operating conditions. This divergence serves as the threshold for determining whether the SPE control limits need to be updated. Then, under varying operating conditions, the KL divergence is calculated using the normal operating data and test data under the varying operating conditions, and compared with the KL divergence threshold. If it is greater than the threshold, the SPE control limits are updated using the data from the new operating conditions, following the method in step 4. If it is less than the threshold, no update is performed. The formula for calculating the KL divergence is as follows:

[0065]

[0066] Where P(i) is the probability of distribution P at point i, and Q(i) is the probability of distribution Q at point i.

[0067] Furthermore, in step 6, during the online detection phase, the variable operating condition data is processed through a trained SFA-Koopa deep prediction neural network to obtain predicted values, and the prediction error is calculated by comparing these predicted values ​​with the actual values, thereby obtaining the real-time SPE value.

[0068]

[0069] The real-time SPE value is compared with the SPE control limit. If the real-time SPE value is greater than the SPE control limit, it is judged as abnormal. If the real-time SPE value is less than the SPE control limit, it is judged as normal.

[0070] The beneficial effects of this invention are as follows: The anomaly detection method for non-stationary industrial processes based on the combination of slow feature decomposition and Koopman high-dimensional space provides an efficient and intelligent solution for anomaly detection in non-stationary time series data during industrial manufacturing. This method not only effectively extracts the fast and slow change features of industrial processes with minimal prior knowledge and human intervention, but also constructs a high-precision linear dynamic prediction model by mapping low-dimensional nonlinear dynamic systems to a high-dimensional linear space. This achieves accurate separation of multi-scale dynamic change characteristics and accurate prediction of the mapping from low-dimensional nonlinear space to high-dimensional linear space. Furthermore, this method exhibits good generalization performance under various industrial operating conditions, enabling rapid and accurate identification of abnormal patterns. It provides technical support for intelligent monitoring and fault early warning of complex industrial processes, while also helping to improve the stability and reliability of equipment operation and reduce industrial operation and maintenance costs.

[0071] Other advantages, objectives, and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination, or may be learned from practice of the invention. The objectives and other advantages of the invention can be realized and obtained through the following description. Attached Figure Description

[0072] To make the objectives, technical solutions, and advantages of the present invention clearer, the preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings, wherein:

[0073] Figure 1 Flowchart of the nonstationary industrial time series anomaly detection method based on slow feature decomposition and Koopman high-dimensional space prediction provided by the present invention;

[0074] Figure 2 The structure diagram of the SFA-Koopa deep prediction neural network model provided by this invention is shown below;

[0075] Figure 3 The following are the prediction results of the SFA-Koopa deep predictive neural network model described in this invention under normal operating conditions: (a) is a comparison between the actual and predicted RMS values ​​of accelerometer 1, (b) is a comparison between the actual and predicted RMS values ​​of accelerometer 2, (c) is a comparison between the actual and predicted current values, (d) is a comparison between the actual and predicted pressure values, (e) is a comparison between the actual and predicted temperature values, (f) is a comparison between the actual and predicted thermocouple values, (g) is a comparison between the actual and predicted voltage values, and (h) is a comparison between the actual and predicted RMS values ​​of volumetric flow rate.

[0076] Figure 4 The following are some abnormal detection results under varying operating conditions in the case described in this invention, where (a) is the exponential result of closing inlet valve 1, (b) is the δ function result of closing inlet valve 1, (c) is the exponential result of closing outlet valve 2, (d) is the δ function result of closing outlet valve 2, (e) is the abrupt behavior of rotor imbalance, (f) is the linear behavior of rotor imbalance, (g) is the step behavior of rotor imbalance, and (h) is the exponential behavior of rotor imbalance. Detailed Implementation

[0077] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. Unless otherwise specified, the following embodiments and features can be combined with each other.

[0078] It should be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. Therefore, the drawings only show the components related to the present invention and are not drawn according to the actual number, shape and size of the components in the actual implementation. In the actual implementation, the form, quantity and proportion of each component can be arbitrarily changed, and the layout of the components may also be more complex.

[0079] In the following description, numerous details are explored to provide a more thorough explanation of embodiments of the invention. However, it will be apparent to those skilled in the art that embodiments of the invention may be practiced without these specific details. In other embodiments, well-known structures and devices are shown in block diagram form rather than in detail to avoid obscuring embodiments of the invention.

[0080] Example 1:

[0081] The SKAB dataset (Skoltech Anomaly Benchmark) is a widely used benchmark dataset for detecting time-series anomalies in industrial equipment, primarily for detecting faults or abnormal behaviors in industrial equipment. The dataset records time-series data for the normal operation of a water circulation industrial process and 34 different types of abnormal operating conditions. Specifically, at a sampling frequency of 1Hz, it records eight variable signals, including vibration acceleration, motor current and voltage, water pump circuit pressure, engine temperature, fluid temperature in the circuit, and fluid flow rate. The abnormal operating condition types include 34 faults such as closing inlet valve 1, closing outlet valve 2, rotor imbalance (abrupt, linear, step, delta function, exponential) behavior, and increasing or decreasing water temperature in the tank.

[0082] This invention provides a method for detecting anomalies in nonstationary industrial processes based on slow eigenvalue decomposition and Koopman high-dimensional space prediction, comprising the following steps:

[0083] Step 1: Collect industrial time-series data from the industrial process, dividing it into normal operating condition data and variable operating condition data. In this embodiment, the normal operating condition data contains a total of 9405 data points, which are divided into 6584 training data points, 1881 test data points, and 940 validation data points in a 7:2:1 ratio. The variable operating condition data includes a different number of test data points for each condition, and the fault occurs after a period of normal operation, followed by a return to normal operation.

[0084] The segmented industrial time-series data is preprocessed. Specifically, min-max normalization is first performed using the maximum and minimum values ​​of the normal operating data.

[0085]

[0086] Where j = 1, 2, ..., m represents the number of variables in the industrial process; in this embodiment, j takes a maximum of 8; x k,j This represents the data at the k-th time step belonging to variable j. x represents k,j The normalized value; Let $\mathbf{j}$ and $\mathbf{j}$ represent the maximum and minimum values ​​of the $j$-th variable in the training set of the normal operating condition dataset, respectively.

[0087] Then, an average sliding window filter with a sliding window length N of 24 is used to remove the influence of noise.

[0088]

[0089] in This represents the data of the j-th variable at time ti after normalization. This represents the data at time t after the j-th variable has undergone average moving filter.

[0090] Step 2: Feature engineering based on slow feature analysis decomposes the data in the normal operating condition training set into two parts: fast-changing and slow-changing features. Specifically, firstly, the training set data from the preprocessed normal operating condition data... Perform nonlinear mapping:

[0091]

[0092] Then calculate the rate of change of the features after nonlinear mapping over time:

[0093]

[0094] Calculate the covariance matrix of the characteristic signal:

[0095] C0=<Φ(X)Φ(X) T >

[0096] Calculate the covariance matrix of the time derivative:

[0097]

[0098] The optimization objective is to minimize the rate of change of the output features over time.

[0099]

[0100] The constraints are:

[0101] W T C0W=I

[0102] Where W is the projection matrix, which is obtained by solving the generalized eigenvalue problem of CW = C0WΛ. The smaller the eigenvalue Λ, the slower the corresponding eigenvalue changes.

[0103] Based on the scree plot, the eigenvectors corresponding to the smaller eigenvalues ​​before the mutation of eigenvalue Λ are selected to form the slow feature projection matrix W1, and the eigenvectors corresponding to the larger eigenvalues ​​after the mutation are selected to form the fast feature projection matrix W2. The slow feature projection matrix W1 is then used to project the nonlinear mapping matrix Φ(X) to output the slow features. The nonlinear mapping matrix Φ(X) is projected using the fast feature projection matrix W2 to output the fast feature. To decompose the parts that change in speed.

[0104] Step 3: Train a deep prediction neural network based on SFA-Koopa using normal operating condition training data. The deep prediction neural network based on SFA-Koopa contains multiple layers of SFA-Koopman prediction blocks. The deep prediction neural network based on SFA-Koopa uses a stackable prediction structure with multi-level residual correction to correct and compensate for the reconstruction residuals of the previous layer. The SFA-Koopman prediction blocks are used to map the fast-changing and slow-changing features to an observable high-dimensional linear space using Koopman theory, and make predictions using local and global Koopman operators respectively. The prediction results of each layer are used as the input features of the next layer. Finally, the prediction values ​​of each layer are added together to obtain the final prediction result.

[0105] Specifically, in each SFA-Kopman prediction block, the slowly changing part Y1 obtained from the decomposition is mapped to the Koopman high-dimensional linear space using the observation function of the MLP encoder structure. The observation function of the MLP encoder structure is as follows:

[0106] Z1 = MLP Encoder(Y1)

[0107] Utilizing the globally learnable Koopman operator K slow Predicting the slowly changing portion:

[0108]

[0109] This represents the predicted value of slow features in high-dimensional space.

[0110] Then, the predicted values ​​in the high-dimensional space are mapped back to the original data space using the inverse observation function of the decoder structure to obtain the predicted values ​​for the slowly changing part:

[0111]

[0112] The rapidly changing portion Y2 obtained from the decomposition is first segmented, with a segment length of S, which is 12 in this embodiment. The input time step number T is 24 in this embodiment. The segmented data is as follows:

[0113] y j = [y1, y2] ∈ R 12×8

[0114] The segmented, rapidly changing data is mapped to the Koopman high-dimensional linear space using the observation function of the MLP encoder structure. The observation function of the encoder structure is:

[0115] z1 = MLP Encoder(y1)

[0116] z2 = MLP Encoder(y2)

[0117] Using the extended dynamic mode decomposition (eDMD) method, the local Koupman operator matrix K, which adapts to rapid dynamic changes, is obtained by solving the least squares method. fast :

[0118] Z back =[z1]

[0119] Z fore =[z2]

[0120]

[0121] Z back This represents the high-dimensional representation of the first T / S-1 segment, Z. fore This represents the high-dimensional representation from the second paragraph to the last paragraph. It is Z back The Moore-Penrose pseudoinverse matrix is ​​obtained by solving for K using the least squares method. fast The value of K; using K fast Reconstructing and predicting known input data yields high-dimensional reconstructed values ​​and high-dimensional predicted values:

[0122]

[0123] in, This represents the high-dimensional reconstruction value of the known input data. This represents the predicted value in a higher dimension.

[0124] Then, the predicted values ​​in the high-dimensional space are mapped back to the original data space using the inverse observation function of the decoder structure to obtain the predicted values ​​for the rapidly changing parts:

[0125]

[0126] The experimental setup for training this SFA-Koopa-based deep prediction neural network is as follows: SFA-Koopa prediction block layers = 3, activation function = ReLU, learning rate = 0.001, optimizer = Adam, batch size = 32, epochs = 50, input sequence length = 24, prediction length = 1. During model training, the Encoder, Decoder, and K... slow These are learnable parameters. To maintain the consistency of the Koopman embedding, the Encoder and Decoder are shared in both the slow-changing and fast-changing parts, denoted as φ. slow and φ fast The model training uses the mean squared error loss function (MSE) to calculate the predicted value Y and the true value Y. gt The mean square error between the two values ​​is used to optimize the learnable parameter K by minimizing this error. slow φ slow and φ fast :

[0127]

[0128] Step 4: Test the trained SFA-Koopa-based deep prediction neural network using the test set data from the normal operating condition data. The prediction results obtained from the normal operating condition data are as follows: Figure 3 As shown, (a) compares the actual and predicted RMS values ​​of accelerometer 1, (b) compares the actual and predicted RMS values ​​of accelerometer 2, (c) compares the actual and predicted current values, (d) compares the actual and predicted pressure values, (e) compares the actual and predicted temperature values, (f) compares the actual and predicted thermocouple values, (g) compares the actual and predicted voltage values, and (h) compares the actual and predicted RMS values ​​of volumetric flow rate. The prediction error between the predicted results and the actual values ​​is calculated, and based on the fact that the SPE control limit approximately follows a weighted chi-square distribution, the SPE control limit for anomaly detection is calculated.

[0129]

[0130] g = v / 2m

[0131] h = 2m 2 / v

[0132]

[0133] Where, x i,k y represents the true value of the k-th variable in the i-th row of the original input data. i,k Let represent the predicted value of the k-th variable in the i-th row, m represent the mean of the SPE for all samples, v represent the variance of the SPE, and h represent the degrees of freedom parameter of the weighted chi-square distribution. It is the critical value of the chi-square distribution, with a confidence level of α, which is specifically taken as 0.05 in this embodiment.

[0134] Step 5: Use the KL divergence between the training data and test data under normal operating conditions as a threshold; calculate the KL divergence between the normal operating data and the test data under normal operating conditions in various variable operating conditions, and compare it with the threshold to determine whether to update the SPE control limit.

[0135] If the KL divergence between the test data under the new operating condition and the test data under the normal operating condition is greater than the KL divergence threshold, then the SPE control limits need to be updated using partial data from the new operating condition according to the method in step 4; if it is less than the threshold, then no update is performed. Specifically, the formula for calculating the KL divergence between the two data distributions is as follows:

[0136]

[0137] Where P(i) is the probability of distribution P at point i, and Q(i) is the probability of distribution Q at point i.

[0138] Step 6, Online Detection Phase: The data of the new operating condition to be detected is input into the trained SFA-Koopa deep prediction neural network. The model will automatically make real-time predictions and calculate the prediction error with the true value, thereby obtaining the real-time SPE value.

[0139]

[0140] The real-time SPE value is compared with the SPE control limit obtained in step 5 to determine in real time whether an anomaly has occurred. If the real-time SPE value is greater than the SPE control limit, it is judged as an anomaly; if the real-time SPE value is less than the SPE control limit, it is judged as normal. This achieves the purpose of detecting anomalies in non-stationary industrial time series. Some of the obtained variable operating condition anomaly detection results are shown below. Figure 4As shown, (a) is the exponential result of closing inlet valve 1, (b) is the δ function result of closing inlet valve 1, (c) is the exponential result of closing outlet valve 2, (d) is the δ function result of closing outlet valve 2, (e) is the abrupt behavior of rotor imbalance, (f) is the linear behavior of rotor imbalance, (g) is the step behavior of rotor imbalance, and (h) is the exponential behavior of rotor imbalance.

[0141] Example 2:

[0142] An electronic device, comprising a memory and a processor;

[0143] The memory is used to store computer programs;

[0144] The processor is configured to implement the method described in Embodiment 1 when executing the computer program.

[0145] Example 3:

[0146] A computer-readable storage medium storing a computer program that, when executed by a processor, implements the method described in Embodiment 1.

[0147] Example 4:

[0148] A computer program product includes a computer program that, when executed by a processor, implements the method described in Example 1.

[0149] In the above embodiments, the reference to "this embodiment" in the specification indicates that a specific feature, structure, or characteristic described in connection with the embodiment is included in at least some embodiments, but not necessarily all embodiments. Multiple appearances of "this embodiment" do not necessarily refer to the same embodiment.

[0150] In the above embodiments, although the invention has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory structures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed. The embodiments of the invention are intended to cover all such substitutions, modifications, and variations falling within the broad scope of the appended claims.

[0151] As will be understood by those skilled in the art, the computer-readable storage medium described in this embodiment allows for the implementation of all or part of the steps in the above method embodiments by computer program-related hardware. The aforementioned computer program can be stored in a computer-readable storage medium. When executed, the program performs the steps of the above method embodiments; and the aforementioned storage medium includes various media capable of storing program code, such as ROM, RAM, magnetic disks, or optical disks.

[0152] The electronic terminal provided in this embodiment includes a processor, a memory, a transceiver, and a communication interface. The memory and the communication interface are connected to the processor and the transceiver and complete communication between them. The memory is used to store computer programs, the communication interface is used to perform communication, and the processor and the transceiver are used to run the computer programs, so that the electronic terminal performs the steps of the above method.

[0153] In this embodiment, the memory may include random access memory (RAM) and may also include non-volatile memory, such as at least one disk storage device.

[0154] The processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0155] This invention can be used in a wide range of general-purpose or special-purpose computing system environments or configurations. Examples include: personal computers, server computers, handheld or portable devices, tablet devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, and distributed computing environments including any of the above systems or devices, etc.

[0156] This invention can be described in the general context of computer-executable instructions, such as program modules, that are executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform a specific task or implement a specific abstract data type. This invention can also be practiced in distributed computing environments where tasks are performed by remote processing devices connected via a communication network. In distributed computing environments, program modules can reside in local and remote computer storage media, including storage devices.

[0157] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for detecting anomalies in non-stationary industrial processes based on slow eigenvalue decomposition and Koopman high-dimensional space prediction, characterized in that: Includes the following steps: Step 1: Collect industrial time-series data from the industrial process, and divide it into normal operating condition data and various variable operating condition data. Each variable operating condition data is first operated normally for a period of time, then an anomaly occurs, and the abnormal operation is maintained for a period of time before returning to normal operation. Divide the normal operating condition data into training set, test set and validation set. Preprocess all collected industrial time-series data. Step 2: Feature engineering based on slow feature analysis, decompose the data in the normal operating condition training set into two parts of features: fast-changing and slow-changing features; Step 3: Train a deep prediction neural network based on SFA-Koopa using normal operating condition training data. The deep prediction neural network based on SFA-Koopa contains multiple layers of SFA-Koopman prediction blocks. The deep prediction neural network based on SFA-Koopa uses a stackable prediction structure with multi-level residual correction to correct and compensate for the reconstruction residuals of the previous layer. The SFA-Koopman prediction blocks are used to map the fast-changing and slow-changing features to an observable high-dimensional linear space using Koopman theory, and make predictions using local Koopman operators and global Koopman operators respectively. The prediction result of each layer is the input feature of the next layer. Finally, the predicted values ​​of each layer are summed to obtain the final prediction result; Step 4: Test the trained SFA-Koopa-based deep prediction neural network using the test set data from the normal operating condition data, calculate the prediction error between the prediction result and the true value, and calculate the SPE control limit for anomaly detection based on the fact that the SPE control limit approximately follows a weighted chi-square distribution. Step 5: Use the KL divergence between the training data and the test data under normal operating conditions as a threshold; Calculate the KL divergence between normal operating data and normal operating condition test data in various variable operating condition operating data, and compare it with the threshold to determine whether to update the SPE control limit; Step 6: In the online detection phase, the data of the new working condition to be detected is input into the trained SFA-Koopa deep prediction neural network to predict the working condition data in real time, calculate the prediction error between the prediction result and the actual value, and then obtain the real-time SPE value. This value is compared with the SPE control limit to determine in real time whether an anomaly has occurred.

2. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 1, characterized in that: In step 1, the data from normal operating conditions are divided into training, testing, and validation sets in a 7:2:1 ratio; assuming the industrial time-series data has J measurement variables and K sampling points, a K×J matrix is ​​formed; using X... n (N n ×J) and X a (N a ×J) respectively describe the normal operating condition dataset and the variable operating condition dataset; N n and N a These represent the number of samples in the normal dataset and the variable operating condition dataset, respectively; for dataset X... n and X a The following maximum-minimum normalization method is used: Where, x k,j This represents the data at the k-th time step of variable j in the dataset, where k∈[1,K] is the sampling time and j∈[1,J] is the variable index. x represents k,j The normalized value, Let $\mathbf{j}$ and $\mathbf{j}$ represent the maximum and minimum values ​​of the $j$-th variable in the training set of the normal operating condition dataset, respectively. Then, an average sliding window filter is used to remove the effects of noise: Where N represents the time step of the sliding window. This represents the data of the j-th variable at time ti after normalization. This represents the data at time t after the j-th variable has undergone average moving filter.

3. The anomaly detection method for non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 1, characterized in that: Step 2, the feature engineering based on slow feature analysis, specifically includes the following steps: Training set data from preprocessed normal operating condition data Perform nonlinear mapping Calculate the rate of change of the features over time after nonlinear mapping. Calculate the covariance matrix C0 of the characteristic signal and the covariance matrix C of the time derivative: C0=<Φ(X)Φ(X) T > The optimization objective is to minimize the rate of change of the output features over time. The constraints are: IN T C0W=I Where W is the projection matrix, which is obtained by solving the generalized eigenvalue problem: CW = C0WΛ The smaller the eigenvalue Λ, the slower the corresponding feature changes; Based on the scree plot, the eigenvectors corresponding to the smaller eigenvalues ​​before the mutation of eigenvalue Λ are selected to form the slow feature projection matrix W1, and the eigenvectors corresponding to the larger eigenvalues ​​after the mutation are selected to form the fast feature projection matrix W2. The slow feature projection matrix W1 is used to project the nonlinear mapping matrix Φ(X) to output the slow feature Y1(t); the fast feature projection matrix W2 is used to project the nonlinear mapping matrix Φ(X) to output the fast feature Y2(t), thus realizing the decomposition of the fast and slow changing parts.

4. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 1, characterized in that: Step 3 specifically includes the following steps: Step 3.1: Using the observation function of the MLP encoder structure, map the slowly changing part Y1 obtained from the decomposition to the Koopman high-dimensional linear space. The observation function of the MLP encoder structure is: Z1 = MLP Encoder(Y1) Utilizing the globally learnable Koopman operator K slow Predicting the slowly changing components obtained from the decomposition: This represents the predicted value of slow features in high-dimensional space. Then, the predicted values ​​in the high-dimensional space are mapped back to the original data space using the inverse observation function of the decoder structure to obtain the predicted values ​​for the slowly changing part: Step 3.2: Divide the rapidly changing part Y2 obtained from the decomposition into segments of length S, with an input time step of T. The segmented data is y. j : The segmented, rapidly changing data is mapped to the Koopman high-dimensional linear space using the observation function of the MLP encoder structure: z i =MLP Encoder(y i ) Using the extended dynamic mode decomposition (eDMD) method, the local Koopman operator K, which adapts to rapid dynamic changes, is calculated by the least squares method. fast Using K fast Reconstructing and predicting known input data yields high-dimensional reconstructed values ​​and high-dimensional predicted values: Where T represents the number of input time steps, This represents the high-dimensional reconstruction value of the known input data. The high-dimensional predicted values ​​are then mapped back to the original data space through the inverse observation function of the decoder structure: Step 3.3: Use the stackable prediction structure with multi-level residual correction to correct and compensate the reconstruction residual of the previous layer, and add the prediction values ​​of the slow-changing part and the fast-changing part of the Koopman predictor to obtain the final prediction output.

5. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 4, characterized in that: Step 3.2 describes using the extended dynamic mode decomposition (eDMD) method to calculate the local Koopman operator K, which adapts to rapid dynamic changes, via the least squares method. fast Specifically, it includes: Z back This represents the high-dimensional representation of the first T / S-1 segment, Z. fore This represents the high-dimensional representation from the second paragraph to the last paragraph. It is Z back The Moore-Penrose pseudoinverse matrix K is obtained by solving for it using the least squares method. fast The value of .

6. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 4, characterized in that: In the training process of a deep predictive neural network based on SFA-Koopa, the Encoder, Decoder, and K... slow These are learnable parameters that share the Encoder and Decoder in both the slowly changing and rapidly changing parts, denoted as φ. slow and φ fast The model training uses the mean squared error loss function (MSE) to calculate the predicted value Y and the true value Y. gt The mean square error between the two values ​​is used to optimize the learnable parameter K by minimizing this error. slow φ slow and φ fast :

7. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 1, characterized in that: In step 4, the trained SFA-Koopa-based deep prediction neural network is tested using test set data from the normal operating condition data. The prediction results obtained from the normal operating condition data are used to calculate the error between the prediction results and the true values. Based on the fact that the SPE control limit approximately follows a weighted chi-square distribution, the SPE control limit used for anomaly detection is calculated. g = v / 2m h=2m 2 / v Where, x i,k y represents the true value of the k-th variable in the i-th row of the original input data. i,k Let represent the predicted value of the k-th variable in the i-th row, m represent the mean of the SPE of all samples, v represent the variance of the SPE, and h represent the degrees of freedom parameter of the weighted chi-square distribution. It is the critical value of the chi-square distribution, with a confidence level of α.

8. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 7, characterized in that: The KL divergence, representing the difference in distribution between the two sets of data, is calculated using training and test data under normal operating conditions. This divergence serves as the threshold for determining whether the SPE control limit needs to be updated. Then, under varying operating conditions, the KL divergence between the two sets of data is calculated using normal operating data and test data under normal operating conditions, and compared with the KL divergence threshold. If the value is greater than the threshold, the SPE control limits are updated using partial data from the new operating condition, following the method in step 4. If the value is less than the threshold, no update is performed. The formula for calculating the KL divergence is as follows: Where P(i) is the probability of distribution P at point i, and Q(i) is the probability of distribution Q at point i.

9. The method for detecting anomalies in non-stationary industrial processes based on slow feature decomposition and Koopman high-dimensional space prediction according to claim 8, characterized in that: In step 6, during the online detection phase, the variable operating condition data is processed through a trained SFA-Koopa deep prediction neural network to obtain predicted values, and the prediction error is calculated by comparing it with the actual values, thereby obtaining the real-time SPE value. The real-time SPE value is compared with the SPE control limit. If the real-time SPE value is greater than the SPE control limit, it is judged as abnormal. If the real-time SPE value is less than the SPE control limit, it is judged as normal.

Images