A Multivariate Time Series Anomaly Detection Method Based on Dual Reconstruction

Abstract

This invention discloses a multivariate time series anomaly detection method based on dual reconstruction, belonging to the field of data processing and analysis technology. The method comprises the following steps: S1, preprocessing the input time series data samples; S2, converting the preprocessed time series data from S1 into sliding window sequence samples; S3, performing a discrimination operation between the sliding window sequence samples from S2 and the time series data samples from S1 using a GAN model; S4, applying the loss functions of an EM distance minimization generator and a maximum discriminator based on the GAN model framework from S3; S5, constructing the loss function of the discriminator network; S6, calculating the error and combining S4 and S5 to obtain the loss function L1 of the GAN model; S7, performing secondary reconstruction to obtain the loss function L2; S8, calculating the anomaly score threshold using the loss functions from S6 and S7; and S9, classifying and labeling the anomaly scores based on the results of S8. This invention employs the above-mentioned multivariate time series anomaly detection method based on dual reconstruction, which can improve the performance and reliability of anomaly detection.

Description

Technical Field

[0001] This invention relates to the field of data processing and analysis technology, and in particular to a multivariate time series anomaly detection method based on dual reconstruction. Background Technology

[0002] With the development of machine learning technology, deep learning methods have received widespread attention in the field of anomaly detection. Among them, reconstruction-based methods have become one of the research hotspots. These methods aim to reconstruct input data by learning normal patterns in the data and determine whether the data is abnormal based on the reconstruction error. For example, autoencoders (AEs) and their variants have been widely used in time series anomaly detection. However, AEs themselves have some limitations, such as being prone to overfitting and overgeneralization, leading to a decrease in the ability to detect anomalous data. Although many researchers have proposed various improvement schemes, such as adding memory modules (e.g., TSMAE) and regularizing the latent space (e.g., AAE), these methods still cannot completely solve the problems of AEs in anomaly detection.

[0003] Following its significant success in image generation and other fields, Generative Adversarial Networks (GANs) have been introduced into time series anomaly detection. GANs learn the distribution of data through adversarial training, generating samples similar to real data. However, GANs suffer from instability during training, such as convergence difficulties and pattern collapse, which severely impacts their performance in anomaly detection. While some research has made progress by improving GAN structure and training methods (e.g., combining MAD-GAN with LSTM to model event sequence distributions), these problems cannot be completely overcome. Summary of the Invention

[0004] The purpose of this invention is to provide a multivariate time series anomaly detection method based on dual reconstruction, which solves the problems of low accuracy, difficulty in detecting complex anomaly data, and model instability in the existing technology when processing multivariate time series anomaly detection, thereby improving the performance and reliability of anomaly detection.

[0005] To achieve the above objectives, this invention provides a multivariate time series anomaly detection method based on dual reconstruction, comprising the following steps:

[0006] S1. Preprocess the input time series data sample X;

[0007] S2. Apply the sliding window technique to convert the preprocessed time series data in S1 into sliding window sequence samples;

[0008] S3. Input the sliding window sequence samples from S2 into the GAN model to obtain the samples generated by the generator, and then pass them to the discriminator to perform a discrimination operation with the time series data samples in S1.

[0009] S4. Based on the GAN model framework of S3, apply the loss function of minimizing the generator and maximizing the discriminator by the EM distance;

[0010] S5. Remove the sigmoid function from the last layer of the discriminator to transform it into a regression task, introduce the Lipschitz continuity constraint, and construct the loss function of the discriminator network;

[0011] S6. Find the generator sample that is most similar to the sliding window sequence sample, calculate the error, and combine the contents of S4 and S5 to obtain the loss function L1 of the GAN model;

[0012] S7. The loss function L2 is obtained by performing secondary reconstruction on the samples generated by the GAN model generator using the VAE encoder.

[0013] S8. Calculate the anomaly score threshold using the loss functions in S6 and S7;

[0014] S9. Distinguish and mark the abnormal scores based on the threshold values ​​obtained in S8.

[0015] Preferably, in S1, data preprocessing involves normalizing the input time series data, and the calculation formula is as follows:

[0016]

[0017] Where min(X) and max(X) are the minimum and maximum vectors in the time series, respectively, x t For the normalized time series, ε is a very small constant vector, and t is 1, 2, ..., n.

[0018] Preferably, in S2, the time series is converted into a sliding window sequence W:

[0019] W = {x t-K+1 ,x t-K+2 ,...,x t}; (2)

[0020] Where K is the length of the sliding window.

[0021] Preferably, the loss function calculation formulas for the generator and discriminator in S4 are as follows:

[0022]

[0023] Where P r (W) is the true sample distribution, P g (z) represents the sample distribution generated by the generator, V is the value function of the entire GAN, D is the discriminator, and G is the generator. This indicates that we want to minimize the generator's loss function. This indicates that we want to maximize the loss function of the discriminator. It is real data sampled from the real data distribution. The data is sampled from the generator's data distribution. D(w) is the discriminator's judgment result on the real data, and D(G(z)) is the discriminator's judgment result on the generated data.

[0024] Preferably, the loss function Pro in S5, which is based on the discriminator to construct the discriminant network, is calculated as follows:

[0025]

[0026] Where f ω (w) is the output of the discriminator.

[0027] Preferably, the formula for calculating the error Res in S6 is as follows:

[0028] Res=|W t -G(Z)|; (5)

[0029] Among them W t G(Z) is the initial input data of the model, and G(Z) is the output data generated by the model.

[0030] Preferably, the L1 loss function of the GAN model in S6 is calculated using the following formula:

[0031] L1 = Res + Pro; (6).

[0032] Preferably, in S7, the Evidence Lower Bound (ELBO) is applied in the VAE encoder to train the neural network to calculate the loss function L2, and the calculation formula is as follows:

[0033]

[0034] Where T i These are the input samples of the encoder, where i is 1, 2, ..., n; σ i μ is the standard deviation of the i-th sample. i It is the mean of the i-th sample, e i It is the latent variable obtained by sampling, μ' i It is the mean of the reconstructed data of the i-th sample output by the decoder.

[0035] Preferably, the process of calculating the outlier score in S8 is as follows:

[0036] S81. Define the anomaly score s = αL1 + βL2, where α and β are weighting coefficients determined experimentally.

[0037] S82, Given n observations X nThe abnormal probability q is determined and an initial threshold c is selected.

[0038] S83. Use the Grimshaw technique to fit the Pareto distribution to infer the possible distribution of extreme values ​​and calculate the outlier score threshold. The calculation formula is as follows:

[0039]

[0040] Where n is the total number of all observations, N c It is the peak number, Z q It is the final threshold. These are the shape parameters estimated in the Grimshaw fitting. It is the standard deviation estimated in the Grimshaw fit.

[0041] Preferably, the process of distinguishing and marking based on the anomaly score threshold in S9 is as follows:

[0042] If the anomaly score of the data is greater than the threshold, the data is marked as an anomaly. The anomaly detection label for each dimension is defined as follows:

[0043]

[0044] Where y i It is the anomaly label for the i-th data, s i It is the corresponding abnormal score, and POT is the threshold calculation method.

[0045] Therefore, the multivariate time series anomaly detection method based on dual reconstruction using the above-described structure has the following advantages compared with the prior art:

[0046] 1. This invention combines GAN and VAE methods for dual reconstruction, which can effectively increase the difference between anomalous and normal data, thereby improving the accuracy of anomaly detection. Experiments on multiple public datasets show that GenRecAD improves AUC and F1 score by 2.27% and 5.11% respectively compared to existing state-of-the-art baseline methods, with an average F1 score of 0.9629;

[0047] 2. By using EM distance to optimize GAN, the common convergence difficulties and mode collapse problems in GAN training are avoided, making the model training more stable, generating better data reconstruction, and improving the reliability of anomaly detection.

[0048] 3. For abnormal data that is similar to normal data, the method of the present invention can amplify its reconstruction error through secondary reconstruction, reduce the possibility of missing abnormal data, and have better detection capability for complex and diverse abnormal data.

[0049] 4. The POT method based on extreme value theory is used to automatically and dynamically select the threshold, which can effectively cope with real-time changing streaming data, overcome the problem that static fixed thresholds cannot adapt to dynamic data changes, and improve the adaptability of the model in practical applications.

[0050] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0051] Figure 1 This is a flowchart of a multivariate time series anomaly detection method based on dual reconstruction according to the present invention.

[0052] Figure 2 This is a schematic diagram of the structural framework of a multivariate time series anomaly detection method based on dual reconstruction according to the present invention. Detailed Implementation

[0053] Example

[0054] In the description of this invention, it should be noted that the terms "upper," "lower," "inner," "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship in which the product of this invention is usually placed when in use. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting this invention.

[0055] like Figures 1-2 As shown, the present invention provides a multivariate time series anomaly detection method based on dual reconstruction, comprising the following steps:

[0056] S1. Preprocess the input time series data sample X by normalizing the input time series data. The calculation formula is as follows:

[0057]

[0058] Where min(X) and max(X) are the minimum and maximum vectors in the time series, respectively, x t For the normalized time series, ε is a very small constant vector, and t is 1, 2, ..., n;

[0059] S2. Apply the sliding window technique to convert the preprocessed time series data in S1 into sliding window sequence samples W as follows:

[0060] W = {x t-K+1 ,x t-K+2 ,...,x t}; (2)

[0061] Where K is the length of the sliding window;

[0062] S3. Input the sliding window sequence samples from S2 into the GAN model to obtain the samples generated by the generator, and then pass them to the discriminator to perform a discrimination operation with the time series data samples in S1.

[0063] GAN data generation based on EM distance involves training a GAN model to map time series samples to a latent space. The latent space sequence serves as the input to the generator. After the generator generates samples, they are passed to the discriminator for discrimination. The discriminator attempts to distinguish between the generated data sequence and the actual normal training data sequence. Through training iterations, the samples generated by the generator and the original data samples are made to behave similarly in the discriminator.

[0064] Based on the standard GAN framework, this invention uses EM distance instead of JS divergence and KL divergence to define the generator loss, in order to solve the gradient vanishing problem in GAN training.

[0065] S4. Based on the GAN model framework in S3, apply the loss function of minimizing the EM distance generator and maximizing the discriminator. The calculation formula is as follows:

[0066]

[0067] Where P r (W) is the true sample distribution, P g (z) represents the sample distribution generated by the generator, V is the value function of the entire GAN, D is the discriminator, and G is the generator. This indicates that we want to minimize the generator's loss function. This indicates that we want to maximize the loss function of the discriminator. It is real data sampled from the real data distribution. The data is sampled from the generator's data distribution. D(w) is the discriminator's judgment result on the real data, and D(G(z)) is the discriminator's judgment result on the generated data.

[0068] S5. Remove the sigmoid function from the last layer of the discriminator to transform it into a regression task. Introduce the Lipschitz continuity constraint and construct the loss function of the discriminator network. The Pro calculation formula is as follows:

[0069]

[0070] Where f ω (w) is the output of the discriminator;

[0071] S6. Find the generator sample that is most similar to the sliding window sequence sample, and calculate the error Res using the following formula:

[0072] Res=|Wt -G(Z)|; (5)

[0073] Among them W t G(Z) represents the initial input data of the model, and G(Z) represents the output data generated by the model.

[0074] Combining the content of S4 and S5, the formula for calculating the L1 loss function of the GAN model is as follows:

[0075] L1 = Res + Pro; (6);

[0076] S7. The loss function L2 obtained by applying the VAE encoder to perform secondary reconstruction on the samples generated by the GAN model generator is calculated as follows:

[0077]

[0078] Where V i These are the input samples of the encoder, where i is 1, 2, ..., n; σ i μ is the standard deviation of the i-th sample. i It is the mean of the i-th sample, e i It is the latent variable obtained by sampling, μ' i It is the mean of the reconstructed data of the i-th sample output by the decoder;

[0079] The encoder uses a neural network to obtain the latent variable e. i Approximate posterior distribution followed The parameter (usually the mean μ of the Gaussian distribution) i and variance ); sample one ∈ i Calculate e i =μ i +σ i ⊙∈ i As input to the decoder, the fitted likelihood distribution p θ (T|e i ), returns the parameter μ' of the distribution it follows. i and .

[0080] S8. Calculate the anomaly score threshold using the loss functions in S6 and S7;

[0081] S81. Define the anomaly score s = αL1 + βL2, where α and β are weighting coefficients determined experimentally.

[0082] S82, Given n observations X n The abnormal probability q is determined and an initial threshold c is selected.

[0083] S83. Use the Grimshaw technique to fit the Pareto distribution to infer the possible distribution of extreme values ​​and calculate the outlier score threshold. The calculation formula is as follows:

[0084]

[0085] Where n is the total number of all observations, N c It is the peak number, Z q It is the final threshold. These are the shape parameters estimated in the Grimshaw fitting. It is the standard deviation estimated in the Grimshaw fit;

[0086] S9. Based on the anomaly score threshold obtained in S8, distinguish and label the data. If the anomaly score of the data is greater than the threshold, then mark the data as an anomaly. The anomaly detection label for each dimension is defined as follows:

[0087]

[0088] Where y i It is the anomaly label for the i-th data, s i It is the corresponding abnormal score, and POT is the threshold calculation method.

[0089] The specific implementation process consists of the following three steps:

[0090] 1. Data Preprocessing Example

[0091] Data normalization: For the dataset used (such as SWaT, SMAP, MSL, UCR, etc.), taking the SWaT dataset as an example, the time series data is first normalized. Assume the time series data for one dimension of this dataset is as follows:

[0092] X = {x1, x2, ..., x} t};

[0093] Where x t This represents the observed value at time t. Normalization is performed using the following formula:

[0094]

[0095] Here, min(X) and max(X) represent the minimum and maximum values ​​of the time series in this dimension, respectively, and ε is a very small constant vector (denoted as ε) to prevent division by zero. This step maps the data to a specific interval, enhancing the stability and comparability of the data, which is helpful for subsequent model training and analysis. For example, if the minimum value of the data in this dimension is 10 and the maximum value is 100, the normalized value for an observation at a certain time would be 0.222.

[0096] Then, the sliding window technique is used to convert the time series into a sliding window sequence. The sliding window size is set to K = 10. For the time series X, the generated sliding window sequence W is as follows:

[0097] W = {x t-K+1 ,x t-K+2 ,...,x t};

[0098] This processing method provides local contextual information for the data points, enabling the model to better capture temporal features and dependencies in the data. For example, for time t=15, the sliding window sequence contains x6, x7, ..., x 15 These 10 data points. These sliding window sequences will serve as input data for subsequent model training and detection.

[0099] 2. Model Training and Parameter Setting

[0100] During model training, in the GAN data generation part based on EM distance, four deconvolutional blocks are used to define the generator and discriminator of the GAN. For the SWAT dataset, the number of channels is set to 1, and the cardinality of the number of features is fixed at 8. The number of output channels of each deconvolutional block is adjusted by stepwise multiplication operations. For example, the first deconvolutional block has 1 output channel, the second has 8×2=16, the third has 16×2=32, and the fourth has 32×5=16. The kernel size is set to (3,3), the stride is (2,2), the padding is 1, and no bias is set. The last layer of the discriminator does not use the sigmoid activation function, but instead adopts a regression task, referencing improvements to WGAN, to avoid the gradient vanishing and mode collapse problems.

[0101] In the VAE-based secondary reconstruction part, the AdamW optimizer is used to train the model. The initial learning rate is set to 0.01, and it is gradually optimized based on the loss on the validation set as training progresses. For example, if the validation set loss does not decrease for several consecutive epochs, the learning rate is multiplied by 0.01 to decay the learning rate. The hidden space dimension is set to 16; for example, when training on the SWaT dataset, the data generated by the generator is mapped to a 16-dimensional latent space for subsequent computation.

[0102] For the POT parameter, the coefficient is set to 10 for all datasets. -4(It is assumed here that this coefficient performs well across the datasets in the experiments). For the SWaT dataset, the initial threshold quantile was set to 99.3% based on experiments. The outlier score was calculated using an α:β ratio of 7:3, which is the optimal ratio obtained through testing on multiple datasets in the experiments. For example, when calculating the outlier score s = αL1 + βL2, for a data point in the SWaT dataset, if L1 = 0.7 and L2 = 0.3, then the outlier score s = 0.7 * 0.5 + 0.3 * 0.3 = 0.44.

[0103] The SWaT dataset was divided into an 80% training dataset and a 20% test dataset, with the training process set to 5 epochs. An early stopping strategy was employed during training; for example, the accuracy on the validation set was calculated after each epoch. Training was stopped if the validation set accuracy did not improve after several consecutive epochs. During training, the changes in training loss and validation loss were recorded. For instance, in the early stages of training, the training loss decreased rapidly, indicating that the model was quickly learning the data features, but there was a potential risk of overfitting. As the number of training epochs increased, the loss gradually stabilized, and the model gradually converged.

[0104] 3. Model Evaluation and Result Analysis

[0105] (1) The detection performance of the GenRecAD model on the SWaT dataset was evaluated using precision, recall, area under the receiver operating characteristic curve (ROC / AUC), and F1 score. Assume that during detection, for a given anomalous event, the model correctly detected 80 anomalous data points (TP), incorrectly detected 20 normal data points (FP), and failed to detect 10 anomalous data points (FN). Then, the precision (P) = 0.8, the recall (R) = 0.889, and the F1 score (F1) = 0.842.

[0106] (2) GenRecAD was compared with other state-of-the-art models (such as TranAD, CAE-M, MAD_GAN, etc.) on the SWaT dataset. Experimental results show that GenRecAD has a higher F1 score than other models on this dataset. For example, TranAD has an F1 score of 0.56 on the SWaT dataset, while GenRecAD has an F1 score of 0.9970. This indicates that GenRecAD has higher accuracy and reliability in detecting outliers and can more effectively identify anomalies in the data. Table 1 shows a comparison of our method with other models on four datasets.

[0107] Table 1. Comparison with other models on four datasets

[0108]

[0109]

[0110] The importance of each module in the model was analyzed through ablation experiments. For example, on the SWAT dataset, when KL divergence calculation was not used (i.e., MSE calculation was used), the F1 score dropped significantly from 0.9970 to 0.6631, a decrease of 36.5%, indicating that KL divergence is essential for VAE loss function calculation. When a simple encoder-decoder structure was used instead of VAE, the F1 score dropped by an average of 10.3%, indicating that VAE has a significant improvement in anomaly detection compared to simple AE. When GAN was not optimized, mode collapse and non-convergence problems occurred during training, and the F1 score dropped significantly on the SWAT dataset, highlighting the importance of optimizing GAN model performance.

[0111] Therefore, the present invention provides a multivariate time series anomaly detection method based on dual reconstruction using the above-mentioned structure, which solves the problems of low accuracy, difficulty in detecting complex anomaly data, and model instability faced by existing technologies in processing multivariate time series anomaly detection, thereby improving the performance and reliability of anomaly detection.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1.A method for detecting anomalies in multivariate time series based on dual reconstruction, characterized by: S1, preprocessing the input time series data sample X; S2, converting the preprocessed time series data in S1 into a sliding window sequence sample using the sliding window technique; S3, inputting the sliding window sequence sample in S2 into the GAN model to obtain the sample generated by the generator and passing it to the discriminator for discrimination operation with the time series data sample in S1; S4, applying the loss function of the generator and the discriminator to minimize the EM distance based on the GAN model framework in S3; The loss function calculation formula of the generator and the discriminator is as follows: where P r (W) is the real sample distribution, P g (z) is the sample distribution generated by the generator, V is the value function of the entire GAN model, D is the discriminator, and G is the generator, indicates that the loss function of the generator is to be minimized, indicates that the loss function of the discriminator is to be maximized, is real data sampled from the real data distribution, is data sampled from the data distribution of the generator, D(w) is the discrimination result of the discriminator on the real data, and D(G(z)) is the discrimination result of the discriminator on the generated data. S5, removing the sigmoid function in the last layer of the discriminator to convert it into a regression task, introducing Lipschitz continuity constraint, and constructing the loss function of the discriminator network; The loss function Pro of the discriminator network is calculated as follows: where f ω (w) is the output of the discriminator; S6, finding the most similar generator sample to the sliding window sequence sample, calculating the error, and combining the contents of S4 and S5 to obtain the loss function L1 of the GAN model; The error Res formula is as follows: Res = |W t -G(Z) | ; (5) where W t is the initial input data of the model, and G(Z) is the output data generated by the model; The loss function L1 of the GAN model is calculated as follows: L1=Res+Pro; (6); S7, applying VAE encoder to the sample generated by the generator of the GAN model to obtain the loss function L2 of the secondary reconstruction; The loss function L2 is calculated by training the neural network using the evidence lower bound ELBO in the VAE encoder, and the calculation formula is as follows: where T i is the input sample of the encoder, i is 1, 2,... n; σ i is the standard deviation of the i-th sample, μ i is the mean of the i-th sample, e i is the latent variable sampled, μ' i is the mean of the reconstructed data of the i-th sample output by the decoder; S8, calculating the anomaly score threshold value through the loss function in S6 and S7; S9, distinguishing and marking according to the anomaly score threshold value obtained in S8. 2.The method for detecting anomalies in multivariate time series based on dual reconstruction according to claim 1, characterized by: The data preprocessing in S1 is to normalize the input time series data, and the calculation formula is as follows: where min(X) and max(X) are the minimum and maximum vectors in the time series, respectively, x t is the normalized time series, and ε is a small constant vector, t = 1, 2,..., n. 3.The method for detecting anomalies in multivariate time series based on dual reconstruction according to claim 1, characterized by: In S2, the time series is converted into a sliding window sequence W: W = {x t-K+1 ,x t-K+2 ,...,x t} ; (2) Where K is the length of the sliding window. 4.The method for detecting anomalies in multivariate time series based on dual reconstruction according to claim 1, characterized by: The process of calculating the anomaly score in S8 is as follows: S81, define the anomaly score s=αL1+βL2, where α and β are weight coefficients determined by experiment; S82, given n observations X n an anomaly probability q and select an initial threshold c; S83, use the Grimshaw technique to fit the Pareto distribution to infer the possible distribution of extreme values and calculate the anomaly score threshold value, and the calculation formula is as follows: where n is the total number of observations, N c is the number of peaks, Z q is the final threshold, is the estimated shape parameter in the Grimshaw fit, is the estimated standard deviation in the Grimshaw fit. 5.The method for detecting anomalies in multivariate time series based on dual reconstruction according to claim 1, characterized by: The process of distinguishing and marking according to the anomaly score threshold value in S9 is as follows: If the anomaly score of the data is greater than the threshold value, the data is marked as abnormal, and the anomaly detection label of each dimension is defined as follows: where y i is the anomaly label of the ith data, s i is its corresponding anomaly score, and POT is the method of calculating the threshold.

Images