A cleaning method for large-scale long-duration abnormal flood data

Abstract

This invention discloses a method for cleaning large-scale, long-duration anomalous flood reporting data. The method establishes an anomaly detection model for each observation station and trains it to identify anomalous rainfall. It also establishes a rainfall data interpolation model and trains it to interpolate anomalous or missing rainfall data. All rainfall monitoring data are detected using the anomaly detection model to identify anomalous data. A sample time period is defined, and a clustering algorithm is used to perform cluster analysis on all anomalous data within that time period to identify anomalous events. The rainfall data interpolation model then interpolates the anomalous or missing data related to these events. This method offers higher identification efficiency and accuracy for large-scale, long-duration anomalous flood reporting data, possesses better anomalous data cleaning capabilities, and can provide high-quality data support for hydrological forecasting, reservoir scheduling, and smart water conservancy construction.

Description

Technical Field

[0001] This invention belongs to the field of hydrological observation data processing technology, specifically relating to a method for cleaning large-scale, long-duration abnormal flood reporting data. Background Technology

[0002] Flood forecasting data, comprised of observational data including various hydrological parameters, serves as a crucial basis for flood warnings, reservoir scheduling, and water resource management. Its quality directly impacts the rationality of water resource allocation and the guarantee of public safety. In recent years, with the development of smart water conservancy projects in my country, the number of hydrological observation stations has increased rapidly. These stations are generating massive amounts of flood forecasting data at an unprecedented density and frequency. However, various interference factors, such as sensor malfunctions and communication interruptions, can affect the quality of flood forecasting data. This impact is amplified rapidly with the increase in data volume, frequently leading to false or inaccurate flood reports. Especially under extreme weather conditions, the probability of long-duration, large-scale abnormal flood forecasting data increases significantly, resulting in a severe decline in the accuracy of flood forecasts.

[0003] To address the issue of anomalies in flood reporting data, it is necessary to clean the data before analysis. To improve data cleaning efficiency, researchers have proposed an integrated cleaning method based on machine learning. For example, Chinese patent ZL202510487396.5, entitled "Method and System for Cleaning Hydrological Data of Hydrological Stations Based on Big Data Analysis," proposes a method that trains a neural network using historical hydrological data to obtain a prediction model. Then, it compares the model's prediction with the observed values ​​to determine the authenticity of the hydrological station's observation data and performs data cleaning based on the judgment results. However, when processing long-sequence data, neural networks are prone to gradient vanishing and feature decay problems, making it impossible to effectively learn data features. In addition, traditional neural networks mainly extract shallow features, making it difficult to capture deep nonlinear variation patterns in rainfall time series data, thus failing to guarantee the reliability of the prediction. For example, a Chinese patent application, titled "A Method for Cleaning Hydrological Remote Sensing Rainfall Data Based on Machine Learning," application number 202511156702.3, proposes a method that analyzes the spatiotemporal correlation of observation stations using the XGBoost algorithm, identifies anomalous data based on the obtained correlations, and then corrects the anomalous data using the Random Forest algorithm. However, this method is only suitable for cleaning single-site data. When large-scale, long-duration anomalous events occur, the missing data is highly concentrated in time and space, and there is a severe shortage of available reference stations in the surrounding area, resulting in limited cleaning effectiveness of this method.

[0004] In summary, existing hydrological data cleaning methods lack the ability to effectively handle large-scale, long-duration anomalous events. This is because traditional neural networks struggle to effectively learn the features of long-sequence data, and once a large-scale, long-duration anomalous event occurs, data is largely missing in both time and space, making data completion difficult. Missing data needs to be supplemented by integrating multi-source hydrological or meteorological data. Therefore, proposing an integrated data cleaning method that can effectively process long-sequence data and integrate multi-source hydrological or meteorological data for supplementation is of great significance for improving the accuracy of flood warnings. Summary of the Invention

[0005] The purpose of this invention is to provide a cleaning method for large-scale, long-duration anomalous flood reporting data. This method first uses an autoencoder model to detect anomalies in single-station time-series data, then uses a density clustering algorithm (DBSCAN) to identify the areas covered by large-scale regional anomalous events, and finally uses a self-attention model (Transformer) combined with multi-source hydrological or meteorological data to replace or fill in the anomalous areas.

[0006] The technical solution adopted in this invention is as follows: A method for cleaning large-scale, long-duration abnormal flood reporting data is disclosed. This method operates based on a hydrological monitoring system comprising multiple hydrological observation stations, each with a unique number, and the locations of any two observation stations are distinct. The method specifically includes the following steps: S1. Set a sampling period, collect meteorological monitoring data measured by all observation stations in the system during the sampling period, determine the location of each observation station, and obtain meteorological data related to the system's monitoring range during the sampling period through another means.

[0007] S2. Establish multiple anomaly detection models, which are neural network models. The number of anomaly detection models is equal to the number of observation stations and they correspond one-to-one. All anomaly detection models are trained using meteorological monitoring data from the corresponding observation stations, so that all anomaly detection models have the ability to identify abnormal rainfall.

[0008] S3. Establish a rainfall data interpolation model, which is a self-attention model. The rainfall data interpolation model is trained using meteorological monitoring data and meteorological related data obtained through another means, so that it has the ability to predict rainfall distribution.

[0009] S4. During the monitoring process, all monitoring data measured at each observation station at all times are detected by their corresponding anomaly discrimination model to determine whether the monitoring data is abnormal.

[0010] S5. Set a sample time period, and use a clustering algorithm to perform cluster analysis on all abnormal data that appear within the sample time period to obtain one or more abnormal events.

[0011] S6. Obtain the rainfall distribution within the spatiotemporal range involved in the abnormal event through the rainfall data interpolation model, and replace or fill the abnormal data based on the obtained rainfall distribution.

[0012] Because large-scale, long-duration anomalous events can cause anomalies or missing data in rainfall monitoring data over a wide area, traditional methods of interpolation using data from nearby observation stations are insufficient for replacing anomalous data and filling in missing data. Therefore, this invention proposes a rainfall data interpolation model for predicting rainfall based on meteorological data obtained through another means. This model is used to replace and fill in the missing data. Furthermore, the anomalies in large-scale rainfall monitoring data are themselves characteristic of anomalous events. Each cluster obtained through clustering algorithms corresponds to an anomalous event. By observing the time-varying characteristics of the shape and size of a cluster, information such as the type and intensity of the anomalous event corresponding to that cluster can be obtained, which helps in the accurate prediction of flood conditions.

[0013] Further optimization involves using an autoencoder model as the anomaly detection model in step S2. The process of training the autoencoder model specifically includes the following steps: S2.1. Arrange all rainfall data in the meteorological monitoring data in chronological order to form an original sequence. Set a window length to establish a sliding window. Divide the original sequence into multiple equal-length subsequences through the sliding window. Then, divide the rainfall process corresponding to the original sequence into a first stage, a second stage, and a third stage. Subsequences corresponding to the first stage are assigned to the training set, and subsequences corresponding to the second and third stages are assigned to the validation set and the test set, respectively. Construct the first, second, and third sample matrices based on the training set, validation set, and test set, respectively. Each subsequence is a row vector of the sample matrix, and all subsequences are arranged in chronological order in the sample matrix.

[0014] S2.2. Set initial weights and initial weight learning rates, establish an initial AutoEncoder model, input the first sample matrix into the initial AutoEncoder model, and perform multiple rounds of training in combination with other meteorological monitoring data from the first stage. After each round of training, the learning rates of the weights of all hidden layers in the AutoEncoder model are updated by the Adam optimizer. The mean squared error function is used as the first loss function. After each round of training, the obtained model is used to predict the rainfall at all times in the first and second stages based on meteorological monitoring data other than rainfall. Then, the training loss is obtained based on the predicted values ​​of the first stage and the corresponding values ​​in the training set using the first loss function, and the validation loss is obtained based on the predicted values ​​of the second stage and the corresponding values ​​in the validation set.

[0015] S2.3. Set termination conditions, the termination conditions including a threshold number of attempts. M With the loss threshold, when it appears continuously during training M In each round, the training of the model is considered to have met the termination condition when the difference between the validation losses of any two adjacent rounds is less than the loss threshold, and the training stops. Find the validation loss with the smallest value and use the model with the smallest validation loss as the sample model.

[0016] S2.4. Establish a comprehensive quality scoring model for evaluating the performance of the sample model, and set a scoring threshold. The comprehensive quality scoring model is expressed as the following formula: (1) in, i Number the observation station, Quality i For observation station i Scoring of the sample model L i train For training loss, L i val To verify the loss, N i The number of rounds in the training process. s i MSE The standard deviation of the error is calculated based on the second-stage predicted values ​​and the corresponding values ​​in the validation set. L To maximize the loss, max N The maximum number of rounds, max s For the maximum standard deviation, W 1. W 2. W 3. W 4 are all weighting coefficients.

[0017] The quality of the sample model is calculated using a comprehensive quality scoring model. i If Quality i If the score is not lower than the set scoring threshold, the sample model is considered to have met the performance standard and will be used as the final anomaly detection model. i If the score is below the threshold, the sample model is considered to have met the performance standard, and the model is retrained according to the procedures described in steps S2.2 and S2.3.

[0018] S2.5. The obtained anomaly discrimination model is used to predict the rainfall at all times in the third stage, and the reconstruction error at each time is obtained based on the predicted value of the third stage and the corresponding value in the test set. Then, the anomaly threshold is determined according to the overall distribution of the reconstruction error. When the detection is performed, the anomaly discrimination model obtains the predicted rainfall value at a certain time through the meteorological monitoring data at a certain time, and compares the rainfall monitoring data and the predicted value. If the error between the two is greater than the anomaly threshold, the monitoring data is considered to be abnormal data.

[0019] The AutoEncoder model employs a symmetrical encoder-decoder structure, learning the temporal features of rainfall through a compression-reconstruction process. The encoder maps the input to a low-dimensional latent space through layer-by-layer compression, while the decoder symmetrically reconstructs the latent space representation back to the original dimension. Each hidden layer is followed by a batch normalization and dropout layer, exhibiting strong generalization ability. The termination condition prevents overfitting. A comprehensive quality scoring model is used to evaluate the reliability of the anomaly detection model corresponding to each observation station from four aspects: training loss, validation loss, number of epochs, and error standard deviation.

[0020] Further optimization involves using a Transformer model for rainfall data interpolation in step S3, and ERA5 reanalysis data obtained through another method. Step S3 specifically includes the following steps: S3.1. Establish a grid covering the overall monitoring range of the hydrological monitoring system, with each grid node corresponding to a unique serial number, and determine the grid unit where each observation station is located; based on the rainfall data measured by each observation station at a certain moment, combined with the grid unit where the observation station is located and the location of the grid node, the rainfall at all grid nodes at that moment is obtained by inverse distance weighted interpolation (IDW), and the rainfall at each grid node at all moments during the sampling period is obtained by following the above process.

[0021] S3.2. Divide the sampling period into a training phase, a verification phase, and a testing phase. Sort the rainfall at all times at the same grid node by time to obtain a sample sequence. Record the part of the sample sequence corresponding to the training phase as the first sequence, the part corresponding to the verification phase as the second sequence, and the part corresponding to the testing phase as the third sequence. Each grid node obtains its corresponding first, second, and third sequences according to the above process.

[0022] S3.3. The ERA5 reanalysis data includes four parameters: convective effective potential energy, temperature, total column water vapor, and total precipitation. The ERA5 reanalysis data is preprocessed, and then bilinear interpolation is used to obtain the ERA5 reanalysis data for each grid node; a time-domain window is set […]. t - t , t + t ], t For the current moment, t The time-domain radius and t ∈Z + A feature matrix is ​​established for each grid node at all times within the sampling period. The feature matrix contains multiple feature vectors and a static vector. The number of feature vectors is equal to and corresponds one-to-one with the number of times contained in the time-domain window. The feature vector contains four elements, which correspond to the four parameters contained in the ERA5 reanalysis data. The static vector contains the corresponding grid node position parameter. The time corresponding to a certain feature matrix is ​​taken as the current time of the time-domain window to determine the time corresponding to the feature vector of the matrix. The value of the feature vector element is determined by combining the time and the corresponding grid node position.

[0023] S3.4. Set the initial learning rate, weight decay coefficient, and learning rate update method, and construct the AdamW optimizer; establish the Transformer model and the second loss function, and train the Transformer model in multiple rounds by combining the first sequence of all grid nodes and the feature matrices of all time points in the training phase. After each round of training, the rainfall at all time points in the validation phase is predicted using the feature matrix of the obtained model set. Then, the validation loss is obtained based on the predicted values ​​in the validation phase and the corresponding time values ​​in the second sequence using the second loss function. After each round of training, the values ​​of all parameters to be learned in the Transformer model are updated by the AdamW optimizer based on the validation loss; set the round threshold. N When training continues N Training stops when the validation loss in the current round does not decrease compared to the previous round. Then, the validation loss with the minimum value is found, and the model with the minimum validation loss is used as the imputation model for the obtained rainfall data.

[0024] S3.5. Based on the feature matrix of the test phase, the obtained rainfall data interpolation model is used to predict the rainfall at all times during the test phase. The prediction ability of the obtained rainfall data interpolation model under various weather conditions is evaluated based on the predicted values ​​of the test phase and the corresponding values ​​in the third sequence.

[0025] ERA5 (The fifth generation ECMWF Re-Analyses) is global atmospheric reanalysis data generated by the European Centre for Medium-Range Weather Forecasts (ECMWF) through the assimilation of multi-source observational data. It contains a large number of climate variables and is widely used in the field of meteorological analysis. After preprocessing the meteorological data contained in the collected ERA5 reanalysis data, including time zone unification, standardization, and elimination of dimensional differences, the distribution of rainfall can be obtained from the relevant data and presented in a grid-like manner.

[0026] Further optimization: The second loss function established in step S3.4 includes two terms: mean squared error (MSE) and mean absolute error (MAE). The second loss function is expressed as follows: (2) in, k For the number of training rounds, j For grid node number, L k for k Validation loss after each training round T To verify the duration of the phase, w j For grid nodes j The sample weights at each location, under rainy and dry weather conditions. w j The values ​​are different. To pass k The model obtained from 1 round of training t Time grid nodes j Rainfall forecast value, For grid nodes j In the second sequence t Rainfall at any given time l MSE The mean square error weights are... l MAE This is the weight for the average absolute error.

[0027] The sample weight values ​​are determined by the rainfall at the corresponding grid node. The sample weight values ​​are higher in rainy weather than in dry weather, which increases the impact of rainfall data on model training and improves the accuracy of model prediction. The specific values ​​of the sample weights can be adaptively set according to the rainfall.

[0028] Further optimization involves using a clustering algorithm based on the DBSCAN algorithm in step S5. Step S5 specifically includes the following steps: S5.1. Establish sample points belonging to each moment within the sample time period. The number of sample points belonging to a certain moment is equal to the number of abnormal data appearing at that moment and they correspond one-to-one. The location of any sample point is the same as the location of its corresponding abnormal data. The location of any abnormal data is the location of the observation station from which the abnormal data was obtained. All sample points are grouped into the same set for cluster analysis.

[0029] S5.2. Set spatial and temporal thresholds. When a sample point is used as the center point, if the straight-line distance between another sample point and the center point does not exceed the set spatial threshold, and the time difference between the other sample point and the center point does not exceed the set temporal threshold, then the sample point is considered to be within the neighborhood radius of the center point; otherwise, it is considered not to be within the neighborhood radius. Set a minimum number of points, and perform cluster analysis on all sample points using the DBSCAN algorithm to obtain one or more clusters. Determine the corresponding abnormal events based on the morphological characteristics and temporal and spatial scales of the clusters.

[0030] The characteristics of anomalous events include the spatial range involved and the duration of time. The traditional DBSCAN algorithm only considers the spatial distance between sample points. In order to achieve effective identification of anomalous events, this invention sets both spatial and temporal thresholds to accurately describe the dynamic characteristics of anomalous events.

[0031] Further optimization, step S6 specifically includes the following steps: S6.1. The sample points to be assigned to the cluster are recorded as points to be replaced. The grid cell where each point to be replaced is located is determined according to the location of the points to be replaced and the grid distribution. The grid nodes contained in the grid cell where a point to be replaced is located are recorded as the associated nodes of that point to be replaced.

[0032] S6.2. Determine the associated nodes at all times within the sample time period. Based on the ERA5 reanalysis data, establish the feature matrix of each associated node by combining the location and time of occurrence of the associated nodes. Input all feature matrices into the rainfall data interpolation model to obtain the rainfall prediction value for each associated node.

[0033] S6.3. Using bilinear interpolation, combined with the location of the point to be swapped and its associated nodes, as well as the predicted rainfall values ​​of the associated nodes, new rainfall data for the point to be swapped is obtained. The new data is then used to replace the abnormal data corresponding to the point to be swapped.

[0034] The beneficial effects of the method of the present invention are as follows: 1. The method of the present invention identifies individual monitoring data of a single station through an anomaly discrimination model, which improves the efficiency and accuracy of anomaly data identification. At the same time, it uses ERA5 reanalysis data as input features of rainfall data interpolation model to replace or fill in anomaly data, which has better cleaning ability for monitoring data of large-scale long-term anomaly events. 2. The method of this invention includes a DBSCAN-based clustering algorithm, which realizes feature identification of abnormal events through continuous time and scale. It can identify multi-scale abnormal events, which is helpful for judging the type and scope of impact of abnormal events and can help with flood forecasting. 3. The training process of the anomaly detection model includes the determination of termination conditions and the judgment of the comprehensive quality scoring model, which improves the accuracy of the anomaly detection model in identifying abnormal data and increases the reliability of the method of the present invention. 4. The second loss function includes sample weights with different values ​​for rainy and dry weather, which helps to improve the quality of the data provided by the rainfall data interpolation model. Attached Figure Description

[0035] Figure 1 A schematic diagram of the overall process of the data cleaning method of this invention.

[0036] Figure 2. Schematic diagram of the anomaly detection effect of the anomaly discrimination model corresponding to observation station 2993.

[0037] Figure 3. Schematic diagram of the anomaly detection effect of monitoring data from multiple observation stations within the same time period.

[0038] Figure 4. Schematic diagram of classification results obtained by clustering algorithm under different settings.

[0039] Figure 5 Scatter plot comparing actual rainfall at grid nodes with Transformer model predictions.

[0040] Figure 6. Spatial distribution comparison of rainfall on the day with the best and worst prediction results of the rainfall data interpolation model.

[0041] Figure 7. Comparison of predicted and actual rainfall at grid nodes. Detailed Implementation

[0042] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below through specific embodiments. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] Example 1: A method for cleaning large-scale, long-duration abnormal flood reporting data. This method operates based on a hydrological monitoring system, which includes multiple hydrological observation stations, each with a unique number, and the locations of any two observation stations are different. The hydrological monitoring system used in this example includes 11,805 hydrological observation stations, and the system's monitoring range covers the area located in the upper reaches of the Yangtze River, between 91°E and 112°E, and between 24°N and 37°N. The specific steps include: S1. Set a sampling period, collect meteorological monitoring data measured by all observation stations in the system during the sampling period, and determine the location of each observation station. In this embodiment, the sampling period is from 00:00 on January 1, 2013 to 23:00 on December 31, 2023. The rainfall data is recorded hourly. The location of the observation station is described by longitude, latitude and elevation. ERA5 reanalysis data is used as another way to obtain meteorological data related to the system's monitoring range during the sampling period.

[0044] S2. Establish multiple anomaly detection models, which are neural network models. The number of anomaly detection models is equal to the number of observation stations and they correspond one-to-one. All anomaly detection models are trained using meteorological monitoring data from the corresponding observation stations, so that all anomaly detection models have the ability to identify abnormal rainfall. In this embodiment, the anomaly detection model is an autoencoder model. The process of training the autoencoder model specifically includes the following steps: S2.1. Arrange all rainfall data in the meteorological monitoring data in chronological order to form an original sequence. Set a window length to establish a sliding window. Divide the original sequence into multiple equal-length subsequences through the sliding window. Then, divide the rainfall process corresponding to the original sequence into a first stage, a second stage, and a third stage. In this embodiment, the window length is 30 h. The first, second, and third stages account for 80%, 10%, and 10% of the total rainfall process duration, respectively. The subsequences corresponding to the first stage are included in the training set, and the subsequences corresponding to the second and third stages are included in the validation set and the test set, respectively. Construct the first, second, and third sample matrices based on the training set, validation set, and test set, respectively. Each subsequence is a row vector of the sample matrix, and all subsequences are arranged in chronological order in the sample matrix.

[0045] S2.2. Set initial weights and initial weight learning rates, establish an initial AutoEncoder model, input the first sample matrix into the initial AutoEncoder model, and perform multiple rounds of training in combination with other meteorological monitoring data from the first stage. After each round of training, the learning rates of the weights of all hidden layers in the AutoEncoder model are updated by the Adam optimizer. The mean squared error function is used as the first loss function. After each round of training, the obtained model is used to predict the rainfall at all times in the first and second stages based on meteorological monitoring data other than rainfall. Then, the training loss is obtained based on the predicted values ​​of the first stage and the corresponding values ​​in the training set using the first loss function, and the validation loss is obtained based on the predicted values ​​of the second stage and the corresponding values ​​in the validation set.

[0046] S2.3. Set termination conditions, the termination conditions including a threshold number of attempts. M In this embodiment, the loss threshold is related to... M =150, loss threshold is 1e-5, when it occurs continuously during training. M In each round, the training of the model is considered to have met the termination condition when the difference between the validation losses of any two adjacent rounds is less than the loss threshold, and the training stops. Find the validation loss with the smallest value and use the model with the smallest validation loss as the sample model.

[0047] S2.4. Establish a comprehensive quality scoring model for evaluating the performance of the sample model, and set a scoring threshold. The comprehensive quality scoring model is shown in formula (1). In this embodiment, the weight coefficients are... W 1 = 0.3, W 2 = 0.3, W 3 = 0.2, W 4 = 0.2.

[0048] The quality of the sample model is calculated using a comprehensive quality scoring model. i If Quality i If the score is not lower than the set scoring threshold, the sample model is considered to have met the performance standard and will be used as the final anomaly detection model. i If the score is below the threshold, the sample model is considered to have met the performance standard, and the model is retrained according to the procedures described in steps S2.2 and S2.3.

[0049] S2.5. The obtained anomaly discrimination model is used to predict the rainfall at all times in the third stage, and the reconstruction error at each time is obtained based on the predicted value of the third stage and the corresponding value in the test set. Then, the anomaly threshold is determined according to the overall distribution of the reconstruction error. When the detection is performed, the anomaly discrimination model obtains the predicted rainfall value at a certain time through the meteorological monitoring data at a certain time, and compares the rainfall monitoring data and the predicted value. If the error between the two is greater than the anomaly threshold, the monitoring data is considered to be abnormal data.

[0050] S3. Establish a rainfall data interpolation model. The rainfall data interpolation model is a self-attention model, trained using meteorological monitoring data and ERA5 reanalysis data to enable it to predict rainfall distribution. In this embodiment, the rainfall data interpolation model is a Transformer model, specifically including the following steps: S3.1. Establish a grid covering the entire monitoring range of the hydrological monitoring system. Each grid node corresponds to a unique serial number, determining the grid unit where each observation station is located. In this embodiment, the grid covering the monitored area has a total of 952,356 nodes. All grid units are the same size and shape, 0.25° × 0.25°. Based on the rainfall data measured by all observation stations at a certain moment, combined with the grid unit where the observation station is located and the location of the grid node, the rainfall at all grid nodes at that moment is obtained by inverse distance weighted interpolation (IDW). The location of the grid node is also described by longitude, latitude, and elevation, and the longitude and latitude are standardized. The rainfall at each grid node at all moments during the sampling period is obtained according to the above process.

[0051] S3.2. The sampling period is divided into a training phase, a verification phase, and a testing phase. In this embodiment, 2013 to 2021 is the training phase, 2022 is the verification phase, and 2023 is the testing phase. The rainfall at all times at the same grid node is sorted by time to obtain a sample sequence. The part of the sample sequence corresponding to the training phase is recorded as the first sequence, the part corresponding to the verification phase is recorded as the second sequence, and the part corresponding to the testing phase is recorded as the third sequence. Each grid node obtains its corresponding first, second, and third sequences according to the above process.

[0052] S3.3. The ERA5 reanalysis data includes four parameters: Convective Available Potential Energy (CAPE), temperature, total column water vapor, and total precipitation. The total precipitation is the hourly cumulative precipitation. Preprocessing of the ERA5 reanalysis data includes unifying the UTC time used in the ERA5 reanalysis data with the Beijing time (UTC+8) used in the meteorological monitoring data, and unifying the rainfall unit to mm. Then, a logarithmic transformation is performed on the total precipitation using the following formula: ,in y orig This represents the original value of total precipitation. y log The total precipitation logarithm was used to standardize the three parameters: CAPE, temperature, and total column water vapor.

[0053] Then, the ERA5 reanalysis data at each grid node was obtained using bilinear interpolation; the time-domain window was set [ t - t , t + t ], t For the current moment, t The time-domain radius and t ∈Z + In this embodiment t = 6. Establish a feature matrix for each grid node at all times within the sampling period. The feature matrix contains multiple feature vectors, and the number of feature vectors is equal to and corresponds one-to-one with the number of times contained in the time-domain window. Each feature vector contains four elements, corresponding to the four parameters contained in the ERA5 reanalysis data. The time corresponding to a certain feature matrix is ​​taken as the current time of the time-domain window to determine the time corresponding to the feature vector of the matrix. The value of the feature vector element is determined by combining the time and the corresponding grid node position. In this embodiment, the feature matrix consists of 13 feature vectors and 1 static vector. The static vector includes the longitude, latitude, and elevation data of the grid node corresponding to the matrix. The 14 vectors in the matrix are numbered by sine and cosine position encoding, with the static vector in the first position, and all subsequent feature vectors arranged in the order of their corresponding times.

[0054] S3.4. Set the initial learning rate, weight decay coefficient, and learning rate update method, and construct the AdamW optimizer. In this embodiment, the initial learning rate is... The weight decay coefficient is 0.01. The learning rate update method is to halve the current learning rate when the validation loss of the current round does not decrease compared with the previous round for 5 consecutive times during the training process. This is specifically implemented through the ReduceLROnPlateau learning rate scheduler. Automatic mixed precision (AMP) is used to accelerate the training process, and the gradient pruning threshold is set to 1.0.

[0055] A Transformer model and a second loss function are established. The Transformer model adopts the Pre-Layer Norm structure, which contains a 6-layer encoder. Each layer contains a self-attention layer and a feedforward neural network. The self-attention layer can perform 8 self-attention calculations simultaneously. The hidden dimension of the feedforward neural network is 1024, and the dropout rate is 0.15. The second loss function is shown in formula (2). When the grid node is in rainy weather, the sample weights are... w j It is 2.0, under rainless weather conditions. w j The value is 1.0, and rain is considered to occur when the hourly rainfall exceeds 0.1 mm. l MSE = 1.0, l MAE = 0.5.

[0056] The Transformer model is trained multiple times by combining the first sequence of all grid nodes and the feature matrices of all time points during the training phase. During training, the four dimensions of the feature vectors are linearly mapped to 256 dimensions through the embedding layer, and the three dimensions of the static vectors are also linearly mapped to 256 dimensions through the embedding layer. After mapping, each of the 14 vectors forms a corresponding label. All labels interact fully through a self-attention mechanism. Then, the label corresponding to the feature vector at the 8th position of the feature matrix is ​​mapped to 1 dimension through the output layer and output as the rainfall prediction value. After each training phase, the rainfall at all time points during the validation phase is predicted using the obtained model set and the feature matrix. Then, the validation loss is derived based on the validation phase prediction value and the corresponding time value in the second sequence using the second loss function. The validation loss is a value directly calculated in the logarithmic domain. After each training phase, the values ​​of all parameters to be learned in the Transformer model are updated based on the validation loss using the AdamW optimizer. A phase threshold is set. N When training continues N Training stops when the validation loss in the current round does not decrease compared to the previous round, in this embodiment. N = 15, and then find the validation loss with the smallest value. The model with the smallest validation loss is used as the imputation model for the obtained rainfall data.

[0057] S3.5. Based on the feature matrix of the test phase, the obtained rainfall data interpolation model predicts rainfall at all times during the test phase. The predictive ability of the obtained rainfall data interpolation model under various weather conditions is evaluated based on the predicted values ​​from the test phase and the corresponding values ​​in the third sequence. The predictive ability of the model is comprehensively evaluated using five parameters: root mean square error (RMSE), mean absolute error (MAE), Nash-Sutcliffe efficiency coefficient (NSE), correlation coefficient (CC), and peak signal-to-noise ratio (PSNR). Before evaluation, the predicted rainfall values ​​need to undergo an inverse logarithmic transformation. The transformation formula is as follows: , of which Log-predicted values The original predicted value, if the obtained If the value is negative, the predicted rainfall will be recorded as 0.

[0058] Based on monitoring data of the region at a certain time in 2022, and combined with the IDW interpolation method, the rainfall at each grid node is calculated. The predicted values ​​at the grid nodes in step S3.4 are compared with the actual rainfall values. Figure 5 As shown in the figure, each data point corresponds to a grid node. The horizontal and vertical axes represent the monitored rainfall data and the predicted value, respectively. The red dashed line is the 1:1 diagonal. The data shown in the figure has an NSE of 0.9579, a CC of 0.9868, an RMSE of only 0.534 mm, a MAE of only 0.242 mm, and a PSNR as high as 33.81 dB. It can be seen from the figure that in the common rainfall range of less than 10 mm, the data points are closely distributed on both sides of the red dashed line, and the model prediction accuracy is high. For rainfall events greater than 10 mm, the overall distribution of data points is slightly lower than the red dashed line, indicating that the model has a certain underestimation trend, but the overall error is controllable.

[0059] Using daily precipitation forecasts as the benchmark, the best forecast day for 2023 was January 18th, and the worst was December 18th. Figures 6(a), 6(b), and 6(c) show the actual rainfall, predicted rainfall, and the difference between actual and predicted rainfall on January 18th, respectively, while Figures 6(d), 6(e), and 6(f) show the actual rainfall, predicted rainfall, and the difference between actual and predicted rainfall on December 18th, respectively. In these six figures, the horizontal axis corresponds to longitude, the vertical axis to latitude, and the color bars on the right correspond to rainfall (unit: mm). The best forecast day had an RMSE of 2.310 mm and a CC of 0.547, successfully reproducing the spatial location and morphology of the heavy rainfall center near 105 °E, with the difference in most areas controlled within ±5 mm. The worst forecast day had an RMSE of 2.703 mm. mm, CC is 0.473. Although the local prediction ability has decreased, the overall spatial structure characteristics have not shown a systematic deviation. The above results show that the rainfall data interpolation model constructed in this invention can maintain stable spatial structure characteristics under different weather conditions and has strong generalization ability.

[0060] S4. All monitoring data from all observation stations at the current time are detected through their corresponding anomaly detection models to determine whether the monitoring data is abnormal. Taking the rainfall monitoring data of a certain observation station 2993 from July 11 to 19, 2024 as an example, the error between the monitoring data and the predicted value is presented in the form of mean square error (MSE). Figure 2(b) is the time-varying curve of rainfall from the 11th to the 19th (unit: mm), and Figure 2(a) is the time-varying curve of error during the same period. The horizontal axis in Figure 2(a) and Figure 2(b) represents time. It can be seen that the curve shapes in Figure 2(a) and Figure 2(b) are roughly the same, and the abnormal data points under different thresholds are concentrated at the time when the rainfall monitoring data changes significantly. This indicates that the detection results have a good correspondence with the actual anomalies, and the model can accurately detect the existence of anomalies in the rainfall data. As can be seen from Figure 2(a), the number of abnormal data points varies under different abnormal thresholds. The number of abnormal data points under the three thresholds of P90, P95 and P99 are 22, 11 and 3 respectively. Therefore, under the three thresholds of P90, P95 and P99, observation station 2993 obtained 22, 11 and 3 abnormal data points respectively.

[0061] A certain region within the system's monitoring range has 430 observation stations. Anomaly detection was performed on the monitoring data of all stations in the same period. The average number of abnormal data points for all monitoring data under three threshold settings (P90, P95, and P99) is shown in Figure 3(a). The average number of abnormal data points shows an upward trend from P90 to P99, indicating that the practice of detecting abnormal data through the anomaly discrimination model has consistency and stability in large-scale detection. The distribution of the anomaly rate of all monitoring data under the P95 threshold is shown in Figure 3(b). In the figure, the horizontal axis represents the anomaly rate and the vertical axis represents the number of abnormal points. The average anomaly rate of the monitoring data of 430 observation stations is 2.6%. As can be seen from Figure 3(b), the anomaly rate is highly concentrated in the low value range, and most data are still judged as normal data. The above situation indicates that the model has a low false alarm rate and can effectively avoid misjudging normal data while ensuring the identification of abnormal data.

[0062] S5. Set a sample time period, and perform cluster analysis on all abnormal data occurring within the sample time period using a clustering algorithm to identify one or more abnormal events; in this embodiment, the clustering algorithm is based on the DBSCAN algorithm, and step S5 specifically includes the following steps: S5.1. Establish sample points belonging to each moment within the sample time period. The number of sample points belonging to a certain moment is equal to the number of outliers that appear at that moment and they correspond one-to-one. The location of any sample point is the same as the location of its corresponding outlier. The location of any outlier is the location of the observation station that obtained the outlier. The location of the observation station is described by two-dimensional coordinates composed of longitude and latitude. All sample points are grouped into the same set for cluster analysis.

[0063] S5.2. Set spatial and temporal thresholds. When a sample point is used as the center point, if the straight-line distance between another sample point and the center point does not exceed the set spatial threshold, and the time difference between the other sample point and the center point does not exceed the set temporal threshold, then the sample point is considered to be within the neighborhood radius of the center point; otherwise, it is considered not to be within the neighborhood radius. Set a minimum number of points, and perform cluster analysis on all sample points using the DBSCAN algorithm to obtain one or more clusters. Determine the corresponding abnormal events based on the morphological characteristics and temporal and spatial scales of the clusters.

[0064] Taking abnormal data from a certain region within the system's monitoring range during a certain period as an example, with a constant time threshold of 30 h, minimum point counts (minpts) of 5, 10, 15, 20, 25, and 30, and spatial thresholds (eps) of 100 km, 120 km, 150 km, and 200 km, the number of clusters is shown in Figure 4(a). In the figure, the horizontal axis represents the minimum point count (minpts), and the vertical axis represents the number of clusters. The noise rate, i.e., the proportion of noise points among all sample points, is shown in Figure 4(b). As can be seen from Figure 4(a), the number of clusters decreases as minpts increases, while the noise rate continuously increases. Considering both clustering effect and meteorological significance, the optimal classification effect was achieved when the number of clusters was between 50 and 80 and the noise rate was between 10 and 20%. By comparing the optimal region with the shape and distribution of each curve in Figures 4(a) and 4(b), a parameter combination of eps = 120 km, minpts = 15, and a time threshold of 30 h was ultimately selected. Based on the clustering algorithm set by these parameters, monitoring data from January to July 2024 within the region extending east to 99.45 °E, west to 90.06 °E, south to 31.12 °N, and north to 36.19 °N were classified. Twenty-five large-scale meteorologically significant anomalous events were successfully identified. The largest anomalous event involved 2296 anomalous data points, lasted 423 h, and had a maximum spatial span of 1026.4 km; the smallest anomalous event involved 21 anomalous data points, lasted 72 h, and had a spatial range of 150.8 km. km, the above results show that the method of the present invention has successfully achieved full-scale identification from regional short-term anomalies to ultra-large persistent anomalies.

[0065] S6. Obtain the rainfall distribution within the spatiotemporal range involved in the abnormal event using the rainfall data interpolation model, and replace or fill in the abnormal data based on the obtained rainfall distribution, specifically including the following steps: S6.1. The sample points to be assigned to the cluster are recorded as points to be replaced. The grid cell where each point to be replaced is located is determined according to the location of the points to be replaced and the grid distribution. The grid nodes contained in the grid cell where a point to be replaced is located are recorded as the associated nodes of that point to be replaced.

[0066] S6.2. Determine the associated nodes at all times within the sample time period. Based on the ERA5 reanalysis data, establish the feature matrix of each associated node by combining the location and time of occurrence of the associated nodes. Input all feature matrices into the rainfall data interpolation model to obtain the rainfall prediction value for each associated node.

[0067] S6.3. Using bilinear interpolation, combined with the location of the point to be swapped and its associated nodes, as well as the predicted rainfall values ​​of the associated nodes, new rainfall data for the point to be swapped is obtained. The new data is then used to replace the abnormal data corresponding to the point to be swapped.

[0068] Taking the 2015 rainfall monitoring data of a certain area within the system's monitoring range as an example, the area contains a total of 249,048 grid nodes. Figure 7(a) shows a comparison between the predicted values ​​at the grid nodes and the actual values ​​derived from the monitoring data. Each data point in the figure corresponds to one grid node. The horizontal axis represents the actual value, and the vertical axis represents the predicted value, both in mm. The red solid line is the 1:1 diagonal. The data shown in the figure has an NSE of 0.841 and a CC of 0.917. The error distribution is shown in Figure 7(b), where the error distribution is concentrated around zero, with a bias of only 0.002 mm. Furthermore, Figure 7(a) shows that the data points are concentrated near the red solid line, indicating almost no systematic bias. These results demonstrate that bilinear interpolation inversion can effectively maintain the discriminative ability of rainfall events, verifying the reliability and accuracy of the data interpolation process from grid node prediction to observation station data in this invention.

Claims

1. A method for cleaning large-scale, long-duration abnormal flood reporting data, the method operating based on a hydrological monitoring system, said hydrological monitoring system comprising multiple hydrological observation stations, each station corresponding to a unique number, and the locations of any two observation stations being different, characterized in that, Specifically, the following steps are included: S1. Set a sampling period, collect meteorological monitoring data measured by all observation stations of the system during the sampling period, determine the location of each observation station, and obtain meteorological-related data of the monitored area of ​​the system during the sampling period through another means; S2. Establish multiple anomaly detection models, which are neural network models. The number of anomaly detection models is equal to the number of observation stations and they correspond one-to-one. All anomaly detection models are trained using meteorological monitoring data from the corresponding observation stations, so that all anomaly detection models have the ability to identify abnormal rainfall. S3. Establish a rainfall data interpolation model, which is a self-attention model. The rainfall data interpolation model is trained using meteorological monitoring data and meteorological related data obtained through another means, so that it has the ability to predict rainfall distribution. S4. During the monitoring process, all monitoring data measured at each observation station at all times are detected by its corresponding anomaly discrimination model to determine whether the monitoring data is abnormal. S5. Set a sample time period, and use a clustering algorithm to perform cluster analysis on all abnormal data that appear within the sample time period to obtain one or more abnormal events; S6. Obtain the rainfall distribution within the spatiotemporal range involved in the abnormal event through the rainfall data interpolation model, and replace or fill the abnormal data based on the obtained rainfall distribution.

2. The method for cleaning large-scale, long-duration abnormal flood reporting data as described in claim 1, characterized in that, The anomaly detection model in step S2 is an autoencoder model. The process of training the autoencoder model specifically includes the following steps: S2.

1. Arrange all rainfall data in the meteorological monitoring data in chronological order to form an original sequence. Set a window length to establish a sliding window. Divide the original sequence into multiple equal-length subsequences through the sliding window. Then, divide the rainfall process corresponding to the original sequence into a first stage, a second stage, and a third stage. Subsequences corresponding to the first stage are assigned to the training set, and subsequences corresponding to the second and third stages are assigned to the validation set and the test set, respectively. Construct the first, second, and third sample matrices based on the training set, validation set, and test set, respectively. Each subsequence is a row vector of the sample matrix, and all subsequences are arranged in chronological order in the sample matrix. S2.

2. Set initial weights and initial weight learning rates, establish an initial AutoEncoder model, input the first sample matrix into the initial AutoEncoder model, and perform multiple rounds of training in combination with other meteorological monitoring data from the first stage. After each round of training, the learning rates of the weights of all hidden layers in the AutoEncoder model are updated by the Adam optimizer. The mean squared error function is used as the first loss function. After each round of training, the obtained model is used to predict the rainfall at all times in the first and second stages based on meteorological monitoring data other than rainfall. Then, the first loss function is used to derive the training loss based on the predicted values ​​of the first stage and the corresponding values ​​in the training set, and the validation loss is derived based on the predicted values ​​of the second stage and the corresponding values ​​in the validation set. S2.

3. Set termination conditions, the termination conditions including a threshold number of attempts. M With the loss threshold, when it appears continuously during training M In each round, the model training is considered to have met the termination condition when the difference between the validation losses of any two adjacent rounds is less than the loss threshold, and the training stops; find the validation loss with the smallest value, and use the model with the smallest validation loss as the sample model; S2.

4. Establish a comprehensive quality scoring model for evaluating the performance of the sample model, and set a scoring threshold. The comprehensive quality scoring model is expressed as the following formula: (1) in, i Number the observation station, Quality i For observation station i Scoring of the sample model L i train For training loss, L i val To verify the loss, N i The number of rounds in the training process. σ i MSE The standard deviation of the error is calculated based on the second-stage predicted values ​​and the corresponding values ​​in the validation set. L To maximize the loss, max N The maximum number of rounds, max σ For the maximum standard deviation, W 1. W 2. W 3. W 4 are all weighting coefficients; The quality of the sample model is calculated using a comprehensive quality scoring model. i If Quality i If the score is not lower than the set scoring threshold, the sample model is considered to have met the performance standard and will be used as the final anomaly detection model. i If the score is below the threshold, the sample model is considered to have met the performance standard, and the model is retrained according to the procedures described in steps S2.2 and S2.

3. 3.S2.

5. The obtained anomaly discrimination model is used to predict the rainfall at all times in the third stage, and the reconstruction error at each time is obtained based on the predicted value of the third stage and the corresponding value in the test set. Then, the anomaly threshold is determined according to the overall distribution of the reconstruction error. When the detection is performed, the anomaly discrimination model obtains the predicted value of the rainfall at a certain time through the meteorological monitoring data at a certain time, and compares the rainfall monitoring data and the predicted value. If the error between the two is greater than the anomaly threshold, the monitoring data is considered to be abnormal data.

4. The method for cleaning large-scale, long-duration abnormal flood reporting data as described in claim 1, characterized in that, In step S3, the rainfall data interpolation model is the Transformer model, and the meteorological data obtained through another method is ERA5 reanalysis data. Step S3 specifically includes the following steps: S3.

1. Establish a grid covering the entire monitoring range of the hydrological monitoring system, with each grid node corresponding to a unique serial number, and determine the grid unit where each observation station is located; based on the rainfall data measured by each observation station at a certain moment, combined with the grid unit where the observation station is located and the location of the grid node, the rainfall at all grid nodes at that moment is obtained by inverse distance weighted interpolation (IDW), and the rainfall at each grid node at all moments during the sampling period is obtained by following the above process; S3.

2. Divide the sampling period into a training phase, a verification phase, and a testing phase. Sort the rainfall at all times at the same grid node by time to obtain a sample sequence. Record the part of the sample sequence corresponding to the training phase as the first sequence, the part corresponding to the verification phase as the second sequence, and the part corresponding to the testing phase as the third sequence. Each grid node obtains its corresponding first, second, and third sequences according to the above process. S3.

3. The ERA5 reanalysis data includes four parameters: convective effective potential energy, temperature, total column water vapor, and total precipitation. The ERA5 reanalysis data is preprocessed, and then bilinear interpolation is used to obtain the ERA5 reanalysis data for each grid node; a time-domain window is set […]. t - τ , t + τ ], t For the current moment, τ The time-domain radius and τ ∈Z + A feature matrix is ​​established for each grid node at all times within the sampling period. The feature matrix contains multiple feature vectors and a static vector. The number of feature vectors is equal to and corresponds one-to-one with the number of times contained in the time-domain window. The feature vector contains four elements, which correspond to the four parameters contained in the ERA5 reanalysis data. The static vector contains the corresponding grid node position parameter. The time corresponding to a certain feature matrix is ​​taken as the current time of the time-domain window to determine the time corresponding to the feature vector of the matrix. The value of the feature vector element is determined by combining the time and the corresponding grid node position. S3.

4. Set the initial learning rate, weight decay coefficient, and learning rate update method, and construct the AdamW optimizer; establish the Transformer model and the second loss function, and train the Transformer model in multiple rounds by combining the first sequence of all grid nodes and the feature matrices of all time points in the training phase. After each round of training, the rainfall at all time points in the validation phase is predicted using the feature matrix of the obtained model set. Then, the validation loss is obtained based on the predicted values ​​in the validation phase and the corresponding time values ​​in the second sequence using the second loss function. After each round of training, the values ​​of all parameters to be learned in the Transformer model are updated by the AdamW optimizer based on the validation loss; set the round threshold. N When training continues N Training stops when the validation loss in the current round does not decrease compared to the previous round. Then, the validation loss with the smallest value is found, and the model with the smallest validation loss is used as the imputation model for the obtained rainfall data. S3.

5. Based on the feature matrix of the test phase, the obtained rainfall data interpolation model is used to predict the rainfall at all times during the test phase. The prediction ability of the obtained rainfall data interpolation model under various weather conditions is evaluated based on the predicted values ​​of the test phase and the corresponding values ​​in the third sequence.

5. A method for cleaning large-scale, long-duration abnormal flood reporting data as described in claim 3, characterized in that, The second loss function established in step S3.4 includes two terms: mean squared error (MSE) and mean absolute error (MAE). The second loss function is expressed as follows: (2) in, k For the number of training rounds, j For grid node number, L k for k Validation loss after each training round T To verify the duration of the phase, w j For grid nodes j The sample weights at each location, under rainy and dry weather conditions. w j The values ​​are different. To pass k The model obtained from 1 round of training t Time grid nodes j Rainfall forecast value, For grid nodes j In the second sequence t Rainfall at any given time λ MSE The mean square error weights are... λ MAE This is the weight for the average absolute error.

6. The method for cleaning large-scale, long-duration abnormal flood reporting data as described in claim 1, characterized in that, The clustering algorithm in step S5 is based on the DBSCAN algorithm, and step S5 specifically includes the following steps: S5.

1. Establish sample points belonging to each moment within the sample time period. The number of sample points belonging to a certain moment is equal to the number of abnormal data appearing at that moment and they correspond one-to-one. The location of any sample point is the same as the location of its corresponding abnormal data. The location of any abnormal data is the location of the observation station from which the abnormal data was obtained. All sample points are grouped into the same set for cluster analysis. S5.

2. Set spatial and temporal thresholds. When a sample point is used as the center point, if the straight-line distance between another sample point and the center point does not exceed the set spatial threshold, and the time difference between the other sample point and the center point does not exceed the set temporal threshold, then the sample point is considered to be within the neighborhood radius of the center point; otherwise, it is considered not to be within the neighborhood radius. Set a minimum number of points, and perform cluster analysis on all sample points using the DBSCAN algorithm to obtain one or more clusters. Determine the corresponding abnormal events based on the morphological characteristics and temporal and spatial scales of the clusters.

7. The method for cleaning large-scale, long-duration abnormal flood reporting data as described in claim 3, characterized in that, Step S6 specifically includes the following steps: S6.

1. The sample points to be assigned to the cluster are recorded as points to be replaced. The grid cell where each point to be replaced is located is determined according to the location of the points to be replaced and the grid distribution. The grid nodes contained in the grid cell where a point to be replaced is recorded as the associated nodes of that point to be replaced. S6.

2. Determine the associated nodes at all times within the sample time period. Based on the ERA5 reanalysis data, establish the feature matrix of each associated node by combining the location and time of occurrence of the associated nodes. Input all feature matrices into the rainfall data interpolation model to obtain the rainfall prediction value for each associated node. S6.

3. Using bilinear interpolation, combined with the location of the point to be swapped and its associated nodes, as well as the predicted rainfall values ​​of the associated nodes, new rainfall data for the point to be swapped is obtained. The new data is then used to replace the abnormal data corresponding to the point to be swapped.

Images