Multivariate time sequence anomaly detection method based on double reconstruction

Abstract

The invention discloses a multivariate time sequence anomaly detection method based on double reconstruction, and relates to the technical field of data processing and analysis. The method comprises the following steps: S1, preprocessing an input time sequence data sample; s2, converting the time sequence data preprocessed in the S1 into a sliding window sequence sample; s3, performing discrimination operation on the sliding window sequence sample in the S2 and the time sequence data sample in the S1 through a GAN model; s4, applying loss functions of an EM distance minimization generator and a maximization discriminator on the basis of the GAN model framework in the S3; s5, constructing a loss function of the discrimination network; s6, calculating an error, and combining S4 and S5 to obtain a loss function L1 of the GAN model; s7, performing secondary reconstruction to obtain a loss function L2; s8, calculating an abnormal score threshold value through the loss functions of S6 and S7; s9, distinguishing and marking the abnormal scores according to the result of S8; by adopting the multivariate time sequence anomaly detection method based on double reconstruction, the anomaly detection performance and reliability can be improved.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing and analysis, and in particular to a multivariate time series anomaly detection method based on double reconstruction. Background Art

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

[0003] After achieving great success in image generation and other fields, Generative Adversarial Networks (GANs) have also been introduced into time series anomaly detection. GANs learn the distribution of data through adversarial training and generate samples similar to real data. However, GANs have unstable problems during training, such as convergence difficulties and mode collapse, which seriously affect their performance in anomaly detection. Although some studies have made some progress by improving the structure and training methods of GANs (such as MAD-GAN combined with LSTM to model event sequence distribution), these problems cannot be completely overcome. Summary of the invention

[0004] The purpose of the present invention is to provide a multivariate time series anomaly detection method based on double reconstruction, so as to solve the problems faced by the prior art in processing multivariate time series anomaly detection, such as low accuracy, difficulty in detecting complex anomaly data, and model instability, and to improve the performance and reliability of anomaly detection.

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

[0006] S1, preprocess the input time series data sample X;

[0007] S2, applying sliding window technology to convert the preprocessed time series data in S1 into sliding window sequence samples;

[0008] S3, input the sliding window sequence samples of S2 into the GAN model to obtain the samples generated by the generator, which are passed to the discriminator and the time series data samples in S1 for discrimination operation;

[0009] S4. Based on the GAN model framework of S3, EM distance is applied to minimize the loss function of the generator and maximize the loss function of the discriminator;

[0010] S5. Remove the sigmoid function of the last layer of the discriminator and transform it into a regression task. Introduce the Lipschitz continuity constraint and construct the loss function of the discriminant 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 L of the GAN model 1 ;

[0012] S7, apply the VAE encoder to reconstruct the samples generated by the GAN model generator to obtain the loss function L 2 ;

[0013] S8, calculating the anomaly score threshold using the loss function in S6 and S7;

[0014] S9. Differentiate and mark according to the abnormal score threshold obtained in S8.

[0015] Preferably, the data preprocessing in S1 is to normalize 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, and x t is the normalized time series, ε is a 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 sliding window length.

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

[0022]

[0023] Where P r (W) is the true sample distribution, P g (z) is the sample distribution generated by the generator, V is the value function of the entire GAN, D is the discriminator, G is the generator, Indicates that the loss function of the generator should be minimized. Indicates that the loss function of the discriminator should be maximized. is real data sampled from the real data distribution, is the data sampled from the data distribution of the generator, D(w) is the result of the discriminator on the real data, and D(G(z)) is the result of the discriminator on the generated data.

[0024] Preferably, the loss function Pro of the discriminant network constructed based on the discriminator in S5 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] Where W t is the initial input data of the model, and G(Z) is the output data generated by the model.

[0030] Preferably, the loss function L of the GAN model in S6 is 1 The calculation formula is as follows:

[0031] L 1 =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 L 2 , the calculation formula is as follows:

[0033]

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

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

[0036] S81. Define anomaly score s = αL 1 +βL 2, where α and β are weight coefficients determined experimentally;

[0037] S82. Given n observations X n The abnormal probability q and select an initial threshold c;

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

[0039]

[0040] Where n is the total number of observations, N c is the number of peaks, Z q is the final threshold, is the shape parameter estimated in the Grimshaw fit, is the estimated standard deviation in the Grimshaw fit.

[0041] Preferably, the process of distinguishing and marking according to the abnormal 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 abnormal. The anomaly detection label for each dimension is defined as follows:

[0043]

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

[0045] Therefore, the present invention adopts a multivariate time series anomaly detection method based on double reconstruction with the above structure, which has the following beneficial effects compared with the prior art:

[0046] 1. The present invention combines GAN and VAE to perform dual reconstruction, which can effectively increase the difference between abnormal data and normal data, thereby improving the accuracy of anomaly detection. Experiments on multiple public datasets show that GenRecAD improves AUC and F1 scores by 2.27% and 5.11% respectively compared with existing advanced baseline methods, and the average F1 score reaches 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 model training more stable, better able to generate data reconstruction, and improving the reliability of anomaly detection;

[0048] 3. For abnormal data 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 capabilities for complex and diverse abnormal data;

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

[0050] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 A method flow chart of a multivariate time series anomaly detection method based on double reconstruction according to the present invention;

[0052] Figure 2 It is a schematic diagram of the structural framework of a multivariate time series anomaly detection method based on double reconstruction of the present invention. DETAILED DESCRIPTION

[0053] Example

[0054] In the description of the present invention, it should be noted that the terms "upper", "lower", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings, or are directions or positional relationships in which the product of the invention is usually placed when in use. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore should not be understood as a limitation on the present invention.

[0055] like Figure 1-Figure 2 As shown, a multivariate time series anomaly detection method based on double reconstruction of the present invention comprises the following steps:

[0056] S1. Preprocess the input time series data sample X and normalize 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, and x t is the normalized time series, ε is a small constant vector, and t is 1, 2, ..., n;

[0059] S2, apply the sliding window technology 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 sliding window length;

[0062] S3, input the sliding window sequence samples of S2 into the GAN model to obtain the samples generated by the generator, which are passed to the discriminator and the time series data samples in S1 for discrimination operation;

[0063] GAN data generation based on EM distance, training the GAN model, mapping time series samples to latent space, the latent space sequence as the input of the generator, through the generator to generate samples after passing to the discriminator for discrimination operation; the discriminator attempts to distinguish the generated data sequence from the actual normal training data sequence, after training iterations, the samples generated by the generator and the original data samples in the discriminator perform similarly.

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

[0065] S4. Based on the GAN model framework of S3, the EM distance is applied to minimize the generator and maximize the loss function of the discriminator. The calculation formula is as follows:

[0066]

[0067] Where P r (W) is the true sample distribution, P g (z) is the sample distribution generated by the generator, V is the value function of the entire GAN, D is the discriminator, G is the generator, Indicates that the loss function of the generator should be minimized. Indicates that the loss function of the discriminator should be maximized. is real data sampled from the real data distribution, is the data sampled from the data distribution of the generator, D(w) is the result of the discriminator on the real data, and D(G(z)) is the result of the discriminator on the generated data.

[0068] S5. Remove the sigmoid function of the last layer of the discriminator and transform it into a regression task. Introduce the Lipschitz continuity constraint and construct the loss function of the discriminant 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 formula as follows:

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

[0073] Where W t is the initial input data of the model, and G(Z) is the output data generated by the model;

[0074] Combined with the contents of S4 and S5, we can get the loss function L of the GAN model. 1 The calculation formula is as follows:

[0075] L 1 =Res+Pro; (6);

[0076] S7, apply the VAE encoder to reconstruct the samples generated by the GAN model generator to obtain the loss function L 2 , the calculation formula is as follows:

[0077]

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

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

[0080] S8, calculating the anomaly score threshold using the loss function in S6 and S7;

[0081] S81. Define anomaly score s = αL1 +βL 2 , where α and β are weight coefficients determined experimentally;

[0082] S82. Given n observations X n The abnormal probability q and select an initial threshold c;

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

[0084]

[0085] Where n is the total number of observations, N c is the number of peaks, Z q is the final threshold, is the shape parameter estimated in the Grimshaw fit, is the estimated standard deviation in the Grimshaw fit;

[0086] S9. Differentiate and mark the data according to the anomaly score threshold obtained in S8. If the anomaly score of the data is greater than the threshold, mark the data as abnormal. The anomaly detection label of each dimension is defined as follows:

[0087]

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

[0089] The specific implementation process is divided into the following three steps:

[0090] 1. Data Preprocessing Example

[0091] Data normalization: For the datasets used (such as SWaT, SMAP, MSL, UCR, etc.), take the SWaT dataset as an example, first normalize its time series data. Assume that the time series data of a certain dimension in the dataset is:

[0092] X={x 1 ,x 2 ,...,x t};

[0093] where x t Represents the observed value at time t. It is normalized by the following formula:

[0094]

[0095] Here min(X) and max(X) are the minimum and maximum values ​​in the time series of this dimension, respectively, and ε is a very small constant vector (set to ) to prevent division by zero. This step can map the data to a specific interval, enhance the stability and comparability of the data, and help with subsequent model training and analysis. For example, if the minimum value of the dimension data is 10 and the maximum value is 100, for the observation at a certain moment, the normalized value is 0.222.

[0096] Then, the sliding window technique is used to convert the time series into a sliding window sequence. The sliding window size K is set to 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 context information for the data points, allowing the model to better capture the temporal characteristics and dependencies in the data. For example, at time t = 15, the sliding window sequence contains x 6 , x 7 , ..., x 15 These 10 data points and sliding window sequences will serve as input data for subsequent model training and testing.

[0099] 2. Model training and parameter setting

[0100] During model training, in the GAN data generation part based on EM distance, four deconvolution blocks are used to define the generator and discriminator of GAN. For the SWaT dataset, the number of channels is set to 1, and the cardinality of the number of features is fixed to 8. The number of output channels of each deconvolution block is adjusted by stepwise multiplication operations. For example, the number of output channels of the first deconvolution block is 8×2=16 for the second, 16×2=32 for the third, and 32×5=16 for the fourth. The convolution kernel size is set to (3,3), the step size 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 adopts a regression task. Referring to the improvement of WGAN, the gradient disappearance and mode collapse problems are avoided.

[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. As the training progresses, the learning rate is gradually optimized according to the loss of the validation set. For example, after each epoch, if the validation set loss has not decreased for consecutive epochs, the learning rate is multiplied by . The latent space dimension is set to 16. For example, when training the SWaT dataset, the data generated by the generator is mapped to a 16-dimensional latent space for subsequent calculations.

[0102] For the POT parameter, the coefficient was set to 10 for all data sets. -4 (Here it is assumed that this coefficient performs well for each dataset in the experiment.) For the SWaT dataset, the initial threshold quantile is set to 99.3% based on the experiment. The anomaly score calculation is set to α:β=7:3, which is the optimal ratio obtained by testing multiple datasets in the experiment. For example, in calculating the anomaly score s=αL 1 +βL 2 When, for a data point in the SWaT dataset, if L 1 =0.7,L 2 =0.3, then the anomaly score s=0.7*0.5+0.3*0.3=0.44.

[0103] The SWaT dataset is divided into 80% training dataset and 20% test dataset, and the training process is set to 5 rounds (epoch). An early stopping strategy is adopted in the training process. For example, during the training process, the accuracy of the validation set is calculated after each round. If the accuracy of the validation set does not improve after consecutive rounds, the training is stopped. During the training process, the changes in the training loss and validation loss are recorded. For example, in the early stage of training, the training loss drops rapidly, indicating that the model is learning data features quickly, but there may be a risk of overfitting; as the number of training rounds increases, the loss gradually stabilizes and the model gradually converges.

[0104] 3. Model evaluation and result analysis

[0105] (1) Use precision, recall, area under the receiver operating characteristic curve (ROC / AUC) and F1 score to evaluate the detection performance of the GenRecAD model on the SWaT dataset. Assume that during the detection process, for a certain abnormal event, the model correctly detects TP=80 abnormal data points, incorrectly detects FP=20 normal data points as abnormal, and FN=10 abnormal data points are not detected. Then the precision P=0.8, the recall R=0.889, and the F1 score F1=0.842.

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

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

[0108]

[0109]

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

[0111] Therefore, the present invention adopts a multivariate time series anomaly detection method based on double reconstruction of the above structure, which solves the problems faced by the prior art in processing multivariate time series anomaly detection, such as low accuracy, difficulty in detecting complex abnormal data, and model instability, and improves the performance and reliability of anomaly detection.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A multivariate time series anomaly detection method based on dual reconstruction, characterized by: S1, preprocess the input time series data sample X; S2, applying sliding window technology to convert the preprocessed time series data in S1 into sliding window sequence samples; S3, input the sliding window sequence samples of S2 into the GAN model to obtain the samples generated by the generator, which are passed to the discriminator and the time series data samples in S1 for discrimination operation; S4. Based on the GAN model framework of S3, EM distance is applied to minimize the loss function of the generator and maximize the loss function of the discriminator; S5. Remove the sigmoid function of the last layer of the discriminator and transform it into a regression task. Introduce the Lipschitz continuity constraint and construct the loss function of the discriminant network. 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; S7, applying the VAE encoder to reconstruct the samples generated by the GAN model generator to obtain the loss function L2; S8, calculating the anomaly score threshold using the loss function in S6 and S7; S9. Differentiate and mark according to the abnormal score threshold obtained in S8.

2. The multivariate time series anomaly detection method based on double reconstruction according to claim 1, characterized in that: The data preprocessing in S1 is to normalize the input time series data. The calculation formula is as follows: Where min(X) and max(X) are the minimum and maximum vectors in the time series, respectively, and x t is the normalized time series, ε is a small constant vector, and t is 1, 2, ..., n.

3. The multivariate time series anomaly detection method based on double reconstruction according to claim 2, characterized in that: 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 sliding window length.

4. The multivariate time series anomaly detection method based on double reconstruction according to claim 3 is characterized by: The loss function calculation formula of the generator and discriminator in S4 is as follows: Where P r (W) is the true 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, G is the generator, Indicates that the loss function of the generator should be minimized. Indicates that the loss function of the discriminator should be maximized. is real data sampled from the real data distribution, is the data sampled from the data distribution of the generator, D(w) is the result of the discriminator on the real data, and D(G(z)) is the result of the discriminator on the generated data.

5. The multivariate time series anomaly detection method based on double reconstruction according to claim 4 is characterized in that: The calculation formula of the loss function Pro of the discriminant network constructed based on the discriminator in S5 is as follows: where f ω (w) is the output of the discriminator.

6. The multivariate time series anomaly detection method based on double reconstruction according to claim 5, characterized in that: The formula for calculating the error Res in S6 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.

7. The multivariate time series anomaly detection method based on double reconstruction according to claim 6, characterized in that: The loss function L1 calculation formula of the GAN model in S6 is as follows: L1=Res+Pro; (6).

8. The multivariate time series anomaly detection method based on double reconstruction according to claim 7, characterized in that: 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: Where T i is the input sample of the encoder, i is 1, 2, ... n; σ i is the standard deviation of the ith sample, μ i is the mean of the ith sample, e i is the latent variable obtained by sampling, μ' i is the mean of the reconstructed data of the decoder output i-th sample.

9. The multivariate time series anomaly detection method based on double reconstruction according to claim 8, characterized in that: The process of calculating the anomaly score in S8 is as follows: S81. Define anomaly score s = αL1 + βL2, where α and β are weight coefficients determined by experiments; S82. Given n observations X n The abnormal 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. 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 shape parameter estimated in the Grimshaw fit, is the standard deviation of the estimate in the Grimshaw fit.

10. The multivariate time series anomaly detection method based on double reconstruction according to claim 9, characterized in that: The process of distinguishing and marking according to the anomaly score threshold in S9 is as follows: If the anomaly score of the data is greater than the threshold, the data is marked as abnormal. The anomaly detection label for each dimension is defined as follows: where y i is the abnormal label of the i-th data, s i is the corresponding anomaly score, and POT is the calculation method of the threshold.

Images