An anomaly detection method for water treatment systems based on adaptive mapping and frequency domain decoupling

Abstract

This invention discloses an anomaly detection method for water treatment systems based on adaptive graph construction and frequency domain decoupling, comprising the following steps: acquiring and processing multivariate time-series data of the water treatment system; constructing a variable correlation graph, calculating the dependencies between variables through mutual information, and adaptively determining the connection relationships between nodes using density clustering; using a graph convolutional network to perform feature propagation and aggregation on the variable correlation graph to obtain a spatial feature representation; mapping the spatial features to the frequency domain, and then performing frequency domain decoupling learning based on spectral energy to obtain frequency domain enhanced features; fusing the spatial features and frequency domain enhanced features through residual fusion to obtain fused features; using a Transformer encoder to perform temporal modeling on the fused features to obtain a feature representation containing global temporal dependencies, and reconstructing the input sequence using a reconstruction module; training the model using historical normal data; and finally using the trained model to detect anomalies in the time-series data of the water treatment system to be analyzed.

Description

Technical Field

[0001] This invention relates to the field of anomaly detection technology in water treatment systems, and specifically to an anomaly detection method for water treatment systems based on adaptive mapping and frequency domain decoupling. Background Technology

[0002] With the development of industrial automation and industrial internet technologies, industrial control systems (ICS) are widely used in water treatment systems, waterworks, wastewater treatment plants, and industrial production processes. Water treatment systems typically deploy a large number of sensors and actuators to collect real-time system operating data, such as water level, flow rate, pressure, temperature, and chemical indicators, thereby forming high-dimensional multivariate time series (MTS) data. By analyzing this data, the system's operating status can be monitored in real time, and equipment failures, system anomalies, or network attacks can be detected promptly.

[0003] In recent years, with the increasing prominence of industrial cybersecurity issues, anomaly detection technologies for industrial control systems have received widespread attention. The iTrust Lab at the Singapore University of Technology and Design (SUTD) has built a real-world water treatment system experimental platform and released the SWaT (Secure Water Treatment) dataset for researching anomaly detection in industrial control systems. These datasets record sensor and actuator data from various stages of the water treatment system and are widely used in industrial anomaly detection research.

[0004] Currently, multivariate time series anomaly detection methods can be mainly divided into three categories: statistical methods, machine learning methods, and deep learning methods. However, most methods still have the following problems in industrial multivariate time series anomaly detection:

[0005] 1. Insufficient ability to model variable relationships;

[0006] In water treatment systems, complex coupling relationships often exist between various sensor variables. For example, there is a strong correlation between flow rate and pressure. However, traditional methods often assume that the variables are independent, making it difficult to effectively capture the structural relationships between them. Traditional correlation analysis methods (such as Pearson correlation coefficient) can only describe linear relationships, while nonlinear dependencies often exist between variables in industrial systems, leading to inaccurate relationship modeling.

[0007] 2. Insufficient multi-scale temporal pattern coupling and frequency feature modeling;

[0008] Monitoring data from industrial water treatment systems typically exhibits significant multi-scale characteristics, meaning that time-series signals simultaneously contain low-frequency trend changes and high-frequency dynamic fluctuations. For example, long-term changes in system operating status usually manifest as low-frequency trends, while equipment control actions or local anomalies may exhibit high-frequency fluctuations, and different frequencies often show different variable interactions.

[0009] Many existing anomaly detection methods primarily model time-series data in the temporal domain, such as using recurrent neural networks (LSTM, GRU) or Transformers to learn temporal dependencies. However, these methods typically model all frequency components together, making it difficult to effectively distinguish dynamic patterns at different frequency scales, resulting in insufficient expressive power for multi-scale time-series features. Furthermore, since different frequency components in industrial time series often have different physical meanings, without explicit modeling of the frequency structure, the model may struggle to accurately identify anomalous patterns when faced with complex dynamic changes.

[0010] 3. Insufficient ability to model long-term dependencies;

[0011] Time-series data generated during the operation of industrial equipment typically exhibit significant long-term time-series dependencies. The operating state of the equipment is often related not only to the current observation but also to historical state changes over a longer period. For example, in complex industrial systems, some abnormal states are often caused by long-term slow changes or multi-stage cumulative effects. Summary of the Invention

[0012] This invention aims to solve the technical problems existing in the prior art, and innovatively proposes a water treatment system anomaly detection method based on adaptive mapping and frequency domain decoupling, so as to analyze the detection data and equipment status information during the operation of the water treatment system in real time and improve the accuracy of anomaly monitoring during the operation of the water treatment system.

[0013] To achieve the above-mentioned objectives of this invention, this invention provides a method for anomaly detection in a water treatment system based on adaptive mapping and frequency domain decoupling, comprising the following steps:

[0014] S1: Data preparation: Collect the detection values ​​or status values ​​of various monitoring variables of the water treatment system as multivariate time series data, standardize the collected multivariate time series data, and perform sliding window sampling on the standardized multivariate time series data.

[0015] S2: Construct a time-series data anomaly detection model for a water treatment system. A graph representation learning module is used to represent the features of multivariate time-series data, obtaining spatial features reflecting the dependencies between various monitored variables. These spatial features are then transformed and enhanced in the frequency domain to obtain enhanced frequency domain features. The spatial features and enhanced frequency domain features are fused to obtain fused features. Time-series learning is performed based on these fused features to obtain a time-series feature representation, which is then reconstructed to obtain a reconstructed feature representation.

[0016] S3: Train the model constructed in step S2 using the training set to obtain the time series data anomaly detection model for the water treatment system;

[0017] S4: Based on the time series data anomaly detection model of the water treatment system obtained in step S3, determine the anomaly judgment threshold, and perform anomaly detection on the multivariate time series data to be analyzed according to the anomaly judgment threshold.

[0018] In the above scheme, step S1 further includes the following steps:

[0019] S1-1: Standardize each monitoring variable in the multivariate time series data along the time dimension;

[0020] Multivariate time series data are represented as follows: ;in, Representing time series Total length, Indicates in Constantly monitor the water treatment system A vector consisting of the observed values ​​of 1 monitored variable, where each component corresponds to an observed value of one monitored variable; and ,in, Indicates the first in the water treatment system One monitoring variable in The detected value or state value at any given time. This represents the total number of monitored variables in the water treatment system;

[0021] The Z-score method is used to standardize time series from different monitoring variables. The monitored variable is at the _ The standardized formula corresponding to time is:

[0022] ;

[0023] in, Indicates the first The monitored variable is at the _ The detected value or state value at any given time; Indicates the first The mean of the time series of each monitored variable is calculated using the following formula: ; Indicates the first The standard deviation of the time series of each monitored variable is calculated using the following formula: , Indicates the first Time series data of 1 monitored variable after Z-score standardization;

[0024] S1-2: Perform sliding window sampling on the standardized multivariate time series data;

[0025] Multivariate time-series data are sampled using a sliding window approach to construct time window samples for model training and anomaly detection; at time... At this location, the construction length is Time window sample: ;in, This indicates the preset window length. Indicated by The time window sample is the end time. Indicates the water treatment system in The observation vector at time t, Indicates the water treatment system in The observation vector at time t.

[0026] In the above scheme, step S2 further includes the following steps:

[0027] S2-1: The graph representation learning module is used to represent the features of multivariate time series data, including constructing a graph structure that reflects the dependencies between the monitoring variables, and performing feature propagation and aggregation based on the graph structure to obtain spatial features that reflect the dependencies between the monitoring variables.

[0028] S2-2: Perform frequency domain transformation on spatial features to map them to the frequency domain space to obtain frequency domain features, and perform frequency domain decoupling learning on the frequency domain features to obtain frequency domain enhanced features; fuse spatial features and frequency domain enhanced features to obtain features that fuse spatial structure information and frequency domain enhanced information;

[0029] S2-3: The features that integrate spatial structure information and frequency domain enhancement information are positionally encoded, and then input into the Transformer encoder for temporal feature learning to obtain a feature representation containing the temporal dependencies of each monitoring variable; the feature representation containing the temporal dependencies of each monitoring variable is then input into the reconstruction module, and the feature representation is reconstructed through nonlinear mapping to obtain the reconstructed feature representation.

[0030] In the above scheme: the graph structure in step S2-1 is represented as follows ,in, Represents a set of nodes. Represents the set of edges between nodes;

[0031] Step S2-1 also includes the following steps:

[0032] S2-1-1: Calculate the mutual information values ​​between the monitored variables to obtain the statistical dependencies between the variables, and construct a mutual information matrix based on the mutual information values;

[0033] The formula for calculating the mutual information value is:

[0034] ;

[0035] in, Indicates the first The monitoring variable and the first Mutual information values ​​between the monitored variables; Indicates the first The monitored variable is at the _ Feature vectors within a time window; Indicates the first Marginal entropy of each monitored variable; Indicates the first The monitored variable is at the _ Feature vectors within a time window; Indicates the first Marginal entropy of each monitored variable; Indicates the first The monitoring variable and the first The joint entropy of the monitored variables; the marginal entropy and joint entropy are calculated statistically based on the value sequence of each monitored variable within the time window through probability distribution. The constructed mutual information matrix is ​​represented as follows: ;

[0036] S2-1-2: Determine the distance matrix between monitored variables based on the mutual information matrix. Based on the distance matrix, the DBSCAN clustering algorithm is used to adaptively partition each monitoring variable to obtain a cluster structure. The monitoring variables that are not assigned to any cluster are recorded as noise points and marked as -1. The adjacency relationship between the monitoring variables is determined according to the cluster structure, and an adjacency matrix is ​​constructed.

[0037] Distance matrix The corresponding calculation formula is:

[0038] ;

[0039] in, Let represent the distance matrix, and Let be a matrix whose elements are all 1, and It is a mutual information matrix, and Its element values ​​are the mutual information values ​​between variables; Represents a matrix Perform element-wise operations, when When an element in the matrix has a negative value, it is set to 0 to ensure the non-negativity of the distance matrix.

[0040] The adjacency matrix between the constructed monitoring variables is represented as follows:

[0041] ;

[0042] in, Indicates the first The first monitoring variable and the first If all monitored variables are in the same cluster, otherwise The adjacency matrix is ​​used to characterize the connection relationships between the monitored variables and to provide graph structure information for subsequent graph convolution operations.

[0043] S2-1-3: Based on the diagram structure A graph convolutional network is used to extract spatial features that reflect the dependencies between monitoring variables, and the data of different monitoring variables are integrated into a unified vector representation.

[0044] The formula for calculating graph convolution is as follows:

[0045] ;

[0046] in, The graph convolutional network represents the first... The first time window The output feature matrix of the layer, and Indicates the embedding dimension; The graph convolutional network represents the first... The first time window The output feature matrix of layer +1; This represents the graph convolutional network. The trainable weight matrix corresponding to the layer; It is a non-linear activation function; The extended adjacency matrix is ​​calculated using the following formula: ,in, Indicates the first The adjacency matrix corresponding to each time window Represents the identity matrix. Representation and extended adjacency matrix The corresponding degree matrix is ​​a diagonal matrix, and the formula for calculating the diagonal elements is: ,in, Represents the th in the extended adjacency matrix Line number Column elements; This represents the result after symmetric normalization of the extended adjacency matrix;

[0047] The spatial feature vectors of each monitored variable obtained through the graph convolutional network are represented as follows:

[0048] ;

[0049] in, Indicates the first The monitored variable is at the _ Spatial feature vector at time step, Indicates the embedding dimension. This indicates the number of monitored variables in the water treatment system;

[0050] All monitored variables in the first The spatial feature vectors at each time step are concatenated to obtain the spatial feature representation:

[0051] ;

[0052] in, For the first The spatial features of each time window will be represented in the subsequent model. Record ,Right now Indicates the first The spatial feature matrix of the nth time window, the th OK Indicates the first The monitored variable is at the _ Spatial feature vectors for each time window.

[0053] In the above scheme, step S2-2 also includes the following steps:

[0054] S2-2-1: Perform a Fast Fourier Transform (FFT) on the spatial features along the time dimension to obtain the corresponding frequency domain feature representation. Then, use an adaptive threshold based on spectral energy to partition the frequency domain feature representation, obtaining low-frequency and high-frequency components. For spatial features... Its change process can be represented as follows:

[0055] ;

[0056] in, Represents the Fast Fourier Transform. Represents frequency domain characteristics, and ,in The length of the frequency domain representation; This represents the total number of monitored variables in the water treatment system;

[0057] The formula for calculating spectral energy is:

[0058] ;

[0059] in, This indicates that the square operation is performed on each element of the matrix; The power spectrum is used to characterize the energy intensity distribution of each frequency component.

[0060] The low-frequency and high-frequency components can be obtained using the following formula:

[0061] ;

[0062] ;

[0063] in, This represents element-wise multiplication; This indicates an indicator function used to perform element-by-element discrimination of the power spectrum. When the condition in parentheses is true, the function at the corresponding position takes the value of 1; otherwise, it takes the value of 0. Indicates the energy adaptive threshold; This indicates that when the corresponding element in the power spectrum matrix is ​​greater than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0. This indicates that when the corresponding element in the power spectrum matrix is ​​less than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0.

[0064] S2-2-2: Model the low-frequency components and high-frequency components separately to obtain low-frequency features and high-frequency features, and then fuse the low-frequency features and high-frequency features;

[0065] The calculation formulas for modeling low-frequency and high-frequency components are as follows:

[0066] ;

[0067] ;

[0068] in, This represents the frequency domain feature modeling function constructed for low-frequency components; Indicates low-frequency characteristics; This represents the frequency domain feature modeling function constructed for high-frequency components; Indicates high-frequency characteristics;

[0069] The formula for fusing low-frequency and high-frequency features is:

[0070] ;

[0071] in, This represents the frequency domain characteristics after fusion; This represents element-wise addition. It is a low-frequency characteristic. High-frequency characteristics;

[0072] S2-2-3: The fused features are mapped back to the time domain space through inverse fast Fourier transform to obtain a frequency-domain enhanced time-series feature representation;

[0073] The formula for calculating the inverse fast Fourier transform is:

[0074] ;

[0075] in, This represents the time-series characteristics after frequency domain enhancement. This represents the inverse fast Fourier transform;

[0076] S2-2-4: Representing the enhanced time-series characteristics in the frequency domain Spatial feature representation The fusion is performed through residual connections to obtain the fused feature representation;

[0077] The calculation formula for the fusion process is as follows:

[0078] ;

[0079] in, This represents the temporal features after fusion; This represents the learnable parameters.

[0080] In the above scheme, step S2-3 also includes the following steps:

[0081] S2-3-1: Representation of fusion features Perform positional encoding and then perform layer normalization on the position-encoded features;

[0082] The formula for calculating the position code is:

[0083] ;

[0084] ;

[0085] in, Indicates the time step position index; Indicates the feature dimension index; This represents the position encoding function; This represents the sinusoidal component of the positional encoding in the even-numbered dimensions. This represents the cosine component of the positional encoding in odd-numbered dimensions; Indicates the feature embedding dimension. Represents the sine function. Represents the cosine function;

[0086] The formula for calculating layer normalization of the position-encoded features is as follows:

[0087] ;

[0088] in, The time-series feature representation after normalization of the representation layer. This represents the temporal features after fusion with step S2-2-4. A dimensionally consistent location encoding matrix is ​​used to represent the location information at each time step. This represents the fused temporal feature representation. This represents the standard layer normalization operation;

[0089] S2-3-2: Representing the time-series features after layer normalization The Transformer encoder is used to learn temporal features to capture the temporal dependencies of various monitoring variables in the water treatment system. The Transformer encoder includes a multi-head self-attention module, a residual connection and layer normalization module, and a feedforward neural network module. The output features of the Transformer encoder are then input into the reconstruction module for reconstruction.

[0090] In the above scheme, step S2-3-2 also includes the following steps:

[0091] S2-3-2-1: Representing the time-series features after layer normalization By transforming the input representations of each attention head in the multi-head attention mechanism into query vectors, key vectors, and value vectors through linear mapping, the results are obtained.

[0092] The linear mapping formula is expressed as follows:

[0093] ;

[0094] in, The first The projection weight matrix corresponding to the attention head query vector, key vector, and value vector; Indicates the first The query vector corresponding to each attention head; Indicates the first The key vector corresponding to each attention head; Indicates the first The value vector corresponding to each attention head; Indicates the feature embedding dimension; Indicates the number of attention heads;

[0095] S2-3-2-2: Based on the query vector Key vector Sum value vector Calculate the output of each attention head;

[0096] The attention calculation process is expressed by the following formula:

[0097] ;

[0098] in, Indicates the first The output of each attention head; This is the activation function used to transform the attention weights into a probability distribution; This represents the attention relevance score matrix. For the first The value vector corresponding to each attention head;

[0099] S2-3-2-3: The outputs of each attention head are concatenated and a linear transformation is used to obtain the output feature representation of the multi-head self-attention mechanism;

[0100] The calculation process is expressed by the following formula:

[0101] ;

[0102] in, This represents the output of the multi-head self-attention mechanism; This indicates a concatenation operation of the output results of each attention head; This indicates the output mapping weight matrix; Indicates the number of attention heads; This represents the output of the first attention head. Indicates the first The output of each attention head;

[0103] S2-3-2-4: Output of bullish self-attention With input feature representation Perform residual connections and layer normalization to obtain feature representations that include the temporal dependencies of each monitored variable;

[0104] The calculation process for feature representation is as follows:

[0105] ;

[0106] in, A feature representation that includes the temporal dependencies of each monitored variable; This represents the standard layer normalization operation; Characteristics representing multi-head self-attention output;

[0107] S2-3-2-5: Feature representation containing the time-series dependencies of each monitored variable By performing a feedforward neural network transformation and combining residual connections with layer normalization, the final temporal feature representation of the Transformer encoder output is obtained.

[0108] The calculation process is expressed by the following formula:

[0109] ;

[0110] in, This represents the feature representation after transformation by a feedforward neural network, residual connections, and layer normalization. This represents a feedforward neural network; This represents the standard layer normalization operation;

[0111] S2-3-2-6: Input the final temporal features output by the Transformer encoder into the reconstruction module for reconstruction to obtain the reconstructed feature representation;

[0112] The reconstruction process can be represented by the following formula:

[0113] ;

[0114] in, Indicates the first The reconstruction results for each time window; Indicates the reconstructed mapping function; This represents the final temporal feature representation obtained by the Transformer encoder in step S2-3-2-5.

[0115] In the above scheme, step S3 further includes the following steps:

[0116] S3-1: Construct a contrastive learning objective based on the similarity relationship between node pairs, and calculate the binary cross-entropy loss between the predicted similarity probability of node pairs and the true label. The binary cross-entropy loss function is expressed as follows:

[0117] ;

[0118] in, Represents the binary cross-entropy loss function; and For positive sample node index, Positive sample node pairs represent node pairs that have a real dependency relationship in the graph structure; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates matrix transpose. This represents a measure of the similarity of vector dot products. Represents the set of positive sample node pairs, which contains node pairs in the graph structure that have real dependencies; This represents the set of negative sample node pairs, obtained through negative sampling. For negative sample node index, Negative sample node pairs represent node pairs that do not have a real dependency relationship in the graph structure. Indicates the first Within the first time window The spatial feature vectors of each monitored variable corresponding to a graph node Indicates the first Within the first time window The spatial characteristic orientation of each monitored variable corresponding to a graph node; This represents the Sigmoid activation function;

[0119] S3-2: Calculate the input time characteristics With reconstruction results The mean square error between the two values ​​is used to obtain the reconstruction loss;

[0120] The formula for calculating the reconstruction loss is as follows:

[0121] ;

[0122] in, Indicates the reconstruction loss; This indicates the reconstruction result for the current time window; This is the temporal feature representation of S2-3-1 after position encoding and layer normalization; Denotes the square of the Frobenius norm;

[0123] S3-3: The following objective function is used to jointly train the model:

[0124] ;

[0125] in, Indicates the total loss; Indicates the reconstruction loss; This represents the binary cross-entropy loss; This represents a hyperparameter used to balance the reconstruction loss and the binary cross-entropy loss.

[0126] In the above scheme, step S4 also includes the following steps:

[0127] S4-1: Based on the reconstruction results With the corresponding input time series features Calculate the reconstruction error for each time window; use the following formula to calculate the reconstruction error and use it as the anomaly score. :

[0128] ;

[0129] in, Indicates the first Abnormal scores within a time window; Denotes the square of the Frobenius norm;

[0130] S4-2: Based on the validation set data, calculate the anomaly score, generate candidate thresholds within the anomaly score range, iterate through these candidate thresholds, and select the candidate threshold that optimizes the evaluation index as the anomaly determination threshold. ;

[0131] S4-3: Determine the status of the current time window based on the anomaly detection threshold;

[0132] The status of the current time window is determined using the following formula:

[0133] ;

[0134] in, Indicates the first Anomaly detection results for each time window; This indicates the abnormal score for the corresponding time window; This indicates the threshold for anomaly detection.

[0135] In summary, due to the adoption of the above technical solutions, the beneficial effects of this invention are as follows: Addressing the strong coupling relationships, nonlinear dependencies, and multi-scale dynamic changes in the time series of multiple monitored variables (such as flow rate, pressure, water level, pH value, turbidity, residual chlorine concentration, etc., sensors, pumps, valves, etc.) within a water treatment system, this invention introduces adaptive graph structure learning and frequency domain decoupling learning methods respectively. These methods separate and model different frequency components, enabling the characterization of variable relationships and the expression of multi-scale features. This improves the accuracy and stability of system operation status representation, enhances the precision of anomaly detection, and increases the reliability of operation monitoring.

[0136] To address the prevalent nonlinear dependencies among multiple variables in water treatment systems, this technical solution proposes an adaptive graph structure learning method based on mutual information and DBSCAN. By utilizing mutual information to measure the dependencies between monitored variables, it can effectively capture nonlinear association characteristics that are difficult to represent using traditional linear methods. Furthermore, the DBSCAN clustering algorithm is introduced to adaptively group variables and construct a variable relationship graph. Compared to existing kNN graph construction methods that require pre-setting the number of neighbors, this adaptive graph structure learning method based on mutual information and DBSCAN does not require pre-determining the number of neighbor nodes. It can adaptively determine the neighborhood range based on data distribution characteristics, allowing adjacency relationships to dynamically adjust with variable density changes, thus better reflecting the actual coupling relationship distribution characteristics among monitored variables in water treatment systems. Simultaneously, this method exhibits robustness to outliers, reducing the interference of abnormal data on the neighborhood structure and improving the stability and accuracy of graph structure learning.

[0137] This technical solution also proposes a frequency domain decoupling learning mechanism, which can map multivariate time series from the time domain to the frequency domain, characterizing the overall spectral features of the signal and thus enhancing the ability to represent global change patterns. Considering that water treatment systems exhibit both slowly changing operational trends and dynamic fluctuations caused by disturbances, and that the change patterns corresponding to different frequency components differ, a uniform modeling approach could easily lead to the mixing of different frequency information. Therefore, this invention adaptively partitions different frequency components based on spectral energy and performs feature learning separately in each frequency subspace to characterize the interaction relationships of variables at different frequency scales. While preserving key information, this technical solution can reduce the impact of noise interference, thereby improving the feature representation capability of multivariate time series. Attached Figure Description

[0138] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0139] Figure 1 This is a diagram of the overall framework structure of the present invention;

[0140] Figure 2 This is a performance comparison chart of the present invention with LSTM-VAE, OmniAnomaly, MTAD-GAT, GDN, InterFusion, Anomaly-Transformer, MEMTO, and MtsCID on the SWaT dataset. Detailed Implementation

[0141] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0142] like Figure 1 As shown, an anomaly detection method for a water treatment system based on adaptive mapping and frequency domain decoupling is applied to an industrial water treatment system. This system includes a raw water inlet unit, a reaction tank, a sedimentation tank, a filtration unit, and a disinfection unit. Each treatment unit is connected by pipelines equipped with flow sensors, pumps, and valves. Each treatment unit is equipped with sensors for pressure, water level, pH value, turbidity, and residual chlorine concentration. The method includes the following steps:

[0143] S1: Data Preparation. Real-time acquisition of monitoring variable data corresponding to each monitoring variable in each treatment unit of the water treatment system. The monitoring variables include sensors such as flow rate, pressure, water level, pH value, turbidity, and residual chlorine concentration, as well as actuators such as pumps and valves. The corresponding monitoring variable data includes the detected values ​​of flow rate, pressure, water level, pH value, turbidity, and residual chlorine concentration, and equipment status values ​​such as pump operating status and valve opening. The acquired data is treated as multivariate time-series data and standardized to convert it into standardized multivariate time-series data. In this embodiment, a portion of the monitoring variable data from the long-term stable operation of the water treatment system is extracted as a training set for model training. A portion of the real-time acquired monitoring variable data from each treatment unit of the water treatment system is divided into a validation set. In this embodiment, the first 35% is divided into the validation set, and the rest is the test set. The validation set is only used to determine the anomaly detection threshold, and then the remaining real-time acquired data is used to detect anomalies using this anomaly detection threshold.

[0144] S1-1: Standardize the monitoring variables in the multivariate time series data along the time dimension to transform the data values ​​to the [0,1] interval;

[0145] A multivariate time-series dataset is constructed. Monitoring variable data corresponding to sensors and actuators of each processing unit in the water treatment system are collected from the original file. The data is then parsed and format-converted using data processing tools. In this embodiment, the pandas library in the Python programming language is used to read and clean the original tabular data, removing irrelevant fields such as timestamps, dates, and row numbers (or any combination thereof), retaining only the status values ​​or detection values ​​reflecting the monitoring variables in the water treatment system. The processed data is then converted into a numerical array format using the numpy library, retaining only the data reflecting the status information or detection data of each component in the water treatment system. This data is stored using pickle serialization for subsequent model training and testing.

[0146] To eliminate the differences between different data formats and obtain structured data in a unified format;

[0147] Structured data is concatenated in chronological order to construct multivariate time-series data of the water treatment system. The multivariate time-series data is represented as follows: ;in, Representing time series Total length, Indicates in Constantly monitor the water treatment system A vector consisting of the observed values ​​of 1 monitored variable, where each component corresponds to an observed value of one monitored variable; and ,in, Indicates the first in the water treatment system One monitoring variable in The detected value or state value at any given time. This represents the total number of monitored variables in the water treatment system;

[0148] The Z-score method is used to standardize time series from different monitoring variables. The monitored variable is at the _ The standardized formula corresponding to time is:

[0149] ;

[0150] in, Indicates the first The monitored variable is at the _ The detected value or state value at any given time; Indicates the first The mean of the time series of each monitored variable is calculated using the following formula: ; Indicates the first The standard deviation of the time series of each monitored variable is calculated using the following formula: , Indicates the first Time series data of 1 monitored variable after Z-score standardization;

[0151] S1-2: Perform sliding window sampling on the standardized multivariate time series data;

[0152] Multivariate time-series data are sampled using a sliding window to construct time window samples for model training and anomaly detection.

[0153] At any moment At this location, the construction length is Time window sample: ,in, This indicates the preset window length, in this embodiment... Set it to 64. Indicated by The time window sample is the end time. This represents the total number of monitored variables in the water treatment system. Represents the space of real numbers; Indicates the water treatment system in The observation vector at time t, Indicates the water treatment system in The observation vector at time t, Indicates the water treatment system in The observation vector at each time step is then converted into a sample set suitable for subsequent model training and testing.

[0154] Time window samples are used to characterize continuous The multivariate monitoring variable data at each time point, the observation vector at each time point contains the detection value or state value of all monitoring variables in each treatment unit of the water treatment system at the corresponding time point. Through the above preprocessing steps, the raw data in the form of industrial time series is transformed into standardized time series data with stable quality and consistent scale, so as to serve as the input for subsequent feature modeling and anomaly detection modules.

[0155] S2: Construct an anomaly detection model for time-series data of a water treatment system. A graph representation learning module is used to represent the features of multivariate time-series data. Mutual information is used to calculate the correlation between monitoring variables, and density clustering is combined to adaptively determine the connection relationships between nodes, constructing a feature map. Based on the constructed feature map, a graph convolutional network is used to extract spatial features reflecting the dependencies between monitoring variables. These spatial features are mapped to the frequency domain, and frequency domain transformation and enhancement are performed based on spectral energy to obtain frequency-enhanced features. The spatial features and frequency-enhanced features are fused to obtain fused features. Time-series learning is performed based on the fused features to obtain a time-series feature representation, which is then reconstructed to obtain a reconstructed feature representation.

[0156] S2-1: The graph representation learning module is used to represent the features of multivariate time series data. This includes constructing a graph structure that reflects the dependencies between the monitored variables, using mutual information to calculate the correlation between variables, and using density clustering to adaptively determine the connection relationship between nodes to construct a feature graph. A graph convolutional network is used to propagate and aggregate the features of the graph structure to obtain spatial correlation features that reflect the dependencies between the monitored variables.

[0157] S2-1-1: Calculate the mutual information values ​​between the monitored variables to obtain the statistical dependencies between the variables, and construct a mutual information matrix based on the mutual information values;

[0158] Based on time window samples in multivariate time series data, the value sequence of each monitored variable within each time window is determined; whereby, for each time window sample, a corresponding feature vector is constructed based on the value sequence of each monitored variable within that time window; then the... The monitored variable is at the _ The feature vector within each time window is represented as follows: ;in, Indicates the first Time of the first The detected value or status value of each monitored variable Indicates the first Time of the first The detected value or status value of each monitored variable;

[0159] Based on the value sequence of each monitored variable within the time window, the probability distribution corresponding to each monitored variable is determined, and the marginal entropy of each monitored variable and the joint entropy between any two monitored variables are calculated based on the probability distribution; wherein, the probability distribution is obtained by discretizing the value sequence.

[0160] Determine the mutual information value between any two monitored variables based on marginal entropy and joint entropy;

[0161] The formula for calculating the mutual information value is:

[0162] ;

[0163] This formula can calculate the nonlinear correlation between any two eigenvectors, where, Indicates the first The monitoring variable and the first The mutual information value between the monitored variables is used to measure the first... The first monitoring variable and the first Non-linear dependencies among the monitored variables; Indicates the first The monitored variable is at the _ Feature vectors within a time window; Indicates the first Marginal entropy of each monitored variable; Indicates the first The monitored variable is at the _ Feature vectors within a time window Indicates the first The marginal entropy of a monitored variable is used to measure the uncertainty of a single feature vector; Indicates the first The monitoring variable and the first The joint entropy of the monitored variables is used to measure the uncertainty of the joint relationship between two eigenvectors.

[0164] Marginal entropy and joint entropy are calculated statistically based on the value sequence of each monitored variable within the time window through probability distribution.

[0165] Construct a mutual information matrix based on the mutual information values ​​between any two monitored variables. ;

[0166] in, Mutual information matrix It is a dimension The square formation, The mutual information matrix represents the total number of monitored variables in the water treatment system. Used to characterize the nonlinear relationships between various monitored variables.

[0167] S2-1-2: Based on the mutual information matrix Determine the distance matrix between the monitored variables. Based on the distance matrix, the DBSCAN clustering algorithm is used to adaptively partition each monitoring variable to obtain a cluster structure. The monitoring variables that are not assigned to any cluster are recorded as noise points and marked as -1. The adjacency relationship between the monitoring variables is determined according to the cluster structure, and an adjacency matrix is ​​constructed.

[0168] Transform the mutual information matrix into a distance matrix The formula is as follows:

[0169] ;

[0170] in, Represents the distance matrix. It is a dimensional space where all elements are 1. The square formation, This represents the total number of monitored variables in the water treatment system. Mutual information matrix; This indicates finding the maximum value, here... Represents a matrix Perform element-wise operations, when When an element in the matrix has a negative value, it is set to 0 to ensure the non-negativity of the distance matrix.

[0171] According to the distance matrix A density-based clustering algorithm was used to cluster each monitoring variable to obtain the cluster structure corresponding to each monitoring variable. The monitoring variables that were not assigned to any cluster were recorded as noise points and marked as -1.

[0172] Based on the cluster structure, construct the adjacency matrix between the monitored variables. : where, adjacency matrix It is The square array, in which The adjacency matrix represents the total number of monitored variables in the water treatment system. Used to characterize the in Structural connectivity relationships between monitored variables within each time window, adjacency matrix Each element in It is a binary indicator variable, specifically defined as follows:

[0173] ;

[0174] in, Representing the adjacency matrix The Middle Line number List the elements, Indicates the first The first monitoring variable and the first If the monitored variables are in the same cluster, otherwise... Adjacency matrix It is used to characterize the structural connectivity and potential spatial dependencies among the monitored variables, and to provide structural information for subsequent graph convolution operations; This represents the total number of monitored variables in the water treatment system;

[0175] Based on the variable relationship matrix, the connection relationships between each monitoring variable are determined, and the monitoring variables are clustered based on the density clustering algorithm to obtain the cluster structure corresponding to each monitoring variable;

[0176] S2-1-3: Based on the graph structure, a graph convolutional network is used to extract spatial features that reflect the dependencies between monitoring variables, and the data of different monitoring variables are integrated into a unified vector representation;

[0177] Let the graph structure at each time step be as follows: , among which, graph structure Represents the set of monitored variable nodes. Indicates based on the adjacency matrix A defined set of edges; based on the graph structure The neighborhood information of the feature vectors corresponding to each monitored variable is aggregated to obtain the initial feature matrix. ; where the initial feature matrix Indicates the first The feature matrix corresponding to the samples in each time window has a dimension of . × , This represents the total number of monitored variables in the water treatment system. Indicates the length of the time window;

[0178] According to the graph structure A graph convolutional network is used to extract spatial features that reflect the dependencies between monitoring variables, and the data of different monitoring variables are integrated into a unified vector representation.

[0179] The formula for calculating graph convolution is as follows:

[0180] ;

[0181] in, The graph convolutional network represents the first... The first time window The output feature matrix of the layer, and Indicates the embedding dimension; The graph convolutional network represents the first... The first time window The output feature matrix of layer +1; This represents the graph convolutional network. The trainable weight matrix corresponding to the layer; It is a non-linear activation function; The extended adjacency matrix is ​​calculated using the following formula: ,in, Indicates the first The adjacency matrix corresponding to each time window Represents the identity matrix. Representation and extended adjacency matrix The corresponding degree matrix is ​​a diagonal matrix, and the formula for calculating the diagonal elements is: ,in, Represents the th in the extended adjacency matrix Line number Column elements; This represents the result after symmetric normalization of the extended adjacency matrix;

[0182] The spatial feature vectors of each monitored variable obtained through the graph convolutional network are represented as follows:

[0183] ;

[0184] in, Indicates the first The monitored variable is at the _ Spatial feature vector at time step, Indicates the embedding dimension. This indicates the number of monitored variables in the water treatment system;

[0185] All monitored variables in the first The spatial feature vectors at each time step are concatenated to obtain the spatial feature representation:

[0186] ;

[0187] in, For the first The spatial features of each time window will be represented in the subsequent model. Record ,Right now Indicates the first The spatial feature matrix of the nth time window, the th OK Indicates the first The monitored variable is at the _ Spatial feature vectors for each time window.

[0188] S2-2: Perform frequency domain transformation on spatial features to map them to the frequency domain space to obtain frequency domain features, and perform frequency domain decoupling learning on the frequency domain features to obtain frequency domain enhanced features; fuse spatial features and frequency domain enhanced features to obtain features that fuse spatial structure information and frequency domain enhanced information;

[0189] S2-2-1: Perform a Fast Fourier Transform (FFT) on the spatial features along the time dimension to obtain the corresponding frequency domain feature representation. Then, use an adaptive threshold based on spectral energy to partition the frequency domain feature representation, obtaining low-frequency and high-frequency components. For spatial features... Its change process can be represented as follows:

[0190] ;

[0191] in, Represents the Fast Fourier Transform. Represents frequency domain features, with dimensions of ,Right now ,in The length of the frequency domain representation; Indicates the total number of monitored variables in the water treatment system; spatial characteristics Each channel undergoes independent frequency domain transformation to obtain complete frequency domain characteristics. This frequency domain feature characterizes the spectral features of the original time series at all channel levels, achieving a holistic modeling of its frequency structure.

[0192] In obtaining the frequency domain features of time series Then, its power spectrum is calculated first to characterize the energy distribution of each frequency component.

[0193] Specifically, the energy of a frequency domain signal is defined as the square of the amplitude at each frequency point, and the formula for calculating the spectral energy is:

[0194] ;

[0195] in, This indicates that the square operation is performed on each element of the matrix; The power spectrum is used to characterize the energy intensity distribution of each frequency component.

[0196] The low-frequency and high-frequency components can be obtained using the following formula:

[0197] ;

[0198] ;

[0199] in, This represents low-frequency components, used to characterize the steady-state variation trend of a water treatment system. High-frequency components are used to characterize abrupt or abnormal fluctuations in water treatment systems. This represents element-wise multiplication. This indicates an indicator function used to perform element-by-element discrimination of the power spectrum. When the condition in parentheses is true, the function at the corresponding position takes the value of 1; otherwise, it takes the value of 0. Indicates the energy adaptive threshold; This indicates that when the corresponding element in the power spectrum matrix is ​​greater than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0. This indicates that when the corresponding element in the power spectrum matrix is ​​less than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0.

[0200] S2-2-2: Model the low-frequency components and high-frequency components separately to obtain low-frequency features and high-frequency features, and then fuse the low-frequency features and high-frequency features;

[0201] Low frequency components Input a low-frequency feature modeling function to perform feature learning, and obtain a low-frequency feature representation; then input high-frequency components... The high-frequency feature modeling function is input to perform feature learning, resulting in a high-frequency feature representation;

[0202] The calculation formulas for modeling low-frequency and high-frequency components are as follows:

[0203] ;

[0204] ;

[0205] in, This represents the frequency domain feature modeling function constructed for low-frequency components; Indicates low-frequency characteristics; This represents the frequency domain feature modeling function constructed for high-frequency components; This represents high-frequency features; the modeling function is implemented using Frequency MLP (FreMLP).

[0206] Further representation of low-frequency features High-frequency feature representation To merge;

[0207] The formula for fusing low-frequency and high-frequency features is: ;

[0208] in, This represents the frequency domain features after fusion; This represents element-wise addition. It is a low-frequency characteristic. For high-frequency features, the fused frequency domain features are represented as follows: It is used to comprehensively characterize low-frequency trend information and high-frequency change information, thereby enhancing the ability to depict the operating status of water treatment systems.

[0209] Low-frequency components of the frequency domain characteristics of the water treatment system were analyzed independently. With high frequency components Decoupling modeling is performed to obtain low-frequency feature representations. High-frequency feature representation Then, by adding element-wise, a frequency-domain enhanced feature representation is generated. This representation is then converted into time-domain features via inverse Fourier transform, providing more accurate feature inputs for subsequent state prediction and anomaly detection of the water treatment system.

[0210] S2-2-3: The fused features are mapped back to the time domain space through inverse fast Fourier transform to obtain a frequency-domain enhanced time-series feature representation;

[0211] Frequency domain feature representation after fusion Perform an inverse fast Fourier transform to obtain the corresponding time-domain feature representation;

[0212] The formula for calculating the time-domain feature representation is:

[0213] ;

[0214] in, This represents the time-series characteristics after frequency domain enhancement. This represents the inverse fast Fourier transform;

[0215] S2-2-4: Representing the enhanced time-series characteristics in the frequency domain Spatial feature representation The fusion is performed through residual connections to obtain the fused feature representation;

[0216] Representing the spatial features output by the graph convolutional network The temporal feature representation after frequency domain enhancement as a residual term The features are then fused through residual connections to obtain a fused feature representation; simultaneously, a learnable scaling parameter is introduced. This is used to adjust the contribution of the enhanced temporal features in the frequency domain to the fusion process;

[0217] The calculation formula corresponding to the residual fusion process is:

[0218] ;

[0219] in, Indicates learnable parameters; The fused temporal feature representation combines structural information extracted by GCN with frequency domain enhancement features. This design enables the model to automatically learn appropriate fusion strengths on different datasets. Based on an adaptive residual mechanism, the original temporal features are fused with the temporal enhancement features to obtain the final temporal features that take into account both basic operational information and frequency domain enhancement discriminative information, providing more robust feature inputs for subsequent tasks such as water treatment system status prediction and anomaly detection.

[0220] S2-3: The features that fuse spatial structure information and frequency domain enhancement information are positionally encoded, and then input into the Transformer encoder for temporal feature learning to obtain a feature representation containing the temporal dependencies of each monitoring variable; the feature representation containing the temporal dependencies of each monitoring variable is then input into the reconstruction module, and the feature representation is reconstructed through nonlinear mapping to obtain the reconstructed feature representation;

[0221] S2-3-1: Representation of fusion features Perform positional encoding and then perform layer normalization on the position-encoded features;

[0222] The formula for calculating the position code is:

[0223] ;

[0224] ;

[0225] in, Indicates the time step position index. Indicates the feature dimension index. This represents the position encoding function. This represents the sinusoidal component of the positional encoding in the even-numbered dimensions. This represents the cosine component of the positional encoding in odd-numbered dimensions; For embedded dimensions, Represents the sine function. Represents the cosine function;

[0226] By using alternating sine and cosine position encoding, temporal position information is injected into the fused features, enabling the model to perceive the temporal sequence and long-term dependencies of the water treatment system's state.

[0227] The temporal feature representation with location information is layer-normalized to obtain the layer-normalized temporal feature representation;

[0228] The calculation formula for layer normalization of the position-encoded features is as follows:

[0229] ;

[0230] in, This represents the standard layer normalization operation. This represents the temporal features after fusion with step S2-2-4. A dimensionally consistent location encoding matrix is ​​used to represent the location information at each time step. The representation of time-series features after layer normalization is the time-series feature representation of the state values ​​or detection values ​​of each monitoring variable of the water treatment system after layer normalization with time step location information. This representation not only injects time-series location information to perceive the time sequence of the state of the water treatment system, but also stabilizes the feature distribution through layer normalization, providing standardized input for subsequent attention mechanism modeling. This represents the temporal features after fusion;

[0231] S2-3-2: Representing the time-series features after layer normalization The Transformer encoder is used to learn temporal features to capture the temporal dependencies of various monitoring variables in the water treatment system. The Transformer encoder includes a multi-head self-attention module, a residual connection and layer normalization module, and a feedforward neural network module. The output features of the Transformer encoder are then input into the reconstruction module for reconstruction.

[0232] S2-3-2-1: Representing the time-series features after layer normalization By mapping linear transformations to query vectors, key vectors, and value vectors, the input representations corresponding to each attention head in the multi-head attention mechanism are obtained;

[0233] The linear mapping formula is expressed as follows:

[0234] ;

[0235] in, The first The projection weight matrix corresponding to the attention head query vector, key vector, and value vector; Indicates the first The query vector corresponding to each attention head; Indicates the first The key vector corresponding to each attention head; Indicates the first The value vector corresponding to each attention head; Indicates the feature embedding dimension; Indicates the number of attention heads;

[0236] S2-3-2-2: Based on the query vector Key vector Sum value vector Calculate the output of each attention head to characterize the correlation between different time steps;

[0237] The formula is as follows:

[0238] ;

[0239] in, For the first The output of each attention head This is the activation function used to transform the attention weights into a probability distribution; This represents the attention relevance score matrix, which measures the correlation between any two time steps or components (such as the correlation strength between "current flow sensor" and "valve status 10 minutes ago"). For the first The value vector corresponding to each attention head;

[0240] Calculate the first based on the scaling dot product mechanism. Head's attention output ,pass By normalizing correlation scores and weighting and fusing the fusion values, we can model the dependencies between different components at different time steps in a water treatment system.

[0241] S2-3-2-3: The outputs of each attention head are concatenated and a linear transformation is used to obtain the output feature representation of the multi-head self-attention mechanism;

[0242] The calculation process is expressed by the following formula:

[0243] ;

[0244] in, This represents the output of the multi-head self-attention mechanism; This indicates the concatenation operation of the outputs of each attention head. To output the mapping weight matrix, For the total number of attention heads; This is the output of the first attention head. For the first The output of each attention head;

[0245] S2-3-2-4: Output of bullish self-attention With input feature representation Perform residual connections and layer normalization to obtain feature representations that include the temporal dependencies of each monitored variable;

[0246] The calculation process for feature representation is as follows:

[0247] ;

[0248] in, This represents the standard layer normalization operation; It is a feature representation that includes the temporal dependencies of each monitored variable, which integrates the original information of the water treatment system, frequency domain enhancement information and multi-head attention dependencies; Characteristics representing multi-head self-attention output;

[0249] It can unify water treatment characteristics of different magnitudes (such as pH value, flow rate, liquid level, pressure, state, etc.) into a similar distribution range, avoiding training instability caused by numerical differences; it accelerates model convergence, improves training efficiency, enhances the model's generalization ability to different working conditions, and allows the model to capture complex system dependencies (such as long-term water quality trends, component linkage, abnormal disturbances) while fully preserving the basic information of the original data, greatly improving the robustness and accuracy of water quality prediction, anomaly detection, and fault diagnosis.

[0250] S2-3-2-5: Feature representation containing the time-series dependencies of each monitored variable By performing a feedforward neural network transformation and combining residual connections with layer normalization, the final temporal feature representation of the Transformer encoder output is obtained.

[0251] In this embodiment, the activation function uses a Gaussian error linear unit (GELU) to improve the model's representational power; the output of the feedforward neural network (MLP) is compared with... After performing residual connections, layer normalization is performed again to obtain the output of a single-layer Transformer encoder. The calculation process is expressed by the following formula:

[0252] ;

[0253] in, This represents the feature representation after transformation by a feedforward neural network, residual connections, and layer normalization. This represents a feedforward neural network; This represents the standard layer normalization operation;

[0254] This technical solution employs the smoother GaussianErrorLinearUnit (GELU) activation function to improve the representational power and performance of the Transformer model, with the final output being... It should be noted that the above structure only represents a single layer structure in the model. In practical applications, this model is usually stacked with at least two layers to ensure that the model has sufficient representational power.

[0255] The final output of the Transformer encoder is for each time window. A vector containing temporal context information, which characterizes the temporal dependencies of multivariate time series within a window, is represented as... .

[0256] S2-3-2-6: Input the final temporal features output by the Transformer encoder into the reconstruction module, and reconstruct the features through nonlinear mapping;

[0257] The feature representation containing global temporal information is input into the reconstruction module, and a nonlinear mapping is performed on the feature representation to obtain the reconstruction result;

[0258] The reconstruction results are dimensionally adjusted to obtain a multivariate time series that corresponds to the input time series features in both time and variable dimensions.

[0259] The reconstruction process can be represented by the following formula:

[0260] ;

[0261] in, This indicates the result of the current time window reconstruction; Indicates the reconstructed mapping function; This represents the final temporal feature representation obtained by the Transformer encoder in step S2-3-2-5; that is, Figure 1 The "Temporal Reconstruction" section;

[0262] S3: Train the model constructed in step S2 using the training set to obtain the time series data anomaly detection model for the water treatment system;

[0263] S3-1: Construct a contrastive learning objective based on the similarity relationship between node pairs, and calculate the binary cross-entropy loss between the predicted probability of node pair similarity and the true label to optimize the parameters of the temporal graph convolutional network model.

[0264] The binary cross-entropy loss is as follows:

[0265] ;

[0266] in, Represents the binary cross-entropy loss function; and For positive sample node index, Positive sample node pairs represent node pairs that have a real dependency relationship in the graph structure; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates matrix transpose. This represents a measure of the similarity of vector dot products. Represents the set of positive sample node pairs, which contains node pairs in the graph structure that have real dependencies; This represents the set of negative sample node pairs, obtained through negative sampling. For negative sample node index, Negative sample node pairs represent node pairs that do not have a real dependency relationship in the graph structure. Indicates the first Within the first time window The spatial feature vectors of each monitored variable corresponding to a graph node Indicates the first Within the first time window The spatial characteristic orientation of each monitored variable corresponding to a graph node; This represents the Sigmoid activation function;

[0267] S3-2: Calculate the input time characteristics With reconstruction results The mean square error between the two values ​​is used to obtain the reconstruction loss;

[0268] The calculation is performed using the following formula:

[0269] ;

[0270] in, Indicates the reconstruction loss; This indicates the reconstruction result for the current time window; This is the temporal feature representation of S2-3-1 after position encoding and layer normalization; Denotes the square of the Frobenius norm;

[0271] The smaller the value, the smaller the input window. The quality of the reconstructed data is higher because the model mainly uses normal samples during the training phase, while the reconstruction module learns the ability to recover the original input.

[0272] Using the detected or state values ​​of the monitored variables as observations, if the reconstruction error of the observation is small, it indicates that it conforms to the normal pattern; conversely, if the reconstruction error of the observation is large, it means that the observation deviates from the normal pattern learned by the model and may be an outlier; in other words, the higher the anomaly score, the more likely the observation is to be an anomaly.

[0273] S3-3: Introducing Hyperparameters To balance the two parts of loss, the binary cross-entropy loss and the reconstruction loss, so as to achieve joint training of the model;

[0274] The following objective function is used to jointly train the model:

[0275] ;

[0276] in, Indicates the total loss; Indicates the reconstruction loss; This represents the binary cross-entropy loss; This represents a hyperparameter used to balance the reconstruction loss and the binary cross-entropy loss;

[0277] S4: Using the water treatment time series data anomaly detection model obtained in step S3, determine the anomaly judgment threshold based on the validation set, and perform anomaly detection on the water treatment time series data to be analyzed based on the anomaly judgment threshold.

[0278] S4-1: Input the test set and use the water treatment time series data anomaly detection model obtained in step S3 to obtain the reconstruction result. With the corresponding input time series features Based on the reconstruction results With the corresponding input time series features Calculate the reconstruction error for each time window; use the following formula to calculate the reconstruction error and use it as the anomaly score. :

[0279] ;

[0280] in, Indicates the first Abnormal scores within a time window; Denotes the square of the Frobenius norm;

[0281] S4-2: Based on the validation set data, calculate the anomaly score, generate candidate thresholds within the anomaly score range, iterate through these candidate thresholds, and select the candidate threshold that optimizes the evaluation index as the anomaly determination threshold. ;

[0282] By using uniform sampling, possible thresholds can be thoroughly searched across the entire score range, and the threshold that maximizes the F1 score is selected as the anomaly detection threshold based on its corresponding F1 score. .

[0283] Specifically, based on the anomaly score calculated in step S3, its value range is uniformly divided, and samples are taken at equal intervals of a fixed number of steps K between the minimum and maximum scores to obtain a candidate threshold set. , where K represents the number of samples, which is set to 700 in this embodiment.

[0284] S4-3: Determine the status of the current time window based on the anomaly detection threshold;

[0285] The status of the current time window is determined using the following formula:

[0286] ;

[0287] in, Indicates the first Anomaly detection results for each time window; This indicates the abnormal score for the corresponding time window; This indicates the anomaly detection threshold, used to detect anomalies in real-time collected data.

[0288] If the anomaly score of the current time window is greater than the anomaly judgment threshold, the data corresponding to the current time window will be sent to the operation and maintenance monitoring center of the water treatment system.

[0289] To better illustrate this technical solution, a comparison is made between this solution and other methods such as LSTM-VAE, OmniAnomaly, MTAD-GAT, GDN, InterFusion, Anomaly-Transformer, MEMTO, and MtsCID, based on the SWaT dataset.

[0290] LSTM-VAE is a Long Short-Term Memory Variational Autoencoder, from the paper "A Multimodal Anomaly Detector for Robot-Assisted Feeding Using an LSTM-based Variational Autoencoder"; OmniAnomaly is a full-dimensional anomaly model, from the paper "Robust Anomaly Detection for Multivariate Time Series through Stochastic Recurrent Neural Network"; MTAD-GAT is a multivariate time-series anomaly detection model based on graph attention networks, from the paper "Multivariate Time-series Anomaly Detection via Graph Attention Network"; GDN is a graph bias network, from the paper "Graph neural network-based anomaly detection in multivariate time series"; InterFusion is an interactive fusion model, from the paper "Multivariate Time Series Anomaly Detection and Interpretation using Hierarchical Inter-Metric and Temporal Embedding"; Anomaly-Transformer is a time-series anomaly Transformer deep learning model, from the paper "AnomalyTransformer: Time Series Anomaly". "Detection with Association Discrepancy"; MEMTO is a memory-guided Transformer deep learning model, from the paper "MEMTO: Memory-guided Transformer for Multivariate Time Series Anomaly Detection"; MtsCID is a multivariate time series cross-interaction dual network, from the paper "Multivariate Time Series Anomaly Detection by Capturing Coarse-Grained Intra- and Inter-Variate Dependencies".

[0291] The experimental results are shown in Table 1:

[0292] Table 1: Performance Comparison of Each Model on the SWaT Dataset

[0293]

[0294] As shown in Table 1, the proposed water treatment system timing anomaly detection method exhibits significant superiority in multiple comparative experiments. Experimental results demonstrate that, under the same test data conditions, the method of this invention outperforms existing mainstream methods in key evaluation indicators such as precision, recall, and F1 score. Specifically, the method of this invention achieves a detection precision of 97.51% and a recall of 98.23% on the test set, corresponding to an F1 score of 0.9787. The overall detection performance is significantly better than the comparative model. The bar chart of the data in Table 1 is shown below. Figure 2 As shown.

[0295] In comparisons with methods such as LSTM-VAE, OmniAnomaly, MTAD-GAT, GDN, InterFusion, Anomaly-Transformer, MEMTO, and MtsCID, the method of this invention achieves the best overall performance. Specifically, compared to the GDN model, the recall rate of this invention increases from 0.6812 to 0.9823, an absolute improvement of approximately 0.3011, significantly reducing the risk of false negatives while improving anomaly coverage. Compared to the MTAD-GAT model, the recall rate increases by approximately 0.2687, indicating a significantly enhanced ability to identify anomalous samples. Compared to the LSTM-VAE model, the precision increases from 0.7600 to 0.9751, an improvement of approximately 0.2151, effectively reducing the false positive rate. Compared to methods such as OmniAnomaly and InterFusion, the method of this invention improves precision by approximately 16%–17% and F1 score by approximately 15%, demonstrating superior overall detection performance.

[0296] Furthermore, compared to high recall methods such as Anomaly-Transformer, MEMTO, and MtsCID, the method of this invention improves detection accuracy while maintaining high recall, thereby achieving a better balance between false alarm control and false negative control. This results in an overall F1 score of 0.9787, which is better than all the comparison methods, indicating that the present invention achieves a better synergistic optimization effect between accuracy and completeness.

Claims

1. A method for anomaly detection in a water treatment system based on adaptive graphing and frequency domain decoupling, characterized in that, Includes the following steps: S1: Data preparation: Collect the detection values ​​or status values ​​of various monitoring variables of the water treatment system as multivariate time series data, standardize the collected multivariate time series data, and perform sliding window sampling on the standardized multivariate time series data. S2: Construct an anomaly detection model for time series data of a water treatment system. Use a graph representation learning module to represent the features of multivariate time series data to obtain spatial features that reflect the dependencies between various monitoring variables. Perform frequency domain transformation and enhancement on the spatial features to obtain frequency domain enhanced features. Spatial features and frequency domain enhancement features are fused to obtain fused features; Temporal learning is performed based on fused features to obtain temporal feature representations, and these temporal feature representations are then reconstructed to obtain reconstructed feature representations. S2-1: The graph representation learning module is used to represent the features of multivariate time series data, including constructing a graph structure that reflects the dependencies between the monitoring variables, and performing feature propagation and aggregation based on the graph structure to obtain spatial features that reflect the dependencies between the monitoring variables. The graph structure is represented as ,in, Represents a set of nodes. Represents the set of edges between nodes; S2-1-1: Calculate the mutual information values ​​between the monitored variables to obtain the statistical dependencies between the variables, and construct a mutual information matrix based on the mutual information values; The formula for calculating mutual information value is: ； in, Indicates the first The monitoring variable and the first Mutual information values ​​between the monitored variables; Indicates the first The monitored variable is at the _ Feature vectors within a time window; Indicates the first Marginal entropy of each monitored variable; Indicates the first The monitored variable is at the _ Feature vectors within a time window; Indicates the first Marginal entropy of each monitored variable; Indicates the first The monitoring variable and the first The joint entropy of the monitored variables; the marginal entropy and joint entropy are calculated statistically based on the value sequence of each monitored variable within the time window through probability distribution. The constructed mutual information matrix is ​​represented as follows: ; S2-1-2: Determine the distance matrix between monitored variables based on the mutual information matrix. Based on the distance matrix, the DBSCAN clustering algorithm is used to adaptively partition each monitoring variable to obtain a cluster structure. The monitoring variables that are not assigned to any cluster are recorded as noise points and marked as -1. The adjacency relationship between the monitoring variables is determined according to the cluster structure, and an adjacency matrix is ​​constructed. Distance matrix The corresponding calculation formula is: ； in, Let represent the distance matrix, and Let be a matrix whose elements are all 1, and It is a mutual information matrix, and Its element values ​​are the mutual information values ​​between variables; Represents a matrix Perform element-wise operations, when When an element in the matrix has a negative value, it is set to 0 to ensure the non-negativity of the distance matrix. The adjacency matrix between the constructed monitoring variables is represented as follows: ； in, Indicates the first The first monitoring variable and the first If all monitored variables are in the same cluster, otherwise The adjacency matrix is ​​used to characterize the connection relationships between the monitored variables and to provide graph structure information for subsequent graph convolution operations. S2-1-3: Based on the diagram structure A graph convolutional network is used to extract spatial features that reflect the dependencies between monitoring variables, and the data of different monitoring variables are integrated into a unified vector representation. The formula for calculating graph convolution is as follows: ； in, The graph convolutional network represents the first... The first time window The output feature matrix of the layer, and Indicates the embedding dimension; The graph convolutional network represents the first... The first time window The output feature matrix of layer +1; This represents the graph convolutional network. The trainable weight matrix corresponding to the layer; It is a non-linear activation function; The extended adjacency matrix is ​​calculated using the following formula: ,in, Indicates the first The adjacency matrix corresponding to each time window Represents the identity matrix. Representation and extended adjacency matrix The corresponding degree matrix is ​​a diagonal matrix, and the formula for calculating the diagonal elements is: ,in, Represents the th in the extended adjacency matrix Line number Column elements; This represents the result after symmetric normalization of the extended adjacency matrix; The spatial feature vectors of each monitored variable obtained through the graph convolutional network are represented as follows: ； in, Indicates the first The monitored variable is at the _ Spatial feature vector at time step, Indicates the embedding dimension. This indicates the number of monitored variables in the water treatment system; All monitored variables in the first The spatial feature vectors at each time step are concatenated to obtain the spatial feature representation: ； in, For the first The spatial features of each time window will be represented in the subsequent model. Record ,Right now Indicates the first The spatial feature matrix of the nth time window, the th OK Indicates the first The monitored variable is at the _ Spatial feature vectors of a time window; S2-2: Perform frequency domain transformation on spatial features to map them to the frequency domain space to obtain frequency domain features, and perform frequency domain decoupling learning on the frequency domain features to obtain frequency domain enhanced features; fuse spatial features and frequency domain enhanced features to obtain features that fuse spatial structure information and frequency domain enhanced information; S2-2-1: Perform a Fast Fourier Transform (FFT) on the spatial features along the time dimension to obtain the corresponding frequency domain feature representation. Then, use an adaptive threshold based on spectral energy to partition the frequency domain feature representation, obtaining low-frequency and high-frequency components. For spatial features... Its change process can be represented as follows: ； in, Represents the Fast Fourier Transform. Represents frequency domain characteristics, and ,in The length of the frequency domain representation; This represents the total number of monitored variables in the water treatment system; The formula for calculating spectral energy is: ； in, This indicates that the square operation is performed on each element of the matrix; The power spectrum is used to characterize the energy intensity distribution of each frequency component. The low-frequency and high-frequency components can be obtained using the following formula: ； ； in, This represents element-wise multiplication; This indicates an indicator function used to perform element-by-element discrimination of the power spectrum. When the condition in parentheses is true, the function at the corresponding position takes the value of 1; otherwise, it takes the value of 0. Indicates the energy adaptive threshold; This indicates that when the corresponding element in the power spectrum matrix is ​​greater than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0. This indicates that when the corresponding element in the power spectrum matrix is ​​less than the energy adaptive threshold... The value is 1 if the condition is met, otherwise it is 0. S2-2-2: Model the low-frequency components and high-frequency components separately to obtain low-frequency features and high-frequency features, and then fuse the low-frequency features and high-frequency features; The calculation formulas for modeling low-frequency and high-frequency components are as follows: ； ； in, This represents the frequency domain feature modeling function constructed for low-frequency components; Indicates low-frequency characteristics; This represents the frequency domain feature modeling function constructed for high-frequency components; Indicates high-frequency characteristics; The formula for fusing low-frequency and high-frequency features is: ； in, This represents the frequency domain characteristics after fusion; This represents element-wise addition. It is a low-frequency characteristic. High-frequency characteristics; S2-2-3: The fused features are mapped back to the time domain space through inverse fast Fourier transform to obtain a frequency-domain enhanced time-series feature representation; The formula for calculating the inverse fast Fourier transform is: ； in, This represents the time-series characteristics after frequency domain enhancement. This represents the inverse fast Fourier transform; S2-2-4: Representing the enhanced time-series characteristics in the frequency domain Spatial feature representation The fusion is performed through residual connections to obtain the fused feature representation; The calculation formula for the fusion process is as follows: ； in, This represents the temporal features after fusion; Indicates learnable parameters; S2-3: The features that fuse spatial structure information and frequency domain enhancement information are positionally encoded, and then input into the Transformer encoder for temporal feature learning to obtain a feature representation containing the temporal dependencies of each monitoring variable; the feature representation containing the temporal dependencies of each monitoring variable is then input into the reconstruction module, and the feature representation is reconstructed through nonlinear mapping to obtain the reconstructed feature representation; S2-3-1: Representation of fusion features Perform positional encoding and then perform layer normalization on the position-encoded features; The formula for calculating the position code is: ； ； in, Indicates the time step position index; Indicates the feature dimension index; This represents the position encoding function; This represents the sinusoidal component of the positional encoding in the even-numbered dimensions. This represents the cosine component of the positional encoding in odd-numbered dimensions; Indicates the feature embedding dimension. Represents the sine function. Represents the cosine function; The formula for calculating layer normalization of the position-encoded features is as follows: ； in, The time-series feature representation after normalization of the representation layer. This represents the temporal features after fusion with step S2-2-4. A dimensionally consistent location encoding matrix is ​​used to represent the location information at each time step. This represents the fused temporal feature representation. This represents the standard layer normalization operation; S2-3-2: Representing the time-series features after layer normalization The Transformer encoder is used to learn temporal features to capture the temporal dependencies of various monitoring variables in the water treatment system. The Transformer encoder includes a multi-head self-attention module, a residual connection and layer normalization module, and a feedforward neural network module. The output features of the Transformer encoder are then input into the reconstruction module for reconstruction. S3: Train the model constructed in step S2 using the training set to obtain the time series data anomaly detection model for the water treatment system; S4: Based on the time series data anomaly detection model of the water treatment system obtained in step S3, determine the anomaly judgment threshold, and perform anomaly detection on the multivariate time series data to be analyzed according to the anomaly judgment threshold.

2. The method for anomaly detection in a water treatment system based on adaptive graphing and frequency domain decoupling according to claim 1, characterized in that, Step S1 also includes the following steps: S1-1: Standardize each monitoring variable in the multivariate time series data along the time dimension; Multivariate time series data are represented as follows: ;in, Representing time series Total length, Indicates in Constantly monitor the water treatment system A vector consisting of the observed values ​​of 1 monitored variable, where each component corresponds to an observed value of one monitored variable; and ,in, Indicates the first in the water treatment system One monitoring variable in The detected value or state value at any given time. This represents the total number of monitored variables in the water treatment system; The Z-score method is used to standardize time series from different monitoring variables. The monitored variable is at the _ The standardized formula corresponding to time is: ； in, Indicates the first The monitored variable is at the _ The detected value or state value at any given time; Indicates the first The mean of the time series of each monitored variable is calculated using the following formula: ; Indicates the first The standard deviation of the time series of each monitored variable is calculated using the following formula: , Indicates the first Time series data of 1 monitored variable after Z-score standardization; S1-2: Perform sliding window sampling on the standardized multivariate time series data; Multivariate time-series data are sampled using a sliding window to construct time window samples for model training and anomaly detection; at time... At this location, the construction length is Time window sample: ;in, This indicates the preset window length. Indicated by The time window sample is the end time. Indicates the water treatment system in The observation vector at time t, Indicates the water treatment system in The observation vector at time t.

3. The method for anomaly detection in a water treatment system based on adaptive mapping and frequency domain decoupling according to claim 1, characterized in that, Step S2-3-2 also includes the following steps: S2-3-2-1: Representing the time-series features after layer normalization By transforming the input representations of each attention head in the multi-head attention mechanism into query vectors, key vectors, and value vectors through linear mapping, the results are obtained. The linear mapping formula is expressed as follows: ； in, The first The projection weight matrix corresponding to the attention head query vector, key vector, and value vector; Indicates the first The query vector corresponding to each attention head; Indicates the first The key vector corresponding to each attention head; Indicates the first The value vector corresponding to each attention head; Indicates the feature embedding dimension; Indicates the number of attention heads; S2-3-2-2: Based on the query vector Key vector Sum value vector Calculate the output of each attention head; The attention calculation process is expressed by the following formula: ； in, Indicates the first The output of each attention head; This is the activation function used to transform the attention weights into a probability distribution; This represents the attention relevance score matrix. For the first The value vector corresponding to each attention head; S2-3-2-3: The outputs of each attention head are concatenated and a linear transformation is used to obtain the output feature representation of the multi-head self-attention mechanism; The calculation process is expressed by the following formula: ； in, This represents the output of the multi-head self-attention mechanism; This indicates a concatenation operation of the output results of each attention head; This indicates the output mapping weight matrix; Indicates the number of attention heads; This represents the output of the first attention head. Indicates the first The output of each attention head; S2-3-2-4: Output of bullish self-attention With input feature representation Perform residual connections and layer normalization to obtain feature representations that include the temporal dependencies of each monitored variable; The calculation process is expressed by the following formula: ； in, A feature representation that includes the temporal dependencies of each monitored variable; This represents the standard layer normalization operation; Characteristics representing multi-head self-attention output; S2-3-2-5: Feature representation containing the time-series dependencies of each monitored variable By performing a feedforward neural network transformation and combining residual connections with layer normalization, the final temporal feature representation of the Transformer encoder output is obtained. The calculation process is expressed by the following formula: ； in, This represents the feature representation after transformation by a feedforward neural network, residual connections, and layer normalization. This represents a feedforward neural network; This represents the standard layer normalization operation; S2-3-2-6: Input the final temporal features output by the Transformer encoder into the reconstruction module for reconstruction to obtain the reconstructed feature representation; The reconstruction process can be represented by the following formula: ； in, Indicates the first The reconstruction results for each time window; Indicates the reconstructed mapping function; This represents the final temporal feature representation obtained by the Transformer encoder in step S2-3-2-5.

4. The method for anomaly detection in a water treatment system based on adaptive mapping and frequency domain decoupling according to claim 1, characterized in that, Step S3 also includes the following steps: S3-1: Construct a contrastive learning objective based on the similarity relationship between node pairs, and calculate the binary cross-entropy loss between the predicted similarity probability of node pairs and the true label. The binary cross-entropy loss function is expressed as follows: ； in, Represents the binary cross-entropy loss function; and For positive sample node index, Positive sample node pairs represent node pairs that have a real dependency relationship in the graph structure; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates the first Within the first time window The spatial feature vectors of the graph nodes corresponding to each monitored variable; Indicates matrix transpose. This represents a measure of the similarity of vector dot products. Represents the set of positive sample node pairs, which contains node pairs in the graph structure that have real dependencies; This represents the set of negative sample node pairs, obtained through negative sampling. For negative sample node index, Negative sample node pairs represent node pairs that do not have a real dependency relationship in the graph structure. Indicates the first Within the first time window The spatial feature vectors of each monitored variable corresponding to a graph node Indicates the first Within the first time window The spatial characteristic orientation of each monitored variable corresponding to a graph node; This represents the Sigmoid activation function; S3-2: Calculate the input time characteristics With reconstruction results The mean square error between the two values ​​is used to obtain the reconstruction loss; The formula for calculating the reconstruction loss is as follows: ； in, Indicates the reconstruction loss; This indicates the reconstruction result for the current time window; This is the temporal feature representation of S2-3-1 after position encoding and layer normalization; Denotes the square of the Frobenius norm; S3-3: The following objective function is used to jointly train the model: ； in, Indicates the total loss; Indicates the reconstruction loss; This represents the binary cross-entropy loss; This represents a hyperparameter used to balance the reconstruction loss and the binary cross-entropy loss.

5. A method for anomaly detection in a water treatment system based on adaptive graphing and frequency domain decoupling according to claim 1, characterized in that, Step S4 also includes the following steps: S4-1: Based on the reconstruction results With the corresponding input time series features Calculate the reconstruction error for each time window; use the following formula to calculate the reconstruction error and use it as the anomaly score. : ； in, Indicates the first Abnormal scores within a time window; Denotes the square of the Frobenius norm; S4-2: Based on the validation set data, calculate the anomaly score, generate candidate thresholds within the anomaly score range, iterate through these candidate thresholds, and select the candidate threshold that optimizes the evaluation index as the anomaly determination threshold. ; S4-3: Determine the status of the current time window based on the anomaly detection threshold; The status of the current time window is determined using the following formula: ； in, Indicates the first Anomaly detection results for each time window; This indicates the abnormal score for the corresponding time window; This indicates the threshold for anomaly detection.

Images