Hydrological time sequence anomaly detection method based on spatial-temporal characteristics

Abstract

The invention discloses a hydrological time sequence anomaly detection method based on spatio-temporal characteristics. The method comprises the following steps: firstly, according to digital elevation model (DEM) data of a drainage basin, extracting spatial information such as elevation and gradient in the drainage basin and a river channel distance between hydrological stations, and constructing a spatial characteristic matrix of the drainage basin according to the spatial information; meanwhile, accumulated convergence time is calculated to preprocess the input data; constructing a variational auto-encoder (VAE) according to the spatial feature matrix and the preprocessed input data, and learning spatial and temporal features in combination with a time domain convolutional network (TCN) and a graph convolutional neural network (GCN); and finally, inputting the test data set into the pre-trained network model, calculating an abnormal score, and judging abnormal data in the hydrological time sequence according to the abnormal score. According to the invention, the problems of difficult detection of multi-site hydrological time series data, difficult learning of spatial-temporal characteristics and the like in the prior art can be solved.

Description

Technical Field

[0001] The invention relates to a hydrological time series anomaly detection method based on spatiotemporal characteristics, and belongs to the technical field of hydrological data mining. Background Art

[0002] With the rapid development of Internet of Things and Big Data technologies, the scale of time series data is growing exponentially. Time series often hide valuable information for users, so anomaly detection in time series has broad prospects in the fields of medicine, hydrology, transportation, etc.

[0003] In the field of hydrology, correct hydrological data can provide knowledge support for flood control and drought relief, water conservancy project management, etc. However, on the one hand, with the vigorous development of water conservancy, more and more data are stored in the database, and data transmission anomalies, manual operation errors, and measurement instrument failures will lead to a large amount of abnormal data; on the other hand, the existing anomaly detection in the field of hydrology is mostly aimed at a single hydrological station [Yu Yufeng, Zhu Yuelong, Wan Dingsheng, et al. Hydrological time series anomaly detection based on sliding window prediction [J]. Computer Applications, 2014, 34(08): 2217-2220+2226.], without considering the spatiotemporal connection between multiple hydrological stations. In addition, the current method of using neural networks to extract spatiotemporal features [Ding C, Sun S, Zhao J. MST-GAT: A multimodal spatial-temporal graph attention network for timeseries anomaly detection [J]. Information Fusion, 2023, 89: 527-536.] does not fully explore the hidden spatiotemporal information between stations.

[0004] In order to solve the above problems, hydrological time series anomaly detection based on spatiotemporal characteristics aims to construct the spatiotemporal characteristics between hydrological stations through digital elevation data (DEM) given the time series data of multiple hydrological stations. And combine the data of surrounding stations to perform time series anomaly detection. It helps to improve data quality and facilitate future research in the field of hydrology. Summary of the invention

[0005] Purpose of the invention: The present invention provides a hydrological time series anomaly detection method based on spatiotemporal characteristics, which combines a digital elevation model (DEM), a graph convolutional neural network (GCN), and a variational autoencoder (VAE) to realize hydrological time series anomaly detection based on spatiotemporal characteristics, and can solve the difficulties in spatiotemporal feature mining in the existing technical solutions, as well as the inability to fully detect hidden anomalies in the hydrological time series.

[0006] Technical solution: The present invention provides a hydrological time series anomaly detection method based on spatiotemporal characteristics, comprising the following steps:

[0007] Step S1, extract the elevation, slope, and geographic distance between stations and other information of stations in the basin from the digital elevation model (DEM) data. And use the graph structure G = (V, E) to represent the rainfall stations in the basin and the relationship between stations, where the vector V represents the rainfall stations in the basin; E represents the edge, and the weight of the edge is the degree of association between the rainfall stations.

[0008] Step S2: First, the grid flow velocity is calculated using the Manning formula according to the river distance, river width and elevation information; then the cumulative confluence time is calculated according to the residence time of the rainfall in the grid; finally, the time series is aligned in the time dimension according to the confluence time from each rainfall station to the station to be tested.

[0009] Step S3: Use an improved variational autoencoder (VAE) to build a network model, where the encoder uses a graph convolutional neural network (GCN) to learn spatial features and a time domain convolutional network (TCN) to learn time series features, and the optimal network model is obtained after training.

[0010] Step S4: Input the test set data into the model trained in step S3 to obtain the predicted value and generated value of the real data, and calculate the anomaly score based on the difference between the predicted value, generated value and the real value.

[0011] Step S5, design a method to maximize the F1 Score to calculate the threshold. First, calculate the anomaly score using the training set data to obtain the initial threshold. Then, iteratively increase the threshold on the validation set to maximize the F1 Score. The threshold corresponding to the maximum F1 Score is the final threshold. Finally, compare the anomaly score in step S4 with the threshold. If the anomaly score is greater than the threshold, it is judged as abnormal data.

[0012] Furthermore, the detailed process of the spatial feature extraction in step S1 is as follows:

[0013] Step S1.1: Use ArcGIS to mark the station locations from the U-Tokyo MERIT DEM digital elevation data and calculate the Euclidean distance between stations and spatial information such as river length. Since there is spatial autocorrelation between stations, the Euclidean distance between stations and the reachable distance to the river can be used to construct a spatial relationship graph between stations. In the graph structure, graph nodes are used to represent hydrological stations, and edges are used to store the relationship between nodes. Therefore, the spatial relationship between hydrological stations is stored using an adjacency matrix, which is defined as shown in Formula 1.

[0014]

[0015] Where a ij represents the element in the i-th row and j-th column of the adjacency matrix, where i and j represent rainfall stations. ij represents the geographical distance (Euclidean distance) between rainfall stations i and j, D ij represents the river distance between rainfall stations i and j. If there is no upstream and downstream relationship between rainfall stations i and j, then a ij is 0, α, β are hyper parameters.

[0016] Step S1.2: Based on the geographical location relationship diagram of the hydrological stations obtained in step S1.1, the hydrological data of a certain station is integrated with the data of the surrounding related hydrological stations, as shown in Formula 2.

[0017] X i+1 =AX i #(Formula 2)

[0018] Where A represents the adjacency matrix obtained in the previous step, X i+1 Represents the fusion matrix of the i+1th layer.

[0019] Furthermore, the detailed process of input data alignment in step S2 is as follows:

[0020] Step S2.1, ArcGIS was used to mark the station locations from the U-Tokyo MERIT DEM digital elevation data, and to extract spatial information such as elevation, slope, and river channel length. The flow velocity was estimated and the confluence time was calculated according to the Manning formula, as shown in Equation 3.

[0021]

[0022] Where n is the surface roughness coefficient, R is the hydraulic radius, and S is the water surface gradient of the river section.

[0023] Based on the grid flow velocity calculated by the above formula and the grid size, the residence time of water flow in each grid can be obtained. When the water flow direction is parallel to one side of the grid, the calculation method is shown in Formula 4.

[0024]

[0025] When the water flow direction is parallel to the grid diagonal, the calculation method is shown in Equation 5.

[0026]

[0027] Finally, the residence time in each grid on the water flow path is added together to obtain the final confluence time.

[0028]

[0029] Where m is the number of grids that the confluence path passes through, ΔT i is the residence time of water flow at grid i.

[0030] Step S2.2: According to the confluence time calculated in step S2.1, the input data is shifted in the time dimension. Assuming that the input data has a sliding window width of 10: (X t-10 , X t-9 , X t-8 , X t-7 …X t-1 ), where X t is the time series data at time t. If the estimated confluence time is 4 hours, the aligned input data is (X t-14 , X t-13 , X t-12 , X t-11 …X t-5 ), the confluence time of each station is different, so the offset of different features is also different.

[0031] Furthermore, the detailed process of using variational autoencoder VAE to train the model in step S3 is as follows:

[0032] Step S3.1: Input the time series data of multiple hydrological stations into the network structure model. The network model is generally a variational autoencoder (VAE) structure, in which the encoder uses a time domain convolutional network (TCN) and a graph convolution neural network (GCN) to extract spatiotemporal features.

[0033] Step S3.2: The time series features X1 and spatial features X2 extracted by the encoder are fused with the original time series data X3 in the feature dimension, and then the latent vector Z is obtained through nonlinear transformation.

[0034] Z=f(X1|X2|X3)#(Formula 7)

[0035] Step S3.3: Input the latent vector Z into the decoder and prediction module to obtain generated data and predicted data.

[0036] X1=f1(Z)#(Formula 8)

[0037] X2=f2(Z)#(Formula 9)

[0038] In Formula 8, f1 is the decoding operation, X1 is the generated data, and in Formula 9, f2 is the prediction model, and X2 is the predicted data.

[0039] Furthermore, the detailed process of calculating the abnormality score in step S4 is as follows:

[0040] The test set data X = (x0, x1, x2...x T) is input into the pre-trained model in step S3, and the model output value and the true value are used to calculate the anomaly score Score. The calculation formula is shown in Formula 10.

[0041]

[0042] Where γ is a hyperparameter, x i represents the true value at the i-th moment, and are the estimated values ​​of the real data, and represent the output values ​​of the prediction model and the generation model at the i-th moment respectively.

[0043] Furthermore, the detailed process of step S5 is as follows:

[0044] Step S5.1: First, calculate the initial threshold γ using the training set data; then, maximize the F1 Score on the validation set by continuously increasing the threshold, as shown in Formula 11.

[0045] max f1(L, score, γ)#(Formula 11)

[0046] Where L is the real data label. Score is the anomaly score calculated using the anomaly score calculation formula in step S4. The threshold corresponding to the maximum F1 Score is the final threshold.

[0047] Step S5.2: Compare the final threshold obtained in the above step S5.1 with the anomaly score calculated in step S4, as shown in formula 12. If the anomaly score is greater than the threshold, it is judged as abnormal data.

[0048]

[0049] Where Y j is the label at time j, 1 represents abnormal data, 0 represents normal data, Score j is the anomaly score at time j, and γ is the threshold.

[0050] Beneficial effects: The present invention uses DEM data to mine the spatiotemporal characteristics of the watershed, and uses an intelligent optimization algorithm to conduct anomaly detection research on hydrological time series in combination with spatiotemporal characteristics. The application scenario of the present invention is anomaly detection of hydrological time series. For the existing DEM data and rainfall data of the watershed, the DEM data is first used to estimate the confluence time, and the input data in the sliding window is aligned according to the confluence time to extract the time characteristics; then, according to the river distance and the Euclidean distance between the stations, a spatial graph between the stations is constructed and used in the graph convolutional neural network; then the aligned input data and the spatial graph structure are used for neural network training, and the differential evolution algorithm is used to select the optimal hyperparameters; finally, the anomaly score is calculated based on the reconstruction error and prediction error of the model, and compared with the threshold to obtain the abnormal data. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 It is a flow chart in an embodiment of the present invention.

[0052] Figure 2 is a spatial relationship matrix of a watershed in an embodiment of the present invention;

[0053] Figure 3 It is a process of aligning time series data of a certain river basin in an embodiment of the present invention;

[0054] Figure 4 is a neural network structure diagram used in the embodiment of the present invention

[0055] Figure 5 1 is a diagram of the model training results and anomaly detection results in an embodiment of the present invention. DETAILED DESCRIPTION

[0056] The technical solution of the present invention is described in detail below, but the protection scope of the present invention is not limited to the embodiments.

[0057] like Figure 1 As shown, the similar rainfall process search method based on deep learning in this embodiment includes the following steps:

[0058] Step S1, extract the elevation, slope, and geographic distance between stations and other information of stations in the basin from the digital elevation model (DEM) data. And use the graph structure G = (V, E) to represent the rainfall stations in the basin and the relationship between stations, where the vector V represents the rainfall stations in the basin; E represents the edge, and the weight of the edge is the degree of association between the rainfall stations.

[0059] Step S2: First, the grid flow velocity is calculated using the Manning formula according to the river distance, river width and elevation information; then the cumulative confluence time is calculated according to the residence time of the rainfall in the grid; finally, the time series is aligned in the time dimension according to the confluence time from each rainfall station to the station to be tested.

[0060] Step S3: Use an improved variational autoencoder (VAE) to build a network model, where the encoder uses a graph convolutional neural network (GCN) to learn spatial features and a time domain convolutional network (TCN) to learn time series features, and the optimal network model is obtained after training.

[0061] Step S4: Input the test set data into the model trained in step S3 to obtain the predicted value and generated value of the real data, and calculate the anomaly score based on the difference between the predicted value, generated value and the real value.

[0062] Step S5, design a method to maximize the F1 Score to calculate the threshold. First, calculate the anomaly score using the training set data to obtain the initial threshold. Then, iteratively increase the threshold on the validation set to maximize the F1 Score. The threshold corresponding to the maximum F1 Score is the final threshold. Finally, compare the anomaly score in step S4 with the threshold. If the anomaly score is greater than the threshold, it is judged as abnormal data.

[0063] In step S1 of the embodiment of the present invention, ArcGIS is used to mark the locations of each rainfall station and Dongshan water level station of Qinhuai River from U-Tokyo MERIT DEM digital elevation data, and calculate the Euclidean distance between the stations, as well as spatial information such as river length. Due to the spatial autocorrelation between the stations, the Euclidean distance between the stations and the reachable distance to the river can be used to construct a spatial relationship diagram between the stations. In the graph structure, the hydrological stations are represented by graph nodes, and the relationships between the nodes are stored by edges. Therefore, the spatial relationship between the hydrological stations is stored by an adjacency matrix, which is defined as shown in Formula 1.

[0064]

[0065] Where a ij represents the element in the i-th row and j-th column of the adjacency matrix, where i and j represent rainfall stations. ij represents the geographical distance (Euclidean distance) between rainfall stations i and j, D ij represents the river distance between rainfall stations i and j. If there is no upstream and downstream relationship between rainfall stations i and j, then a ij is 0, α and β are hyperparameters. The adjacency matrix of the spatial relationship graph of the Qinhuai River Basin is as follows: Figure 2 And use the adjacency matrix to merge the surrounding site data.

[0066] X i+1 =AX i #(Formula 2)

[0067] Where A represents the adjacency matrix obtained in the previous step, X i+1 Represents the fusion matrix of the i+1th layer.

[0068] In step S2 of the embodiment of the present invention, spatial information such as elevation, slope, and river length are extracted from the GDEMV2 30M resolution digital elevation data, and the flow velocity is estimated and the confluence time is calculated according to the Manning formula, as shown in Formula 3.

[0069]

[0070] Where n is the surface roughness coefficient, R is the hydraulic radius, and S is the water surface gradient of the river section. Based on the grid flow velocity calculated by the above formula and the grid size, the retention time of water flow in each grid can be obtained. When the water flow direction is parallel to one side of the grid, the calculation method is shown in Formula 4.

[0071]

[0072] When the water flow direction is parallel to the grid diagonal, the calculation method is shown in Equation 5.

[0073]

[0074] Finally, the residence time in each grid on the water flow path is added together to obtain the final confluence time.

[0075]

[0076] Where m is the number of grids that the confluence path passes through, ΔT i is the residence time of water flow at grid i.

[0077] In step S2 of the embodiment of the present invention, the input data is aligned in time dimension according to the confluence time. The schematic diagram of the aligned data is as follows: Figure 3 Assume that the input data of a sliding window width of 10: (X t-10 , X t-9 , X t-8 , X t-7 …X t-1 ), where X t is the time series data at time t. If the estimated confluence time is 4 hours, the aligned input data is (X t-14 , X t-13 , X t-12 , X t-11 …X t-5), the confluence time of each station is different, so the offset of different features is also different.

[0078] In step S3 of the embodiment of the present invention, the time series data of multiple hydrological stations are input into the network structure model. The network model is generally a variational autoencoder (VAE) structure, in which the encoder uses a time domain convolutional network (TCN) and a graph convolution neural network (GCN) to extract spatiotemporal features.

[0079] In the embodiment of the present invention, the time series features X1 and the spatial features X2 extracted by the encoder are fused with the original time series data X3 in the feature dimension, and then a latent vector Z is obtained through nonlinear transformation.

[0080] Z=f(X1|X2|X3)#(Formula 7)

[0081] And input the hidden vector Z into the decoder and prediction module to obtain generated data and predicted data.

[0082] X1=f1(Z)#(Formula 8)

[0083] X2=f2(Z)#(Formula 9)

[0084] In Formula 8, f1 is the decoding operation, X1 is the generated data, and in Formula 9, f2 is the prediction model, and X2 is the predicted data.

[0085] In step S4 of the embodiment of the present invention, the test set data X=(x0, x1, x2...x T ) is input into the pre-trained model in step S3, and the model output value and the true value are used to calculate the anomaly score Score. The calculation formula is shown in Formula 10.

[0086]

[0087] Where γ is a hyperparameter, x i represents the true value at the i-th moment, and are the estimated values ​​of the real data, and represent the output values ​​of the prediction model and the generation model at the i-th moment respectively.

[0088] In step S5 of the embodiment of the present invention, the initial threshold γ is first calculated using the training set data; then the F1 Score is maximized on the validation set by continuously increasing the threshold, as shown in Formula 11.

[0089] max f1(L, score, γ)#(Formula 11)

[0090] Where L is the real data label. Score is the anomaly score calculated using the anomaly score calculation formula in step S4. The threshold corresponding to the maximum F1 Score is the final threshold.

[0091] In the embodiment of the present invention, the final threshold obtained in step S5 is compared with the abnormal score calculated in step S4, as shown in formula 12. If the abnormal score is greater than the threshold, it is judged as abnormal data.

[0092]

[0093] Where Y j is the label at time j, 1 represents abnormal data, 0 represents normal data, Score j is the anomaly score at time j, and γ is the threshold.

Claims

1. A hydrological time series anomaly detection method based on spatiotemporal characteristics, characterized by: The following steps are involved: Step S1, extract the elevation, slope, and geographic distance between stations and other information of stations in the basin from the digital elevation model (DEM) data. And use the graph structure G = (V, E) to represent the rainfall stations in the basin and the relationship between the stations, where the vector y represents the rainfall stations in the basin; E represents the edge, and the weight of the edge is the degree of association between the rainfall stations. Step S2: First, the grid flow velocity is calculated using the Manning formula according to the river distance, river width and elevation information; then the cumulative confluence time is calculated according to the residence time of the rainfall in the grid; finally, the time series is aligned in the time dimension according to the confluence time from each rainfall station to the station to be tested. Step S3: Use an improved variational autoencoder (VAE) to build a network model, where the encoder uses a graph convolutional neural network (GCN) to learn spatial features and a time domain convolutional network (TCN) to learn time series features, and the optimal network model is obtained after training. Step S4: Input the test set data into the model trained in step S3 to obtain the predicted value and generated value of the real data, and calculate the anomaly score based on the difference between the predicted value, generated value and the real value. Step S5, design a method to maximize the F1 Score to calculate the threshold. First, calculate the anomaly score through the training set data to obtain the initial threshold; then maximize the F1 Score on the validation set by iteratively increasing the threshold, and the threshold corresponding to the maximum F1 Score is the final threshold; finally, compare the anomaly score in step S4 with the threshold. If the anomaly score is greater than the threshold, it is judged as abnormal data.

2. The hydrological time series anomaly detection method based on spatiotemporal characteristics according to claim 1 is characterized by: The detailed process of step S1 is as follows: Step S1.1: Use ArcGIS to mark the station locations from the U-Tokyo MERIT DEM digital elevation data and calculate the Euclidean distance between stations and spatial information such as river length. Since there is spatial autocorrelation between stations, the Euclidean distance between stations and the reachable distance to the river can be used to construct a spatial relationship graph between stations. In the graph structure, graph nodes are used to represent hydrological stations, and edges are used to store the relationship between nodes. Therefore, the spatial relationship between hydrological stations is stored using an adjacency matrix, which is defined as shown in Formula 1. Where a ij represents the element in the i-th row and j-th column of the adjacency matrix, where i and j represent rainfall stations. ij represents the geographical distance (Euclidean distance) between rainfall stations i and j, D ij represents the river distance between rainfall stations i and j. If there is no upstream and downstream relationship between rainfall stations i and j, then a ij is 0, α, β are hyper parameters. Step S1.2: Based on the geographical location relationship diagram of the hydrological stations obtained in step S1.1, the hydrological data of a certain station is integrated with the data of the surrounding related hydrological stations, as shown in Formula 2. X i+1 =AX i #(Formula 2) Where A represents the adjacency matrix obtained in the previous step, X i+1 Represents the fusion matrix of the i+1th layer.

3. The hydrological time series anomaly detection method based on spatiotemporal characteristics according to claim 1 is characterized by: The detailed process of step S2 is as follows: Step S2.1, ArcGIS was used to mark the station locations from the U-Tokyo MERIT DEM digital elevation data, and to extract spatial information such as elevation, slope, and river channel length. The flow velocity was estimated and the confluence time was calculated according to the Manning formula, as shown in Equation 3. Where n is the surface roughness coefficient, R is the hydraulic radius, and S is the water surface gradient of the river section. Based on the grid flow velocity calculated by the above formula and the grid size, the residence time of water flow in each grid can be obtained. When the water flow direction is parallel to one side of the grid, the calculation method is shown in Formula 4. When the water flow direction is parallel to the grid diagonal, the calculation method is shown in Equation 5. Finally, the residence time in each grid on the water flow path is added to obtain the final confluence time, which is calculated as shown in Equation 6. Where m is the number of grids that the confluence path passes through, ΔT i is the residence time of water flow at grid i. Step S2.2: According to the confluence time calculated in step S2.1, the input data is shifted in the time dimension. Assuming that the input data has a sliding window width of 10 Where X t is the time series data at time t. If the estimated confluence time is 4 hours, the aligned input data is (X t-14 , X t-13 ,X t-12 ,X t-11 …X t-5 ), the confluence time of each station is different, so the offset of different features is also different.

4. The hydrological time series anomaly detection method based on spatiotemporal characteristics according to claim 1 is characterized by: The process of model training in step S3 is as follows: Step S3.1: Input the time series data of multiple hydrological stations into the network structure model. The network model is generally a variational autoencoder (VAE) structure, in which the encoder uses a time domain convolutional network (TCN) and a graph convolution neural network (GCN) to extract spatiotemporal features. Step S3.2: The time series features X1 and spatial features X2 extracted by the encoder are fused with the original time series data X3 in the feature dimension, and then the latent vector Z is obtained through nonlinear transformation. Z=f(X1|X2|X3)#(Formula 7) Step S3.3: Input the latent vector Z into the decoder and prediction module to obtain generated data and predicted data. X1=f1(Z)#(Formula 8) X2=f2(Z)#(Formula 9) In Formula 8, f1 is the decoding operation, X1 is the generated data, and in Formula 9, f2 is the prediction model, and X2 is the predicted data.

5. The hydrological time series anomaly detection method based on spatiotemporal characteristics according to claim 1 is characterized by: The detailed process of calculating the anomaly score in step S4 is as follows: The test set data X = (x0, x1, x2...x T ) is input into the pre-trained model in step S3, and the model output value and the true value are used to calculate the anomaly score Score. The calculation formula is shown in Formula 10. Where γ is a hyperparameter, x i represents the true value at the i-th moment, and are the estimated values ​​of the real data, and represent the output values ​​of the prediction model and the generation model at the i-th moment respectively.

6. The hydrological time series anomaly detection method based on spatiotemporal characteristics according to claim 1 is characterized by: The detailed process of threshold determination and abnormality determination in step S5 is as follows: Step S5.1: First, calculate the initial threshold y using the training set data; then, maximize the F1 Score on the validation set by continuously increasing the threshold, as shown in Formula 11. max f1(L,score,γ)#(Equation 11) Where L is the real data label. Score is the anomaly score calculated using the anomaly score calculation formula in step S4. The threshold corresponding to the maximum F1 Score is the final threshold. Step S5.2: Compare the final threshold obtained in the above step S5.1 with the anomaly score calculated in step S4, as shown in formula 12. If the anomaly score is greater than the threshold, it is judged as abnormal data. Where Y j is the label at time j, 1 represents abnormal data, 0 represents normal data, Score j is the anomaly score at time j, and γ is the threshold.

Images