Physical information constraint embedded non-stationary industrial process anomaly detection method

Abstract

The invention relates to a physical information constraint embedded non-stationary industrial process anomaly detection method, and belongs to the technical field of industrial process time sequence anomaly detection. According to the method, a dynamic coupling relation between key variables is extracted through frequency domain lagging correlation, original time sequence data is decomposed into two subspaces, namely a slowly-varying trend subspace and a quickly-varying disturbance subspace, through slow feature analysis, and long-term stable change and transient disturbance change are modeled respectively. Furthermore, a high-dimensional linear evolution model is constructed in a fast space and a slow space by adopting a Kupman operator, so that the multi-scale dynamic modeling precision is effectively improved. In addition, two types of physical information supervision mechanisms are introduced in the model training process: based on a water pump flow-pressure second-order dynamic equation and a motor electric power conservation law, the physical consistency and engineering interpretability of the prediction process are significantly enhanced. According to the invention, abnormal working condition identification and early warning under variable working conditions in a non-stationary industrial process can be realized.

Description

Technical Field

[0001] This invention belongs to the field of industrial process time series anomaly detection technology, and relates to a multivariate time series modeling method that combines physical information constraints and Koopman space transformation, which is particularly suitable for anomaly perception and fault early warning under non-stationary operating conditions. Background Technology

[0002] As the complexity of industrial manufacturing systems continues to increase, the multivariate time-series data collected in real time during industrial processes exhibits highly dynamic and non-stationary characteristics, posing new challenges to anomaly detection. Accurately identifying system anomalies under varying operating conditions, high-frequency disturbances, and nonlinear coupling has become a crucial research direction for ensuring production safety and improving operational reliability. Traditional anomaly detection methods typically assume data stationarity and simple distribution patterns, making them ill-suited to the complex dynamic evolution of industrial processes. Especially in the early stages of a fault, anomalous signals are often masked by normal fluctuations, lacking in-depth modeling tools for the system's dynamic laws and physical mechanisms, leading to frequent false alarms and missed alarms. Current methods for detecting anomalies in industrial time series data can be broadly categorized into three types: The first type is based on statistical indicators, such as the 3σ criterion and box plots, which rely on normal distribution and fixed thresholds for judgment, but cannot handle non-stationary structures and lag dependencies in time series data. The second type is clustering methods, such as K-means, DBSCAN, and SVM, which are suitable for outlier detection in static feature spaces, but struggle to capture dynamic interactions and evolutionary trends between variables. The third type is based on deep learning prediction and reconstruction models, such as LSTM, Transformer, and AutoEncoder, which have the ability to learn complex dynamic processes, but often neglect physical mechanisms, have weak robustness to varying operating conditions, are insensitive to distribution drift, and are highly dependent on sample size and annotation quality. This makes model training difficult to converge in the context of scarce industrial anomaly data, resulting in unsatisfactory actual deployment effects. In non-stationary industrial scenarios, anomaly detection tasks face the following typical challenges: (1) There are cross-time delay and nonlinear coupling characteristics between variables, and it is difficult to accurately capture system behavior by relying solely on static correlation or a single time scale; (2) Multivariate time series data often contain both slow-changing trends and fast-changing disturbances. Unified modeling will mask their respective characteristics, making it difficult for the model to accurately characterize key dynamic information; (3) Industrial systems generally have clear physical structures and conservation relationships, such as the flow-pressure dynamics of water pumps and the power conservation of motors. Failure to include them in model supervision will weaken the physical rationality of the prediction. Summary of the Invention

[0003] In view of this, the purpose of this invention is to provide an anomaly detection method for non-stationary industrial processes with embedded physical information constraints. Addressing the problems of dynamic coupling between variables, mixed fast and slow features, nonlinear relationships, and lack of physical constraints in non-stationary industrial time series, this invention proposes an anomaly detection method that integrates variable correlation attention mechanism, slow feature analysis, and Koopman high-dimensional linear modeling. Furthermore, physical information supervision based on pump flow-pressure dynamics and the electric power conservation equation is introduced during model training to improve the physical consistency of the prediction process and the accuracy and adaptability of anomaly detection, thereby effectively addressing the needs for anomaly identification and early warning under complex and variable operating conditions in non-stationary industrial processes.

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

[0005] A method for detecting anomalies in non-stationary industrial processes with embedded physical information constraints, specifically including the following steps:

[0006] S1: Acquire multivariate time series data generated during operation in the industrial system, divide it into normal operating condition data and abnormal operating condition data containing faults or changing factors; divide the normal operating condition data into training set, validation set and test set according to the proportion, and perform normalization preprocessing on the original multivariate time series data;

[0007] S2: After input data preprocessing, a variable correlation attention mechanism is introduced. By performing a fast Fourier transform on multivariate time series data, the frequency domain cross-correlation coefficient between variables is calculated, and an attention weight matrix is ​​generated to reflect the correlation strength between variables, thereby enhancing the temporal interaction expression between important variables.

[0008] S3: Apply the slow feature analysis method to the attention-enhanced representation to decouple the multivariate time series data into subspaces with different response speeds, and extract the slow features that represent the long-term trend of the system and the fast features that characterize rapid disturbances respectively.

[0009] S4: Using the Koopman operator theory, fast and slow features are mapped to a high-dimensional observable linear space respectively. Slow features are modeled by the slow-changing Koopman operator, and fast features are modeled by extracting the fast-changing Koopman operator through the extended dynamic mode decomposition (eDMD) method. The model performance is then iteratively optimized layer by layer using a stacked residual structure, and finally the whole multivariate prediction results are output.

[0010] S5: During the training of the prediction model, a physical information equation is introduced as a supervision term to construct two types of physical residual losses: one is based on the second-order flow-pressure dynamic relationship of the pump system; the other is based on the energy conservation equation between motor voltage, current and output power, which is introduced into the total loss function to enhance the physical consistency and interpretability of the prediction.

[0011] S6: In the offline modeling stage, the trained prediction model is evaluated using a test set under normal working conditions. The squared prediction error (SPE) index is calculated based on the difference between the predicted value and the actual value, and an anomaly detection control limit based on its statistical characteristics is established using a weighted chi-square distribution.

[0012] S7: To address the distribution drift problem under changing operating conditions, KL divergence is used to measure the distribution difference between the current operating condition and the original normal operating condition, and based on this, it is determined whether to update the SPE control limit to improve the adaptability of anomaly detection in changing operating conditions.

[0013] S8: During the online detection phase, the real-time collected data is input into the model to obtain the predicted value and calculate the SPE value. The value is then compared with the current control limit to determine whether the current state is abnormal, thus achieving rapid fault detection and dynamic early warning.

[0014] Furthermore, in step S1, the collected multivariate time series data includes time series variables from multiple operating conditions. The normal operating condition data is divided into training set, validation set, and test set in a 7:2:1 ratio. Let the multivariate time series data have J measurement variables and T time steps, forming a T×j matrix. The normal operating condition data and the abnormal operating condition data are denoted as X... n (N n ×J) and X a (N a ×J), where N n and N a These represent the corresponding sample numbers; for dataset X n and X a Perform the following maximum and minimum value normalization:

[0015]

[0016] Where, x k,j This represents the original value of the j-th variable at time k. These are the maximum and minimum values ​​of variable j in the training set, respectively; This is the normalized value.

[0017] Furthermore, in step S2, a variable-related attention mechanism is introduced, specifically including: for the normalized sequence Perform inverted embedding, and then map it to a latent space using a multilayer perceptron:

[0018]

[0019] Where E represents the feature matrix, T represents the time step, J represents the number of variables, D represents the latent space dimension of the mapping, Permute represents dimension transpose, and MLP represents multilayer perceptron;

[0020] The obtained feature matrix E is reused to represent the query matrix Q, key matrix K, and value matrix V; its hysteresis cross-correlation in the frequency domain is calculated:

[0021]

[0022] Among them, S QK (f) represents the correlation between the Q and K frequency domains. denoted by Fourier transform, * denotes complex conjugate, and ⊙ denotes Hadamard product;

[0023] The inverse Fourier transform is applied to recover the time-domain cross-correlation matrix:

[0024]

[0025] Among them, R Q,K (τ) represents the correlation between Q and K in the time domain. Indicates the inverse Fourier transform;

[0026] Weighted summation of correlations over all lag steps τ:

[0027]

[0028] Wherein, COR(q) i k j ) represents the variable q i k j Total correlation across all lag steps Represents variable q i k j The correlation at lag step τ, λ τ The weighting coefficients for the learned lag step size;

[0029] The final attention weight matrix A is obtained by Softmax normalization. i，j :

[0030] A i，j =Softmax(COR(q) i k j ))

[0031] The value matrix features are fused using an attention weight matrix to obtain an enhanced representation, which is then used as input for the subsequent slow feature analysis module.

[0032] X VCA (t)=A i，j ·V

[0033] Among them, X VCA (t) represents the enhanced transvariate lag correlation representation based on attention mechanism weighting, and V represents the value matrix.

[0034] Furthermore, in step S3, the slow feature analysis method specifically includes: enhancing the variable representation output by the variable correlation attention mechanism module. By zero-meaning along the time axis, we obtain:

[0035] X c =X VCA -μ t (X VCA )

[0036] Among them, X c This represents the augmented representation of the centered variable, μ. t (·) indicates the mean along the time axis;

[0037] Calculate the sample covariance matrix ∑ x And perform feature decomposition on it:

[0038]

[0039] Among them, V x Represents the covariance matrix ∑ x eigenvector matrix, Λ x Represents the covariance matrix ∑ x The eigenvalue diagonal matrix;

[0040] Construct the whitening matrix and obtain the whitened data:

[0041]

[0042] Among them, W w Represents the whitening transformation matrix. This represents the whitened input.

[0043] Next, calculate the time derivative difference matrix:

[0044]

[0045] Where D represents the time derivative difference matrix;

[0046] And define the derivative covariance matrix C based on the time derivative difference matrix. d And perform feature decomposition:

[0047]

[0048] Among them, Λ d =diag(λ1,...,λ J The eigenvalues ​​are arranged in ascending order; V d Represents the time derivative covariance matrix;

[0049] The slow feature dimension k is determined based on the cumulative variance contribution rate, and the first k feature vectors are used to construct the slow feature projection matrix W. s Based on this, the slow feature Y is calculated by projection. slow :

[0050]

[0051] Then, the slow features are restored from the whitening space to the original input space using the inverse mapping, so as to perform fast and slow separation consistent with the original sequence:

[0052]

[0053] Among them, W o W represents the mapping matrix projected from the slow feature space back to the original input space. w Let X represent the whitening transformation matrix. slow Indicates slow change characteristics;

[0054] The rapidly changing residual X is then obtained by subtracting it from the original sequence. fast :

[0055] X fast =X VCA -X slow

[0056] To achieve effective separation of the fast and slow dynamic responses of the system.

[0057] Furthermore, step S4 specifically includes:

[0058] (1) For the slowly changing feature X slow First, it is mapped to the Koopman high-dimensional observation space by a parameter-learnable encoder:

[0059] Z slow =φ s (X slow )

[0060] Among them, Z slow A high-dimensional representation of slowly changing characteristics, φ s (·) represents the encoder function for slow features;

[0061] The globally shared learnable, slowly changing Koopman operator K s Linear progression yields predicted values ​​in high-dimensional space.

[0062]

[0063] The decoder then performs inverse mapping to obtain the prediction results for the slowly changing portion of the original data space.

[0064]

[0065] (2) For rapid change characteristics First, divide the sequence into several non-overlapping subsequences with a segment length of S and an input time step size of T.

[0066]

[0067] Where, x i This represents the i-th time window subsequence of fast features;

[0068] The segmented, rapidly changing data is mapped to the Koopman high-dimensional observation space using an observation function with the same learnable encoder structure as the slowly changing part:

[0069]

[0070] Among them, z i A high-dimensional representation of rapidly changing characteristics, φ f (·) denotes the fast feature local encoder function;

[0071] Using the extended dynamic mode decomposition (eDMD) method, the fast dynamic change Koopman operator K, which adapts to local fitting, is calculated by the least squares method. f :

[0072]

[0073] Among them, Z former This represents the high-dimensional representation of the first T / S-1 segment, Z. later This represents the high-dimensional representation from the second paragraph to the last paragraph. This represents the Moore-Penrose pseudoinverse, specifically obtained by solving for K using the least squares method. f The value;

[0074] Using K f Reconstructing and predicting known input data yields high-dimensional reconstructed values ​​and high-dimensional predicted values:

[0075]

[0076] Where T represents the number of input time steps. This represents the high-dimensional reconstruction value of the known input data. This represents the high-dimensional predicted value; then, the decoder maps the high-dimensional predicted value back to the original data space:

[0077]

[0078] (3) The residual results of each level of reconstruction are iteratively fed back to the next level to form a residual correction structure, thereby improving the modeling accuracy layer by layer. Finally, the predicted values ​​of the slow-changing and fast-changing parts of all levels are integrated as the multivariate prediction results of the model as a whole.

[0079] Furthermore, step S5 specifically includes:

[0080] (1) The physical information equation constructed based on the second-order dynamic conservation relationship of the flow rate and pressure of the water pump:

[0081]

[0082] Where Q represents volumetric flow rate, P represents pump pressure, M represents equivalent inertia coefficient, D represents system damping coefficient from pipe friction, and K represents elastic coefficient from pump body and pipe. The three coefficients are set as learnable neural network parameters.

[0083] (2) The physical information equation constructed based on the principle that the output electrical power of the motor should be approximately equal to the mechanical power consumed by the water pump:

[0084] V·I=η·Q·P

[0085] Where V represents the motor voltage, I represents the motor current, and η represents the pump's transmission efficiency.

[0086] (3) During the training of the deep predictive neural network of the VCSFA-Koopa prediction model based on slow feature analysis of variable correlation with physical information constraints, the encoder φ s Decoder φ f and the slow-changing Koppman operator K s These are learnable parameters; the total loss function for model training consists of the mean squared error (MSE) loss function, the second-order dynamic physical information loss function based on pump flow-pressure, and the physical information loss function based on the conservation of the motor's output power.

[0087]

[0088] in, This represents the prediction loss term, which is the mean squared error between the model's predicted value and the actual value. This represents the second-order kinetic loss term of the pump flow rate versus pressure. α1 and α2 are the physical loss terms for the conservation of electrical power output of the motor, and are weighting coefficients.

[0089] By minimizing the total loss To optimize the learnable parameter K s , φ s and φ f :

[0090] Furthermore, in step S5, the second-order kinetic physical loss term of the pump flow rate versus pressure... The expression is:

[0091]

[0092] Q(t) represents the predicted flow rate of the model. This represents its first-order time derivative. The second time derivative is represented by P(t), which is the predicted pump pressure value of the model. T represents the time step.

[0093] Conservation of physical losses in motor output power The expression is:

[0094]

[0095] Where V(t) represents the model's predicted motor voltage, I(t) represents the model's predicted motor current, Q(t) represents the model's predicted flow rate, and P(t) represents the model's predicted pump pressure; η represents the learnable pump transmission efficiency parameter obtained from the sigmoid mapping.

[0096] Furthermore, step S6 specifically includes: using the test set under normal operating conditions, calculating the error between the predicted value and the true value through the trained VCSFA-Koopa model based on physical information, and constructing control limits for anomaly detection based on the squared prediction error.

[0097]

[0098] Among them, SPE i Let e ​​represent the squared prediction error at the i-th time step. i Let x represent the prediction residual vector at the i-th time step. i,k y represents the true value of the k-th variable at the i-th time step. i,k For the corresponding predicted value;

[0099] The SPE statistic for all samples is estimated using the following formula:

[0100]

[0101] Where m represents the mean of SPE for all samples, and v represents the variance of SPE.

[0102] Further calculation of control limit distribution parameters:

[0103] g = v / 2m

[0104] h = 2m 2 / v

[0105] Where v represents the variance of the SPE, g represents the weighting factor of the chi-square distribution, and h is the degree of freedom parameter of the weighted chi-square distribution;

[0106] SPE control limit SPE T It is used to determine the threshold of real-time prediction error, so as to achieve accurate identification and early warning of system anomalies. Its calculation method is as follows:

[0107]

[0108] in, This represents the critical value of the chi-square distribution with a confidence level of α and degrees of freedom of h.

[0109] Furthermore, step S7 specifically includes:

[0110] (1) Construct the corresponding empirical probability distribution based on the training set data and test set data under normal working conditions, and calculate the KL divergence between training and test. As a criterion for judgment:

[0111]

[0112] Where δ represents the reference threshold, used to measure whether there is a significant distribution drift under the new operating condition; P train (i) represents the empirical probability distribution corresponding to the training set data, P test (i) represents the empirical probability distribution corresponding to the test set data;

[0113] (2) When the operating conditions change, extract test data from the new operating conditions and compare the probability distribution with the original test data under normal operating conditions to obtain the KL divergence of the actual operating condition deviation.

[0114]

[0115] in, P represents the actual KL divergence value between the new operating condition and the normal operating condition. new (i) represents the empirical probability distribution corresponding to the new working condition data, P normal (i) represents the empirical probability distribution corresponding to the original normal operating condition data;

[0116] (3) When When the value is greater than δ, it indicates that the system distribution has changed, and the SPE control limit of the prediction error needs to be recalculated based on the new operating condition data according to step S6; if If the value is less than δ, it is considered that the operating condition has not drifted significantly, and the original control limit is maintained to reduce the risk of false alarms.

[0117] Furthermore, step S8 specifically includes: the online model detection stage, which involves receiving input data for new operating conditions, obtaining predicted values ​​through a pre-trained VCSFA-Koopa model based on physical information constraints, calculating the prediction error with the actual values, and thus obtaining the real-time SPE value.

[0118]

[0119] Real-time SPE value SPE new With SPE control limit SPE T Compare, if SPE new >SPE T If it is, it is judged as abnormal; if SPE new <SPE T If it is normal, it will be judged as normal, realizing abnormal alarm and dynamic early warning.

[0120] The beneficial effects of this invention are as follows: The proposed anomaly detection method for non-stationary industrial processes, combining physical information constraints, variable correlation attention mechanisms, slow feature analysis, and Koopman high-dimensional linear modeling, addresses challenges such as fast-slow dynamic mixing, multi-condition distribution drift, and insufficient physical consistency in multivariate non-stationary time series in industry. It constructs a unified modeling framework with multi-scale decoupling, physical perception, and high-dimensional linear prediction capabilities. This method guides attention weight generation through frequency domain cross-correlation analysis, enhancing the ability to model temporal coupling between variables. Furthermore, it utilizes slow feature analysis to dynamically decouple fast and slow data in industrial time series, extracting stable evolution trends and local disturbance features. Then, it leverages Koopman operator theory to achieve linear evolution modeling of dynamics at different scales in high-dimensional space, effectively improving the model's interpretability and prediction accuracy. In addition, during model training, a physical residual loss function based on the second-order flow-pressure dynamic equation of a pump system and the power conservation formula of a motor is introduced, effectively embedding prior information about the system's structural response and energy transfer mechanism, enhancing the model's physical consistency and generalization ability. By combining the distribution drift adaptive detection mechanism of KL divergence, the anomaly detection threshold can be dynamically adjusted to achieve robust fault identification and early warning under changing operating conditions. This method is applicable to anomaly monitoring of multivariable industrial systems under complex operating conditions, providing a theoretical basis and technical support for real-time fault detection and intelligent operation and maintenance in intelligent manufacturing scenarios.

[0121] This method ensures both detection accuracy and physical interpretability, making it valuable for practical deployment and suitable for key application scenarios such as intelligent manufacturing, industrial control, and predictive maintenance of equipment.

[0122] 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

[0123] 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:

[0124] Figure 1 Flowchart of the Koopman high-dimensional space nonstationary industrial process time series anomaly detection method with embedded physical information constraints provided by the present invention;

[0125] Figure 2 The overall framework diagram of the VCSFA-Koopa prediction model with embedded physical information constraints provided by this invention;

[0126] Figure 3 The structural diagram of the variable correlation calculation module and the slow feature decomposition module of the VCSFA-Koopa prediction model with embedded physical information constraints provided by the present invention;

[0127] Figure 4 The structure diagram of the fast and slow dual-channel Koopman predictor for the VCSFA-Koopa prediction model with embedded physical information constraints provided by the present invention;

[0128] Figure 5 The VCSFA-Koopa prediction model with embedded physical information constraints provided by this invention predicts the results of multiple variables under normal operating conditions in the SKAB dataset; where (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 motor current, (d) is a comparison between the actual and predicted pressure, (e) is a comparison between the actual and predicted temperature, (f) is a comparison between the actual and predicted thermocouple values, (g) is a comparison between the actual and predicted motor voltage, and (h) is a comparison between the actual and predicted RMS values ​​of volumetric flow rate.

[0129] Figure 6 The following are some abnormal detection results under varying operating conditions in embodiments of the present invention, wherein (a) is the step behavior of closing inlet valve 1, (b) is the δ function behavior of closing inlet valve 1, (c) is the step behavior of closing inlet valve 2, (d) is the δ function behavior of closing inlet 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

[0130] 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.

[0131] Please see Figures 1-6 This invention provides a method for detecting anomalies in non-stationary industrial processes that combines physical information with Koopman spatial transform, specifically including the following steps:

[0132] Step 1: Collect multivariate time series data from the SKAB industrial process operation and classify it into normal operating condition data and various variable operating condition data. In this implementation scheme, the normal operating condition data contains 9405 time points, which are divided into three subsets in a 7:2:1 ratio: 6584 data points are used for model training, 1881 for testing, and 940 for validation. For the abnormal operating condition part, the test data lengths vary depending on the type of fault, and the abnormal state usually occurs after a period of stable operation, after which the system gradually recovers to normal operation.

[0133] The SKAB dataset (Skoltech Anomaly Benchmark) is a standardized dataset widely used in research on time-series anomaly detection in industrial processes. Its main purpose is to identify potential faults and abnormal states during the operation of industrial equipment. This dataset originates from an industrial experimental platform of a water circulation system, covering operational data under normal conditions and 34 typical abnormal conditions. The data sampling frequency is 1Hz, recording a total of 8 variables, covering key indicators such as vibration acceleration, motor voltage, motor current, water pump circuit pressure, engine temperature, circuit fluid temperature, and fluid velocity. Abnormal operating conditions encompass various types, such as inlet valve 1 closing, outlet valve 2 closing, abrupt, linear, step, delta, and exponential responses caused by rotor imbalance, and various disturbances such as changes in water tank temperature (increase or decrease), comprehensively simulating 34 typical abnormal system behaviors.

[0134] The partitioned multivariate time series data undergoes normalization preprocessing. Specifically, firstly, min-max normalization is performed using the maximum and minimum values ​​of the normal operating condition data.

[0135]

[0136] Where, x k,j This is the original value of the j-th variable at the k-th time; j = 1, 2, ..., m, and in this embodiment, j can be as high as 8. These are the maximum and minimum values ​​of variable j in the training set under normal operating conditions, respectively. This is the normalized value.

[0137] Step 2: After input data preprocessing, a variable correlation attention mechanism is introduced. By performing a Fast Fourier Transform on the multivariate time series data, the frequency domain cross-correlation coefficient between variable pairs is calculated, generating an attention matrix to reflect the correlation strength between variables, thereby enhancing the temporal interaction expression of important variables. Specifically, firstly, the normalized sequence... Perform inverted embedding, and then map it to a latent space using a multilayer perceptron:

[0138]

[0139] Where J represents the number of variables and D represents the potential spatial dimension of the mapping.

[0140] The obtained feature matrix E is reused in the representation of the query matrix, key matrix, and value matrix. Its hysteresis cross-correlation in the frequency domain is calculated:

[0141]

[0142] in, denoted by Fourier transform, * denotes complex conjugate, and ⊙ denotes Hadamard product.

[0143] The inverse Fourier transform is applied to recover the time-domain cross-correlation matrix:

[0144]

[0145] in, This represents the inverse Fourier transform.

[0146] Weighted summation of correlations over all lag steps τ:

[0147]

[0148] Where, λ τ The weighting coefficients for the learned lag step size.

[0149] The final attention weight matrix is ​​obtained by Softmax normalization:

[0150] A i,j =Softmax(COR(q) i ,k j))

[0151] The attention weight matrix is ​​used to fuse the value matrix features to obtain an enhanced representation, which is then used as input for the subsequent slow feature analysis module.

[0152] X VCA (t)=A i,j ·V

[0153] Step 3: Apply slow feature analysis to the attention-enhanced representation to decouple the multivariate time series data into subspaces with different response speeds, and extract slow features representing the long-term trend of the system and fast features characterizing rapid perturbations respectively; specifically, first, the variable-enhanced representation output by the variable correlation attention mechanism module is analyzed. By zero-meaning along the time axis, we obtain:

[0154] X c =X VCA -μ t (X VCA )

[0155] Calculate the sample covariance matrix Σ x And perform feature decomposition on it:

[0156]

[0157] Construct the whitening matrix and obtain the whitened data:

[0158]

[0159] Calculate the time derivative difference matrix:

[0160]

[0161] And define the derivative covariance matrix C based on the time derivative difference matrix. d And perform feature decomposition:

[0162]

[0163] Among them, Λ d =diag(λ1,…,λ) J The eigenvalues ​​are arranged in ascending order. The slow feature dimension k is determined based on the cumulative variance contribution rate, and the first k eigenvectors are used to construct the slow feature projection matrix W. s .

[0164] In this example, the cumulative variance contribution rate is set to 30%, the dimension k of the slow feature is determined to be 1, and the slow feature is calculated based on this projection:

[0165]

[0166] Then, the slow features are restored from the whitening space to the original input space using the inverse mapping, so as to perform fast and slow separation consistent with the original sequence:

[0167]

[0168] The rapid change residual is then obtained by subtracting it from the original sequence:

[0169] X fast =X VCA -X slow

[0170] To achieve effective separation of the fast and slow dynamic responses of the system.

[0171] Step 4: Using Koopman operator theory, fast and slow features are mapped to a high-dimensional observable linear space. Slow features are modeled using the slow-changing Koopman operator, while fast features are modeled piecewise by extracting the fast-changing Koopman operator using the extended dynamic mode decomposition (eDMD) method. A stacked residual structure is then used to iteratively optimize the model performance layer by layer, ultimately outputting the entire multivariate prediction result. Specifically, for the slow-changing feature X... slow First, it is mapped to the Koopman high-dimensional observation space by a parameter-learnable encoder:

[0172] Z slow =φ s (X slow )

[0173] The globally shared learnable, slowly changing Koopman operator K s Linear progression yields predicted values ​​in high-dimensional space:

[0174]

[0175] The decoder then performs inverse mapping to obtain the prediction results for the slowly changing portion of the original data space:

[0176]

[0177] For rapidly changing characteristics First, divide the sequence into several non-overlapping subsequences with a segment length of S and an input time step size of T. In this example, T is 24 and S is 12.

[0178]

[0179] The segmented, rapidly changing data is mapped to the Koopman high-dimensional observation space using an observation function with the same learnable encoder structure as the slowly changing part:

[0180] z i =φ f (x i ), i = 1, 2

[0181] Using the extended dynamic mode decomposition (eDMD) method, the fast dynamic change Koopman operator K, which adapts to local fitting, is calculated by the least squares method. f :

[0182] Z former =[z1]

[0183] Z later =[z2]

[0184]

[0185] Among them, Z former This represents the high-dimensional representation of the first T / S-1 segment, Z. later This represents the high-dimensional representation from the second paragraph to the last paragraph. This represents the Moore-Penrose pseudoinverse, specifically obtained by solving for K using the least squares method. f The value of .

[0186] Using K f Reconstructing and predicting known input data yields high-dimensional reconstructed values ​​and high-dimensional predicted values:

[0187]

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

[0189] The decoder then maps the predicted values ​​from the high-dimensional space back to the original data space.

[0190]

[0191] Finally, the reconstruction residuals at each level are iteratively fed back to the next level to form a residual correction structure, thereby improving the modeling accuracy layer by layer. In the end, the predicted values ​​of the slow-changing and fast-changing parts at all levels are integrated as the overall multivariate prediction results of the model.

[0192] Step 5: During model training, a physical information equation is introduced as a supervision term to construct two types of physical residual losses: one based on the second-order flow-pressure dynamic relationship of the pump system; the other based on the energy conservation equation between motor voltage, current, and output power, which is introduced into the total loss function to enhance the physical consistency and interpretability of the predictions. Specifically, the physical information equation is first constructed based on the second-order flow-pressure dynamic conservation relationship of the pump:

[0193]

[0194] Where Q represents volumetric flow rate, P represents pump pressure, M represents equivalent inertia coefficient, D represents system damping coefficient from pipe friction, and K represents elastic coefficient from pump body and pipe. These three coefficients are set as learnable neural network parameters.

[0195] Secondly, based on the physical information equation that the output electrical power of the motor should be approximately equal to the mechanical power consumed by the water pump:

[0196] V·I=η·Q·P

[0197] Where V represents the motor voltage, I represents the motor current, and η represents the pump's transmission efficiency.

[0198] In the training process of the deep predictive neural network of the VCSFA-Koopa model based on physical information constraints, the encoder φ s Decoder φ f and the slow-changing Koppman operator K s These are learnable parameters; the total loss function for model training consists of the mean squared error (MSE) loss function, the second-order dynamic physical information loss function based on pump flow-pressure, and the physical information loss function based on the conservation of the motor's output power.

[0199]

[0200] in, This represents the prediction loss term, which is the mean squared error between the model's predicted values ​​and the actual values. This refers to the physical loss term in pump dynamics. α1 and α2 are the physical loss terms for motor power conservation, and α1 and α2 are weighting coefficients.

[0201] In this example, α1 is set to 0.10 and α2 to 0.02. The goal is to minimize the total loss. To optimize the learnable parameter K s , φ s and φ f :

[0202]

[0203] The experimental setup for training the VCSFA-Koopa deep prediction neural network based on physical information constraints is as follows: VCSFA-Koopa prediction block layer number = 5, activation function is ReLU, learning rate = 0.001, optimizer is Adam, batch size = 32, epochs = 100, input sequence length = 24, prediction length = 6.

[0204] Step 6: In the offline modeling stage, the trained prediction model is evaluated for error using test data under normal operating conditions. The prediction results obtained from the normal operating condition data are as follows: Figure 5 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 motor current, (d) compares the actual and predicted pressure, (e) compares the actual and predicted temperature, (f) compares the actual and predicted thermocouple values, (g) compares the actual and predicted motor voltage, and (h) compares the actual and predicted RMS values ​​of volumetric flow rate. The squared prediction error (SPE) index is calculated based on the difference between the predicted and actual values, and a chi-square distribution anomaly detection control limit is established based on its statistical characteristics. Using test data under normal operating conditions, the error between the predicted and actual values ​​is calculated through a trained VCSFA-Koopa model based on physical information constraints, and a control limit for anomaly detection is constructed based on the squared prediction error.

[0205]

[0206] Where, x i,k y represents the true value of the k-th variable at the i-th time step. i,k For the corresponding predicted value.

[0207] The SPE statistic for all samples is estimated using the following formula:

[0208]

[0209] Where m represents the mean of the SPE for all samples, and v represents the variance of the SPE. Further calculation of the control limit distribution parameters is then performed.

[0210] g = v / 2m

[0211] h = 2m 2 / v

[0212] Where v represents the variance of the SPE, and h is the degree of freedom parameter of the weighted chi-square distribution.

[0213] The final SPE control limits for anomaly detection are:

[0214]

[0215] in, This represents the critical value of the chi-square distribution with a confidence level of α and degrees of freedom of h. Specifically, in this embodiment, α is set to 0.05. The SPE control limit is used to perform threshold discrimination on real-time prediction errors, thereby achieving accurate identification and early warning of system anomalies.

[0216] Step 7: To address the distribution drift issue caused by the evolution of operating conditions over time, KL divergence is used to measure the distribution difference between the current operating condition and the original normal operating condition. Based on this, it is determined whether to update the SPE control limits to improve the adaptability of anomaly detection in changing operating conditions. Specifically, firstly, empirical probability distributions are constructed based on the training and test sets of normal operating conditions, and the KL divergence between training and testing is calculated as the discrimination criterion.

[0217]

[0218] Where δ represents the reference threshold, used to measure whether there is a significant distribution drift under the new operating condition; P train (i) represents the empirical probability distribution corresponding to the training set data, P test (i) represents the empirical probability distribution corresponding to the test set data.

[0219] When the operating conditions change, test data is extracted from the new operating conditions and compared with the original test data under normal operating conditions using a probability distribution to obtain the KL divergence of the actual operating condition deviation.

[0220]

[0221] in, P represents the actual KL divergence value between the new operating condition and the normal operating condition. new (i) represents the empirical probability distribution corresponding to the new working condition data, P normal (i) represents the empirical probability distribution corresponding to the original normal operating condition data.

[0222] Step 8: In the online detection phase, real-time collected data is input into the model to obtain predicted values ​​and calculate the SPE value. This value is then compared with the current control limits to determine if the current state is abnormal, achieving rapid fault detection and dynamic early warning. Specifically, input data from new operating conditions is received, and a pre-trained VCSFA-Koopa model based on physical information is used to obtain predicted values. The prediction error is then calculated against the actual values ​​to obtain the real-time SPE value.

[0223]

[0224] Compare the real-time SPE value with the SPE control limit. If the SPE... new >SPE T If it is, it is judged as abnormal; if SPE new <SPE T If it is normal, it will be judged as normal, realizing abnormal alarm and dynamic early warning. Figure 6The following are some abnormal detection results under varying operating conditions in this embodiment: (a) is the step behavior of closing inlet valve 1, (b) is the δ function behavior of closing inlet valve 1, (c) is the step behavior of closing inlet valve 2, (d) is the δ function behavior of closing inlet 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.

[0225] 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 with embedded physical information constraints, characterized in that, The method specifically includes the following steps: S1: Acquire multivariate time series data generated during operation in the industrial system, divide it into normal operating condition data and abnormal operating condition data containing faults or changing factors; divide the normal operating condition data into training set, validation set and test set according to the proportion, and perform normalization preprocessing on the original multivariate time series data; S2: After input data preprocessing, a variable correlation attention mechanism is introduced. By performing a fast Fourier transform on multivariate time series data, the frequency domain cross-correlation coefficient between variables is calculated, and an attention weight matrix is ​​generated to reflect the correlation strength between variables, thereby enhancing the temporal interaction expression between important variables. S3: Apply the slow feature analysis method to the attention-enhanced representation to decouple the multivariate time series data into subspaces with different response speeds, and extract the slow features that represent the long-term trend of the system and the fast features that characterize rapid disturbances respectively. S4: Using the Koopman operator theory, fast and slow features are mapped to a high-dimensional observable linear space respectively. Slow features are modeled by the slow-changing Koopman operator, and fast features are modeled by extracting the fast-changing Koopman operator through the extended dynamic mode decomposition method. The model performance is then iteratively optimized layer by layer using a stacked residual structure, and finally the whole multivariate prediction results are output. S5: During the training of the prediction model, a physical information equation is introduced as a supervision term to construct two types of physical residual losses: one is based on the second-order flow-pressure dynamic relationship of the pump system; the other is based on the energy conservation equation between motor voltage, current and output power, which is introduced into the total loss function to enhance the physical consistency and interpretability of the prediction. S6: In the offline modeling stage, the training prediction model is evaluated using the test set under normal working conditions. The squared prediction error (SPE) index is calculated based on the difference between the predicted value and the actual value, and an anomaly detection control limit based on its statistical characteristics is established using a weighted chi-square distribution. S7: To address the distribution drift problem under changing operating conditions, KL divergence is used to measure the distribution difference between the current operating condition and the original normal operating condition, and based on this, it is determined whether to update the SPE control limit to improve the adaptability of anomaly detection in changing operating conditions. S8: During the online detection phase, the real-time collected data is input into the model to obtain the predicted value and calculate the SPE value. The value is then compared with the current control limit to determine whether the current state is abnormal, thus achieving rapid fault detection and dynamic early warning.

2. The method for detecting anomalies in non-stationary industrial processes according to claim 1, characterized in that, In step S1, the collected multivariate time series data includes time series variables from multiple operating conditions. The normal operating condition data is divided into training, validation, and test sets in a 7:2:1 ratio. Let the multivariate time series data have J measurement variables and T time steps, forming a T×J matrix. The normal operating condition data and the abnormal operating condition data are denoted as X... n (N n ×J) and X a (N a ×J), where N n and N a These represent the corresponding sample numbers; for dataset X n and X a Perform the following maximum and minimum value normalization: Where, x k,j This represents the original value of the j-th variable at time k. These are the maximum and minimum values ​​of variable j in the training set, respectively; This is the normalized value.

3. The method for detecting anomalies in non-stationary industrial processes according to claim 2, characterized in that, In step S2, a variable-related attention mechanism is introduced, specifically including: for the normalized sequence Perform inverted embedding, and then map it to a latent space using a multilayer perceptron: Where E represents the feature matrix, T represents the time step, J represents the number of variables, D represents the latent space dimension of the mapping, Permute represents dimension transpose, and MLP represents multilayer perceptron; The obtained feature matrix E is reused to represent the query matrix Q, key matrix K, and value matrix V; its hysteresis cross-correlation in the frequency domain is calculated: Among them, S QK (f) represents the correlation in the Q and K frequency domains. denoted by Fourier transform, * denotes complex conjugate, and ⊙ denotes Hadamard product; The inverse Fourier transform is applied to recover the time-domain cross-correlation matrix: Among them, R Q,K (τ) represents the correlation between Q and K in the time domain. Indicates the inverse Fourier transform; Weighted summation of correlations over all lag steps τ: Where τ represents the lag time step, COR(q) i ,k j ) represents the variable q i ,k j Total correlation across all lag steps Represents variable q i ,k j The correlation at lag step τ, λ τ The weighting coefficients for the learned lag step size; The final attention weight matrix A is obtained by Softmax normalization. i,j : A i,j =Softmax(COR(q i ,k j )) The value matrix features are fused using an attention weight matrix to obtain an enhanced representation, which is then used as input for the subsequent slow feature analysis module. X VCA (t)=A i,j ·V Among them, X VCA (t) represents the enhanced transvariate lag correlation representation based on attention mechanism weighting, and V represents the value matrix.

4. The method for detecting anomalies in non-stationary industrial processes according to claim 3, characterized in that, In step S3, the slow feature analysis method specifically includes: augmenting the variable representation output by the variable correlation attention mechanism module. By zero-meaning along the time axis, we obtain: X c =X VCA -m t (X VCA ) Among them, X c This represents the augmented representation of the centered variable, μ. t (·) indicates the mean along the time axis; Calculate the sample covariance matrix Σ x And perform feature decomposition on it: Among them, V x Represents the covariance matrix Σ x eigenvector matrix, Λ x Represents the covariance matrix Σ x The eigenvalue diagonal matrix; Construct the whitening matrix and obtain the whitened data: Among them, W w Represents the whitening transformation matrix. This represents the whitened input. Next, calculate the time derivative difference matrix: Where D represents the time derivative difference matrix; And define the derivative covariance matrix C based on the time derivative difference matrix. d And perform feature decomposition: Among them, Λ d =diag(λ1,…,λ) J ), the eigenvalues ​​are arranged in ascending order; V d Represents the time derivative covariance matrix; The slow feature dimension k is determined based on the cumulative variance contribution rate, and the first k feature vectors are used to construct the slow feature projection matrix W. s Based on this, the slow feature Y is calculated by projection. slow : Then, the slow features are restored from the whitening space to the original input space using the inverse mapping, so as to perform fast and slow separation consistent with the original sequence: Among them, W o W represents the mapping matrix projected from the slow feature space back to the original input space. w Let X represent the whitening transformation matrix. slow Indicates slow change characteristics; The rapidly changing residual X is then obtained by subtracting it from the original sequence. fast : X fast =X VCA -X slow To achieve effective separation of the fast and slow dynamic responses of the system.

5. The method for detecting anomalies in non-stationary industrial processes according to claim 4, characterized in that, Step S4 specifically includes: (1) For the slowly changing feature X slow First, it is mapped to the Koopman high-dimensional observation space by a parameter-learnable encoder: Z slow =φ s (X slow ) Among them, Z slow A high-dimensional representation of slowly changing characteristics, φ s (·) represents the encoder function for slow features; The globally shared learnable, slowly changing Koopman operator K s Linear progression yields predicted values ​​in high-dimensional space. The decoder then performs inverse mapping to obtain the prediction results for the slowly changing portion of the original data space. (2) For rapid change characteristics First, divide the sequence into several non-overlapping subsequences with a segment length of S and an input time step size of T. Where, x i This represents the i-th time window subsequence of fast features; The segmented, rapidly changing data is mapped to the Koopman high-dimensional observation space using an observation function with the same learnable encoder structure as the slowly changing part: Among them, z i A high-dimensional representation of rapidly changing characteristics, φ f (·) denotes the fast feature local encoder function; Using the extended dynamic mode decomposition method, the fast dynamic change Koopman operator K, which adapts to local fitting, is calculated by the least squares method. f : Among them, Z former This represents the high-dimensional representation of the first T / S-1 segment, Z. later This represents the high-dimensional representation from the second paragraph to the last paragraph. This represents the Moore-Penrose pseudoinverse, specifically obtained by solving for K using the least squares method. f The value; Using K f 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. Represents the predicted value in a higher dimension; The decoder then maps the predicted values ​​from the high-dimensional space back to the original data space. (3) The residual results of each level of reconstruction are iteratively fed back to the next level to form a residual correction structure, thereby improving the modeling accuracy layer by layer. Finally, the predicted values ​​of the slow-changing and fast-changing parts of all levels are integrated as the multivariate prediction results of the model as a whole.

6. The method for detecting anomalies in non-stationary industrial processes according to claim 5, characterized in that, Step S5 specifically includes: (1) The physical information equation constructed based on the second-order dynamic conservation relationship of the flow rate and pressure of the water pump: Where Q represents volumetric flow rate, P represents pump pressure, M represents equivalent inertia coefficient, D represents system damping coefficient from pipe friction, and K represents elastic coefficient from pump body and pipe. The three coefficients are set as learnable neural network parameters. (2) The physical information equation constructed based on the principle that the output electrical power of the motor should be approximately equal to the mechanical power consumed by the water pump: V·I=η·Q·P Where V represents the motor voltage, I represents the motor current, and η represents the pump's transmission efficiency; (3) During the training process of the deep predictive neural network of the Koopman prediction model based on slow feature analysis of variable correlations constrained by physical information, the encoder φ s Decoder φ f and the slow-changing Koppman operator K s These are learnable parameters; the total loss function for model training consists of the mean squared error (MSE) loss function, the second-order dynamic physical information loss function based on pump flow-pressure, and the physical information loss function based on the conservation of the motor's output power. in, This represents the prediction loss term, which is the mean squared error between the model's predicted value and the actual value. This represents the second-order kinetic loss term of the pump flow rate versus pressure. α1 and α2 are the physical loss terms for the conservation of electrical power output of the motor, and are weighting coefficients. By minimizing the total loss To optimize the learnable parameter K s , φ s and φ f :

7. The method for detecting anomalies in non-stationary industrial processes according to claim 6, characterized in that, In step S5, the second-order kinetic physical loss term of the pump flow rate versus pressure is... The expression is: Q(t) represents the predicted flow rate of the model. This represents its first-order time derivative. The second time derivative is represented by P(t), which is the predicted pump pressure value of the model. T represents the time step. Conservation of physical losses in motor output power The expression is: Where V(t) represents the model's predicted motor voltage, I(t) represents the model's predicted motor current, Q(t) represents the model's predicted flow rate, and P(t) represents the model's predicted pump pressure; η represents the learnable pump transmission efficiency parameter obtained from the sigmoid mapping.

8. The method for detecting anomalies in non-stationary industrial processes according to claim 6, characterized in that, Step S6 specifically includes: using the test set under normal operating conditions, calculating the error between the predicted value and the true value through the trained VCSFA-Koopa model based on physical information, and constructing control limits for anomaly detection based on the squared prediction error. Among them, SPE i Let e ​​represent the squared prediction error at the i-th time step. i Let x represent the prediction residual vector at the i-th time step. i,k y represents the true value of the k-th variable at the i-th time step. i,k For the corresponding predicted value; The SPE statistic for all samples is estimated using the following formula: Where m represents the mean of SPE for all samples, and v represents the variance of SPE. Further calculation of control limit distribution parameters: g = v / 2m h=2m 2 / v Where v represents the variance of the SPE, g represents the weighting factor of the chi-square distribution, and h is the degree of freedom parameter of the weighted chi-square distribution; SPE control limit SPE T It is used to determine the threshold of real-time prediction error, so as to achieve accurate identification and early warning of system anomalies. Its calculation method is as follows: in, This represents the critical value of the chi-square distribution with a confidence level of α and degrees of freedom of h.

9. The method for detecting anomalies in non-stationary industrial processes according to claim 8, characterized in that, Step S7 specifically includes: (1) Construct the corresponding empirical probability distribution based on the training set data and test set data under normal working conditions, and calculate the KL divergence between training and test. As a criterion for judgment: Where δ represents the reference threshold, used to measure whether there is a significant distribution drift under the new operating condition; P train (i) represents the empirical probability distribution corresponding to the training set data, P test (i) represents the empirical probability distribution corresponding to the test set data; (2) When the operating conditions change, extract test data from the new operating conditions and compare the probability distribution with the original test data under normal operating conditions to obtain the KL divergence of the actual operating condition deviation. in, P represents the actual KL divergence value between the new operating condition and the normal operating condition. new (i) represents the empirical probability distribution corresponding to the new working condition data, P normal (i) represents the empirical probability distribution corresponding to the original normal operating condition data; (3) When When the value is greater than δ, it indicates that the system distribution has changed, and the SPE control limit of the prediction error needs to be recalculated based on the new operating condition data according to step S6; if If the value is less than δ, it is considered that the operating condition has not drifted significantly, and the original control limit is maintained to reduce the risk of false alarms.

10. The method for detecting anomalies in non-stationary industrial processes according to claim 9, characterized in that, Step S8 specifically includes: the online model detection stage, which receives input data for new working conditions, obtains predicted values ​​through a trained VCSFA-Koopa model based on physical information constraints, calculates the prediction error with the true values, and thus obtains the real-time SPE value. Real-time SPE value SPE new With SPE control limit SPE T Compare, if SPE new SPE T If it is, it is judged as abnormal; if SPE new <SPE T If it is normal, it will be judged as normal, realizing abnormal alarm and dynamic early warning.

Images