AI-based meteorological station data quality control method integrated with numerical forecast

Abstract

This paper provides an AI-powered meteorological station data quality control method integrated with numerical forecasts. This method involves standardizing multivariate meteorological observation and numerical forecast data, constructing an LSTM-AE physics-driven model to extract observed and forecasted physical features, reconstructing the observed time series using an LSTM decoder, training the model using the deviation between the reconstructed and original sequences as an indicator, assessing abnormalities in the data to be detected, determining quality ratings, and providing feedback to optimize the model. This method achieves high-precision quality control of meteorological station data, improves anomaly detection performance, and ensures the reliability of meteorological data.

Description

Technical Field

[0001] The present invention relates to data control technology, and in particular to an AI meteorological station data quality control method integrated with numerical forecasting. Background Art

[0002] Meteorological observation data is a crucial foundation for weather forecasting, climate research, and disaster prevention and mitigation. With the expansion of meteorological observation networks, the volume of observation data has increased dramatically, making data quality control a critical component of meteorological operations. Traditional meteorological data quality control methods primarily rely on statistical and physical rules for anomaly detection, such as extreme value testing, internal consistency testing, and temporal continuity testing. However, with the increasing complexity of observation systems and the increase in data dimensionality, traditional methods struggle to fully exploit the complex physical connections and spatiotemporal evolutionary characteristics of multivariate meteorological elements.

[0003] In recent years, with the development of artificial intelligence (AI) technology, deep learning methods have been increasingly applied to meteorological data quality control. Recurrent neural networks, in particular, effectively capture the dynamic characteristics of time series data, providing a new technical approach for meteorological data quality control. Long short-term memory (LSTM) networks, with their superior time series modeling capabilities, have become a crucial tool for processing meteorological time series data. Furthermore, autoencoder (AE) structures possess the ability to extract key features from high-dimensional data, enabling dimensionality reduction and reconstruction while preserving the data's essential information.

[0004] Existing quality control methods often separate observation data from numerical forecast data, failing to fully utilize the physical background knowledge contained in the numerical forecast results to guide the quality control process of observation data, resulting in inefficient identification of abnormal data.

[0005] Traditional methods mostly use fixed thresholds and preset rules for anomaly detection. They lack the ability to learn the complex nonlinear relationships between meteorological elements, are difficult to adapt to the differences in meteorological characteristics in different regions and seasons, and are prone to misjudgments and missed judgments.

[0006] Existing machine learning-based methods generally have the "black box" problem. The model decision-making process is difficult to explain, and there is a lack of a continuous optimization mechanism. Model parameters and detection standards cannot be automatically adjusted as new data accumulates, affecting the reliability and practicality of quality control results. Summary of the Invention

[0007] The embodiments of the present invention provide an AI meteorological station data quality control method integrated with numerical forecasts, which can solve the problems in the prior art.

[0008] A first aspect of an embodiment of the present invention provides an AI meteorological station data quality control method integrated with numerical forecasting, comprising:

[0009] Acquiring multivariate meteorological observation data and numerical forecast data from a meteorological station, and performing standardization processing on the multivariate meteorological observation data and the numerical forecast data to form an observation time series dataset and a numerical forecast time series dataset;

[0010] Constructing a physical driving model based on the LSTM-AE structure, extracting the physical relationships and temporal evolution characteristics between meteorological elements in the observation time series dataset through a long short-term memory network encoder based on the physical driving model, parsing the observation time series dataset and the numerical forecast time series dataset to obtain observed physical characteristics and predicted physical characteristics;

[0011] Processing the observed physical features and the forecast physical features based on the LSTM decoder in the physical drive model, completing feature reconstruction according to the physical laws and spatiotemporal evolution characteristics of meteorological elements, and generating a reconstructed observation time series;

[0012] The observation time series dataset is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation indicator to complete the training and parameter optimization of the LSTM codec in the physical drive model;

[0013] Based on the physical driving model, an abnormal state evaluation is performed on the data to be detected, and a quality grade label is determined according to the result of the abnormal state evaluation. The labeling result is fed back to optimize the feature extraction capability and anomaly detection performance of the LSTM-AE structure, thereby realizing quality control and continuous optimization of meteorological station data.

[0014] A physical driving model based on the LSTM-AE structure is constructed. Based on the physical driving model, the physical relationship and temporal evolution characteristics between meteorological elements in the observation time series dataset are extracted through a long short-term memory network encoder. The observation time series dataset and the numerical forecast time series dataset are parsed to obtain the observed physical characteristics and the predicted physical characteristics, including:

[0015] Calculating the observation data of the observation time series data set through a long short-term memory network encoder to obtain an observation mean vector and an observation variance vector, constructing an observation feature normal distribution based on the observation mean vector and the observation variance vector, and sampling from the observation feature normal distribution to obtain a latent variable of the observation data;

[0016] Introducing a priori normal distribution constraints on the network parameters of the physical driving model, and calculating the physical association weights of meteorological elements at different time scales through a Bayesian neural network based on the prior normal distribution constraints;

[0017] Performing weighted reconstruction on the latent variables of the observation data according to the physical association weights to obtain a reconstructed observation feature vector, constructing a conditional probability distribution based on the reconstructed observation feature vector, and performing multiple sampling from the conditional probability distribution to obtain an observed physical feature;

[0018] The numerical forecast time series data set is input into the long short-term memory network encoder to obtain a forecast mean vector and a forecast variance vector of the numerical forecast time series data set, a forecast feature normal distribution is constructed based on the forecast mean vector and the forecast variance vector, latent variables of the forecast data are sampled from the forecast feature normal distribution, and the latent variables of the forecast data are weightedly reconstructed based on the physical association weights to obtain forecast physical characteristics that characterize the physical relationship of meteorological elements.

[0019] The LSTM decoder in the physical drive model processes the observed physical features and the predicted physical features, completes feature reconstruction according to the physical laws and spatiotemporal evolution characteristics of meteorological elements, and generates a reconstructed observation time series including:

[0020] Performing initial reconstruction of the observed physical features and the predicted physical features based on the LSTM decoder to generate initial reconstruction features, calculating a reconstruction error between the initial reconstruction features and the observed time series, and combining the reconstruction error with a hidden state to obtain an error propagation feature;

[0021] Combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term;

[0022] combining the dynamic error compensation term with the initial reconstructed feature to obtain a corrected reconstructed feature, calculating the meteorological physical correlation between air pressure, temperature, humidity, wind speed, and precipitation in the corrected reconstructed feature, and comparing the meteorological physical correlation with the physical correlation between corresponding meteorological elements in the observation time series to obtain a deviation value between the corrected reconstructed feature and the observation time series in terms of physical laws;

[0023] The forgetting weight, input weight and output weight in the LSTM decoder are adjusted according to the deviation value, the cell state and hidden state are updated by the adjusted weights, and a new round of reconstruction of the observed physical characteristics and the predicted physical characteristics is performed based on the updated cell state and the hidden state to generate a reconstructed observation time series that conforms to the physical laws and spatiotemporal evolution characteristics of meteorological elements.

[0024] Combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term includes:

[0025] Mapping the error propagation feature and the initial reconstruction feature to a feature space of the same dimension to obtain a projected error propagation feature and a projected initial reconstruction feature; dividing the projected initial reconstruction feature into a plurality of local regions, calculating a mean and a standard deviation for each of the local regions, constructing a local feature statistic based on the mean and the standard deviation, and combining the local feature statistic with the projected error propagation feature to generate a local compensation coefficient;

[0026] Calculating a correlation matrix between the projected initial reconstructed features and the projected error propagation features, calculating a global feature correlation based on the correlation matrix, and generating a global compensation coefficient according to the global feature correlation;

[0027] Performing feature splicing on the projected initial reconstructed features and the projected error propagation features, generating a weight coefficient according to the result of the feature splicing, and performing a weighted combination of the local compensation coefficient and the global compensation coefficient based on the weight coefficient to obtain a final compensation coefficient;

[0028] The final compensation coefficient is subjected to element-wise multiplication operation with the projected error propagation feature, including: taking each dimension feature of the final compensation coefficient as a weight, weighting the feature of the corresponding dimension of the projected error propagation feature, combining the weighted features of each dimension to obtain a compensation feature, fusing the compensation feature with the initial reconstruction feature to generate a corrected reconstruction feature, and obtaining a dynamic error compensation term according to the difference between the corrected reconstruction feature and the initial reconstruction feature.

[0029] The observation time series data set is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation index, including:

[0030] Performing feature reconstruction on the observation time series data set to obtain a reconstructed observation time series, wherein the reconstructed observation time series has the same temporal distribution characteristics as the original observation time series;

[0031] The timing deviation between the reconstructed observation time series and the original observation time series at corresponding moments is calculated, and the timing deviation is used as an evaluation index to evaluate the reconstruction effect.

[0032] Performing an abnormality evaluation on the data to be detected based on the physical drive model, and determining a quality grade mark according to a result of the abnormality evaluation includes:

[0033] Input the data to be detected into the physical drive model to generate a state feature vector, construct a state transition diagram based on the state feature vector, and generate an initial abnormal state prediction value by analyzing the matching degree between the state transition diagram and the physical constraint conditions;

[0034] A prior distribution is constructed based on the initial abnormal state prediction value, and the Markov chain Monte Carlo method is used to perform stratified repeated sampling on the prior distribution. The node weights and edge probabilities of the state transition diagram are iteratively updated according to the sampling results, and the parameter posterior distribution of the optimal parameter configuration of the state transition diagram is obtained through multiple rounds of iterative optimization;

[0035] Performing multiple stratified sampling on the parameter posterior distribution using the Markov Chain Monte Carlo method, combining the results of the multiple stratified sampling with the optimal parameter configuration to construct a new state transition diagram, and obtaining an abnormal state score by comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram;

[0036] Performing a probability density estimation on the node distribution in the new state transition graph, analyzing the fluctuation characteristics of the node weights and edge probabilities based on the Markov chain Monte Carlo method, and calculating an uncertainty measure based on the fluctuation characteristics combined with the confidence interval of the physical constraints;

[0037] A segmented mapping relationship is established according to the correlation between the numerical distribution of the abnormal state score and the topological structure of the new state transition diagram, and a threshold of the segmented mapping relationship is adaptively adjusted based on the uncertainty metric to determine a quality level mark.

[0038] Combining the results of the multiple stratified samplings with the optimal parameter configuration to construct a new state transition diagram, and comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram to obtain an abnormal state score includes:

[0039] Constructing a node distribution probability matrix based on the results of the multiple stratified samplings, performing a convolution operation on the node distribution probability matrix and the optimal parameter configuration to obtain node weights, constructing a node distribution of a state transition diagram based on the node distribution probability matrix, performing Fourier transform on the node distribution to obtain frequency domain features, optimizing the node distribution probability matrix based on the frequency domain features to obtain edge probabilities, and combining the node weights and the edge probabilities to construct a new state transition diagram;

[0040] The node weights of the new state transition diagram and the original state transition diagram are subjected to the Fourier transform to obtain the node weight frequency domain features, the log likelihood ratio of the node weight frequency domain features is calculated to obtain the node weight difference value, the edge probability of the new state transition diagram and the original state transition diagram is subjected to the Fourier transform to obtain the edge probability frequency domain features, the log likelihood ratio of the edge probability frequency domain features is calculated to obtain the edge probability difference value, and the node weight difference value and the edge probability difference value are weightedly fused to obtain the abnormal state score.

[0041] According to a second aspect of an embodiment of the present invention, an electronic device is provided, including:

[0042] processor;

[0043] a memory for storing processor-executable instructions;

[0044] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0045] According to a third aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0046] The beneficial effects of this application are as follows:

[0047] The AI ​​meteorological station data quality control method integrated with numerical forecasts provided by the present invention can effectively extract the physical relationships and temporal evolution characteristics between meteorological elements by constructing a physical driving model based on the LSTM-AE structure, so that the model has stronger anomaly detection capabilities and improves the accuracy of meteorological station data quality control.

[0048] This method incorporates numerical forecast data into the model as external physical driving information, enhancing the model's ability to understand meteorological physical processes, enabling more accurate distinction between normal meteorological changes and abnormal data, reducing the misjudgment rate, and improving the reliability and scientific nature of quality control.

[0049] The present invention implements a feedback optimization mechanism for quality control results. By feeding back the quality marking results into the model training process, the feature extraction capability and anomaly detection performance of the LSTM-AE structure are continuously optimized, forming a self-improving closed-loop system. This makes the quality control of meteorological station data more intelligent and adaptive, and adapts to the meteorological characteristics of different regions and seasons. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 This is a flow chart of a method for quality control of AI meteorological station data integrated with numerical forecasts according to an embodiment of the present invention;

[0051] Figure 2 This is a histogram comparing the accuracy of meteorological data reconstruction based on the LSTM decoder in an embodiment of the present invention;

[0052] Figure 3 A bar chart comparing the accuracy of different error compensation methods according to an embodiment of the present invention;

[0053] Figure 4 This is a flow chart of abnormal state assessment based on a physical drive model according to an embodiment of the present invention;

[0054] Figure 5 This is a comparison chart of abnormal state scoring accuracy in an embodiment of the present invention. DETAILED DESCRIPTION

[0055] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0056] The technical solution of the present invention is described in detail below with reference to specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0057] Figure 1 FIG. 1 is a flow chart of a method for controlling the quality of AI meteorological station data integrated with numerical forecasts according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0058] Acquiring multivariate meteorological observation data and numerical forecast data from a meteorological station, and performing standardization processing on the multivariate meteorological observation data and the numerical forecast data to form an observation time series dataset and a numerical forecast time series dataset;

[0059] Constructing a physical driving model based on the LSTM-AE structure, extracting the physical relationships and temporal evolution characteristics between meteorological elements in the observation time series dataset through a long short-term memory network encoder based on the physical driving model, parsing the observation time series dataset and the numerical forecast time series dataset to obtain observed physical characteristics and predicted physical characteristics;

[0060] Processing the observed physical features and the forecast physical features based on the LSTM decoder in the physical drive model, completing feature reconstruction according to the physical laws and spatiotemporal evolution characteristics of meteorological elements, and generating a reconstructed observation time series;

[0061] The observation time series dataset is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation indicator to complete the training and parameter optimization of the LSTM codec in the physical drive model;

[0062] Based on the physical driving model, an abnormal state evaluation is performed on the data to be detected, and a quality grade label is determined according to the result of the abnormal state evaluation. The labeling result is fed back to optimize the feature extraction capability and anomaly detection performance of the LSTM-AE structure, thereby realizing quality control and continuous optimization of meteorological station data.

[0063] The input data is preprocessed, and historical observation data from meteorological observation stations is organized into a time series dataset X with dimensions n × m, where n represents the sequence length and m represents the number of features. Features include meteorological elements such as temperature, humidity, and wind speed. The specific values ​​are normalized to a range between 0 and 1. Simultaneously, a numerical forecast dataset X' is compiled for the corresponding time and station, maintaining the same data structure as the observation data.

[0064] A dual-path encoder architecture was designed. The first encoder used a two-layer LSTM network with 128 and 64 hidden layer nodes, respectively. It took observation data X as input and output a 64-dimensional latent space representation Z1. The second encoder used the same network architecture, took numerical forecast data X' as input, and output a 64-dimensional latent space representation Z2. Z1 and Z2 were concatenated along the feature dimension to obtain a 128-dimensional fused feature representation Z.

[0065] The decoder uses a two-layer LSTM network structure with 128 and 256 hidden layer nodes, respectively. It takes as input the fused feature representation Z and outputs reconstructed data with the same dimensions as the original observations. During training, mean squared error is used as the loss function, and the Adam optimizer is used for parameter optimization. The learning rate is set to 0.001, and the number of training epochs is 100.

[0066] In the specific implementation, meteorological station data from January to December 2020 in a certain city was used as an example, involving a total of 20 observation stations. The observation data included three elements: temperature, relative humidity, and wind speed, with a sampling interval of one hour. 80% of the data was used for training and 20% for testing. Numerical forecast data came from the Central Meteorological Observatory's forecast products, with a forecast validity period of 24 hours. After model training and testing, the correlation coefficient between the reconstructed data and the original observation data reached above 0.92, and the root mean square error was kept within 0.15 of the standardized data.

[0067] This method effectively fuses observational data and numerical forecasts through a dual-path encoder structure, allowing the decoder to accurately reconstruct the original data features. This method avoids the subjectivity of manually setting weights in traditional fusion methods, improving the objectivity and accuracy of the fusion results.

[0068] In an optional embodiment, a physical driving model based on an LSTM-AE structure is constructed. Based on the physical driving model, the physical relationship and temporal evolution characteristics between meteorological elements in the observation time series dataset are extracted using a long short-term memory network encoder. The observation time series dataset and the numerical forecast time series dataset are parsed to obtain the observed physical characteristics and the predicted physical characteristics, including:

[0069] Calculating the observation data of the observation time series data set through a long short-term memory network encoder to obtain an observation mean vector and an observation variance vector, constructing an observation feature normal distribution based on the observation mean vector and the observation variance vector, and sampling from the observation feature normal distribution to obtain a latent variable of the observation data;

[0070] Introducing a priori normal distribution constraints on the network parameters of the physical driving model, and calculating the physical association weights of meteorological elements at different time scales through a Bayesian neural network based on the prior normal distribution constraints;

[0071] Performing weighted reconstruction on the latent variables of the observation data according to the physical association weights to obtain a reconstructed observation feature vector, constructing a conditional probability distribution based on the reconstructed observation feature vector, and performing multiple sampling from the conditional probability distribution to obtain an observed physical feature;

[0072] The numerical forecast time series data set is input into the long short-term memory network encoder to obtain a forecast mean vector and a forecast variance vector of the numerical forecast time series data set, a forecast feature normal distribution is constructed based on the forecast mean vector and the forecast variance vector, latent variables of the forecast data are sampled from the forecast feature normal distribution, and the latent variables of the forecast data are weightedly reconstructed based on the physical association weights to obtain forecast physical characteristics that characterize the physical relationship of meteorological elements.

[0073] The physics-driven model of the LSTM-AE architecture consists of two parts: an encoder and a decoder. The encoder compresses high-dimensional input data into a low-dimensional latent representation, while the decoder restores this latent representation to the original data dimensions. The model inputs are observation time series datasets and numerical forecast time series datasets, both of which contain meteorological elements such as air pressure, temperature, humidity, wind speed, and precipitation. The data shape is [N, T, F], where N is the number of samples, T is the length of the time series, and F is the number of meteorological elements. In practice, N = 1000, T = 168 (representing a week of hourly observation data), and F = 5 (representing five meteorological elements).

[0074] The network architecture of the Long Short-Term Memory (LSTM) network encoder consists of an input layer, an LSTM layer, and a fully connected layer. The input layer receives observed time series data; the LSTM layer consists of three stacked layers of LSTM units, each containing 128 hidden neurons, to capture long-term and short-term dependencies in the time series; the fully connected layer maps the LSTM output into two vectors: a mean vector and a variance vector, each with a dimension of 64. To enhance physical constraints, an attention mechanism is incorporated into the encoder, allowing the model to focus more on the physical connections between different meteorological elements. The attention weight matrix is ​​of size [5,5], representing the distribution of attention between each of the five meteorological elements.

[0075] The specific steps of the encoder to process the observed time series data are to input the standardized observation data into the LSTM layer. The LSTM layer controls the information flow through the gating mechanism (including input gate, forget gate and output gate) to capture the time series features. Taking the observation data of a meteorological station from July 1 to July 7, 2022 as an example, the data is standardized and then enters the encoder. The time series features are extracted through the LSTM layer, and the observation mean vector μ is finally output. obs and the observation variance vector σ 2 obs , both are 64-dimensional vectors. Observation mean vector μ obs The first five elements of are [0.32, -0.15, 0.48, 0.21, -0.37], and the observed variance vector σ 2 obs The first 5 elements of are [0.02, 0.04, 0.03, 0.05, 0.02].

[0076] Based on the observed mean vector μ obs and the observation variance vector σ 2 obs Construct the normal distribution of observation characteristics N(μ obs ,σ 2 obs ). This distribution represents the probability distribution of the observed data in the latent space, and each dimension follows a mean of μ obs [i], variance is σ2 obs [i] is an independent normal distribution. The latent variable z of the observed data is sampled from the observed feature normal distribution obs , the sampling process uses the reparameterization technique, that is, z obs =μ obs +ε×sqrt(σ 2 obs ), where ε is a random sample value from the standard normal distribution N(0,1). For the above case, if the first 5 elements of ε are [0.5,-0.2,0.1,-0.4,0.3], then z obs The first 5 elements of are [0.32+0.5×sqrt(0.02),-0.15+(-0.2)sqrt(0.04),0.48+0.1sqrt(0.03),0.21+(-0.4)sqrt(0.05),-0.37+0.3sqrt(0.02)]=[0.39,-0.19,0.49,0.12,-0.34].

[0077] A prior normal distribution constraint is introduced for the network parameters of the physical drive model, transforming the model into a Bayesian neural network. The weights in the network are no longer fixed values ​​but instead follow a prior distribution. In implementation, two parameters, mean and variance, are defined for each weight parameter. The mean is initialized to 0, and the variance is initialized to 0.1. The prior distribution is a normal distribution N(0,0.1), which means that the weight parameters fluctuate around 0 with a variance of 0.1. This constrains the range of weight values ​​and prevents overfitting.

[0078] The physical correlation weights of meteorological elements at different time scales are calculated using a Bayesian neural network, considering three time scales: short-term (hourly), medium-term (daily), and long-term (weekly). The calculation method is to input the observation data of different time scales into the Bayesian network respectively to obtain the corresponding physical correlation weight matrix. The short-term physical correlation weight matrix W short The size is [5,5], which means the mutual influence of the five meteorological elements on the hourly scale; the medium-term physical correlation weight matrix W medium The size is [5,5]; the long-term physical association weight matrix W long Also [5,5].

[0079] The short-term physical association weight matrix W short For example, its element W short [i,j] represents the impact intensity of the i-th meteorological factor on the j-th meteorological factor. For example, W short [0,1]=0.7 means that air pressure (index 0) has a stronger influence on temperature (index 1); W short[1,2]=0.6 means that temperature (index 1) has a moderate effect on humidity (index 2).

[0080] The process of weighted reconstruction of latent variables of observation data according to physical association weights is to combine the physical association weight matrices of three time scales with the latent variables z obs Combined to generate the reconstructed observation feature vector. The specific implementation method is to transform the latent variable z obs It is divided into three sub-vectors, representing short-term, medium-term and long-term features respectively. The dimension of each sub-vector is 64 / 3≈21. The three sub-vectors are multiplied by the corresponding physical association weight matrix respectively, and then the results are concatenated to obtain the reconstructed observation feature vector r obs , the dimension is 64. For example, suppose z obs The first 21 dimensions represent short-term features, so these 21-dimensional features are reorganized into a matrix of [5,21 / 5]=[5,4] (each row represents a meteorological element), and W short After multiplication, we get the weighted short-term features.

[0081] Based on the reconstructed observation feature vector r obs Construct the conditional probability distribution p(x|r obs ), which represents the given reconstructed feature r obs The probability distribution of the original observation data x under the condition of . The implementation method is to convert r obs The input is a small neural network, and the output is the parameters of the conditional distribution. The neural network contains two fully connected layers and the hidden layer size is 128. The network output is the mean vector μ of the conditional distribution cond and variance vector σ 2 cond , both are vectors of the same dimension as the original observation data, that is, [T,F]=[168,5].

[0082] From the conditional probability distribution p(x|r obs ) to obtain the observed physical characteristics. The sampling times are set to 50 times, and each sampling obtains a vector with the same dimension as the original observation data. The specific sampling method is to start with a mean of μ cond , the variance is σ 2 cond The results of 50 samples are randomly drawn from the normal distribution. The average value is taken to obtain the final observed physical characteristics o feat , the dimension is [T,F]=[168,5]. For the above case, if the conditional distribution parameter at a certain time point t=24 (representing the last hour of the first day) is μ cond

[24] =[1.0,22.5,65.0,3.2,0.0],σ 2 cond

[24] =[0.01, 0.25, 1.0, 0.04, 0.0001] (corresponding to the standardized air pressure, temperature, humidity, wind speed and precipitation respectively), then the observed physical characteristics at this time point are o feat

[24] is [1.01, 22.6, 64.8, 3.25, 0.0], which represents the physical characteristics of each meteorological element at that time point.

[0083] The processing of inputting the numerical forecast time series dataset into the LSTM encoder is similar to that of the observation data, except that the input data comes from a different source. The numerical forecast data also contains air pressure, temperature, humidity, wind speed, and precipitation, and the data shape is [N, T, F]. After the data is processed by the encoder, the forecast mean vector μ is obtained. pred and the forecast variance vector σ 2 pred , both are 64-dimensional vectors.

[0084] Based on the predicted mean vector μ pred and the forecast variance vector σ 2 pred Construct the normal distribution of prediction characteristics N(μ pred ,σ 2 pred ). Sample the latent variable z of the forecast data from this distribution pred , the sampling method also uses the reparameterization technique. For example, if the forecast mean vector μ pred The first five elements of are [0.35,-0.18,0.50,0.25,-0.40], and the predicted variance vector σ 2 pred The first 5 elements of are [0.03, 0.05, 0.04, 0.06, 0.03], and the first 5 elements of the random sampling value ε are [0.4, -0.3, 0.2, -0.5, 0.2], then the latent variable z of the forecast data is pred The first 5 elements of are [0.35+0.4×sqrt(0.03),-0.18+(-0.3)sqrt(0.05),0.50+0.2sqrt(0.04),0.25+(-0.5)sqrt(0.06),-0.40+0.2sqrt(0.03)]=[0.42,-0.25,0.54,0.13,-0.37].

[0085] Based on the physical association weights, the latent variables of the forecast data are weighted and reconstructed to obtain the forecast physical characteristics that characterize the physical relationship of meteorological elements. The weighted reconstruction process is the same as that of the observation data, except that the input latent variables become z pred . predIt is divided into three sub-vectors, representing short-term, medium-term and long-term characteristics respectively, multiplied by the corresponding physical association weight matrix, and then the results are spliced ​​to obtain the reconstructed prediction feature vector r pred Based on r pred Construct the conditional probability distribution p(x|r pred ), multiple samples are taken from the distribution and averaged to obtain the final predicted physical characteristics p feat , the dimension is [T,F]=[168,5].

[0086] The observed physical characteristics o obtained by the above method feat and predicted physical characteristics p feat This method not only captures the statistical properties of the raw data but also incorporates the physical relationships and temporal evolution characteristics between meteorological elements, providing a reliable foundation for subsequent feature reconstruction and quality control. Experiments show that compared with traditional methods, the physical features extracted by this method significantly improve the accuracy of anomaly detection while maintaining data integrity. It is particularly capable of identifying anomalous data caused by violations of physical laws.

[0087] In an optional embodiment, the observed physical features and the forecast physical features are processed based on the LSTM decoder in the physical drive model, and feature reconstruction is completed according to the physical laws and spatiotemporal evolution characteristics of meteorological elements to generate a reconstructed observation time series, including:

[0088] Performing initial reconstruction of the observed physical features and the predicted physical features based on the LSTM decoder to generate initial reconstruction features, calculating a reconstruction error between the initial reconstruction features and the observed time series, and combining the reconstruction error with a hidden state to obtain an error propagation feature;

[0089] Combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term;

[0090] combining the dynamic error compensation term with the initial reconstructed feature to obtain a corrected reconstructed feature, calculating the meteorological physical correlation between air pressure, temperature, humidity, wind speed, and precipitation in the corrected reconstructed feature, and comparing the meteorological physical correlation with the physical correlation between corresponding meteorological elements in the observation time series to obtain a deviation value between the corrected reconstructed feature and the observation time series in terms of physical laws;

[0091] The forgetting weight, input weight and output weight in the LSTM decoder are adjusted according to the deviation value, the cell state and hidden state are updated by the adjusted weights, and a new round of reconstruction of the observed physical characteristics and the predicted physical characteristics is performed based on the updated cell state and the hidden state to generate a reconstructed observation time series that conforms to the physical laws and spatiotemporal evolution characteristics of meteorological elements.

[0092] The LSTM decoder receives as input the observed and predicted physical features. Both features are 256-dimensional vectors, representing the key physical information extracted from the observed and predicted data, respectively. The decoder structure consists of an input gate, a forget gate, an output gate, and a cell state. Initially, the forget weight Wf is set to 0.3, the input weight Wi is set to 0.4, and the output weight Wo is set to 0.3. The LSTM decoder first concatenates the observed and predicted physical features to form a 512-dimensional input vector. It then processes this input vector through a gating mechanism to generate the initial reconstructed features.

[0093] After receiving the input vector, the decoder transforms it using its internal parameter matrix to obtain an intermediate state vector. The cell state Ct is initially set to an all-zero vector with a dimension of 128; the hidden state ht is also initially set to an all-zero vector with a dimension of 128. The forget gate controls the retention ratio of the previous cell state through Wf; the input gate determines the acceptance of the current input information through Wi; and the output gate controls the amount of information output from the cell state to the hidden state through Wo. The decoding process is iterative, updating the cell state and hidden state at each time step, ultimately outputting the initial reconstructed features with the dimension of the original input time series, that is, including the five meteorological elements: air pressure, temperature, humidity, wind speed, and precipitation.

[0094] For example, the original observation data from a meteorological station consists of hourly observations from May 1st to May 5th, representing 120 meteorological element data points. After processing with an LSTM encoder, the observed physical features are obtained. The numerical forecast data for the same period is also processed with the same encoder to obtain the predicted physical features. These two features are then fed into a decoder to generate the initial reconstructed features. These features have the same shape as the original observation data: [120, 5], representing the values ​​of the five meteorological elements at the 120 time points.

[0095] The mean squared error (MSE) metric is used to calculate the reconstruction error between the initial reconstructed features and the observed time series. For each meteorological element j at each time point t, the square of the difference between the reconstructed value and the observed value is calculated. This is then averaged across all time points and meteorological elements to obtain the overall reconstruction error. In the example above, the mean reconstruction error for the initial reconstruction process is 0.08, indicating an average deviation of 8% between the reconstructed and observed values.

[0096] The reconstruction error and hidden state are combined to form the error propagation feature. The reconstruction error is represented as a matrix E with a shape of [120, 5]. The hidden state ht at the final time step is represented as a vector H with a dimension of 128. A fully connected layer is used to convert E into a vector E' of the same dimension as H. E' is then concatenated with H to form a 256-dimensional error propagation feature EP. This feature captures the error distribution during the reconstruction process and the internal state information of the LSTM decoder, providing a basis for subsequent correction.

[0097] The error propagation feature EP is combined with the initial reconstruction feature IF to calculate the compensation coefficient. EP is mapped to the same dimensional space as IF through a fully connected layer to obtain the mapped EP'; then the correlation coefficients of EP' and each dimension of IF are calculated to form a correlation matrix CM; based on CM, dimension pairs with significant correlation (dimensional pairs with absolute values ​​of correlation coefficients greater than 0.6) are selected, the covariance of these dimension pairs is calculated and normalized to obtain the preliminary compensation coefficient; finally, the preliminary compensation coefficient is weighted and adjusted according to the confidence level of the error propagation feature (evaluated by the stability of the hidden state) to obtain the final compensation coefficient matrix C , the shape is the same as the initial reconstructed feature, which is [120,5].

[0098] In a real-world example, for the temperature data at 2:00 PM on May 3rd, the initial reconstruction value was 28.6°C, while the actual observed value was 27.2°C, resulting in a reconstruction error of 1.4°C. Using the above method, the compensation coefficient for the temperature dimension at that time point was calculated to be -0.42, indicating that a negative adjustment was required to the initial reconstruction value.

[0099] Based on the compensation coefficient, the initial reconstruction feature is corrected to obtain the dynamic error compensation term, and the compensation coefficient matrix C Perform element-wise multiplication with the error propagation feature map EP' to obtain the base compensation value. Then, weight factors are set based on the physical characteristics of different meteorological elements, such as 0.8 for temperature and pressure, 0.7 for humidity, 0.5 for wind speed, and 0.6 for precipitation. Multiplying the base compensation value with the corresponding weight factors yields the final dynamic error compensation term DC, which has a shape of [120, 5]. For the temperature data in the above example, the calculated dynamic error compensation term is -1.1°C.

[0100] The dynamic error compensation term DC is added to the initial reconstructed feature IF to obtain the corrected reconstructed feature CF. For the above example, the corrected temperature value is 28.6°C + (-1.1°C) = 27.5°C, which is closer to the actual observed value of 27.2°C than the initial reconstructed value of 28.6°C. The reconstruction error is reduced from 1.4°C to 0.3°C.

[0101] When calculating the physical correlation between meteorological elements in the corrected reconstructed features CF, for the pressure-temperature relationship, the correlation of the change rates between the two is calculated; for the temperature-humidity relationship, its physical consistency is verified based on the dew point temperature formula; for the pressure-wind speed relationship, whether the relationship between the pressure gradient and wind speed conforms to meteorological laws is verified; for the temperature-precipitation relationship, whether the temperature change before and after precipitation conforms to common meteorological phenomena is verified. Through these physical law tests, the physical correlation matrix PM between each pair of meteorological elements is obtained. CF , the shape is [5,5], which indicates the degree of physical correlation between the five meteorological elements.

[0102] The same method is used to calculate the physical correlation matrix PM between the corresponding meteorological elements in the observation time series OBS , PM CF With PM OBS Compare and calculate the difference between the two to get the deviation matrix based on physical laws D The smaller the deviation, the more accurately the reconstruction matches the actual laws of meteorological physics. In the above example, the deviation between the reconstructed and observed temperature-humidity relationships was 0.05, indicating a very close physical correlation between the two.

[0103] According to the deviation matrix D Adjust the weights in the LSTM decoder. Meteorological elements with larger deviations require larger weight adjustments. This adjustment is performed by calculating the average deviation of each meteorological element in the deviation matrix and converting the deviation into a weight adjustment factor. For example, if the average deviation of the temperature element is 0.07, the corresponding adjustment factor is -0.03, indicating that the weight associated with temperature needs to be slightly reduced. On the other hand, if the average deviation of precipitation is 0.15, the corresponding adjustment factor is -0.08, indicating that the weight associated with precipitation needs to be significantly reduced.

[0104] The original weights are added to the adjustment factor to obtain the adjusted weights. For example, the forgetting weight associated with temperature processing originally had a value of 0.3, but after adjustment, it becomes 0.3 + (-0.03) = 0.27. The input weight associated with precipitation processing originally had a value of 0.4, but after adjustment, it becomes 0.4 + (-0.08) = 0.32. This weight adjustment method based on physical bias can improve the reconstruction accuracy of specific meteorological elements in a targeted manner.

[0105] The process of updating the cell state and hidden state using the adjusted weights is as follows: using the new forgetting weights to control the degree of forgetting of the previous cell state; using the new input weights to control the degree of acceptance of the current input information; and using the new output weights to control the output of the hidden state. The update formula remains unchanged, but the parameter values ​​change, affecting the way information flows and the reconstruction results.

[0106] Based on the updated cell states and hidden states, a new round of reconstruction of the observed and predicted physical features is performed to generate the final reconstructed observation time series. This new round of reconstruction uses the same decoder structure and process as the initial reconstruction, but with adjusted weight parameters. Experimental results show that after the weight adjustment and the new round of reconstruction, the reconstruction accuracy of meteorological elements is significantly improved, with the average reconstruction error reduced from the initial 0.08 to 0.03. Specifically for each element, the reconstruction error of temperature is reduced by 78%, humidity by 72%, air pressure by 85%, wind speed by 65%, and precipitation by 58%.

[0107] The reconstructed observation time series finally generated is not only numerically close to the original observation data, but more importantly, it is highly consistent with the real observation data in terms of the physical correlation between meteorological elements, ensuring that the reconstruction results conform to the physical laws and spatiotemporal evolution characteristics in meteorology, and providing a reliable basis for subsequent meteorological data quality control.

[0108] Figure 2This is a bar chart comparing the accuracy of meteorological data reconstruction based on the LSTM decoder in an embodiment of the present invention. The figure compares the performance of three prediction models (traditional LSTM reconstruction, dynamic error compensation reconstruction, and physical law constraint reconstruction) on five meteorological elements. Among them, the traditional LSTM reconstruction represents the prediction model of the basic long-short-term memory network, and its performance is relatively weak in various indicators, with temperature prediction accuracy of 76.2%, humidity prediction of 72.5%, air pressure prediction of 81.3%, wind speed prediction of 61.8%, and precipitation prediction of only 57.4% at the lowest; the dynamic error compensation reconstruction represents a model that introduces a dynamic error correction mechanism, which has shown significant improvements compared with the traditional model, with temperature prediction reaching 83.5%, humidity prediction reaching 78.3%, air pressure prediction improved to 88.2%, wind speed prediction reached 70.6%, and precipitation prediction increased to 68.1%; the physical law constraint reconstruction represents an optimization model combined with physical law constraints, achieving the best results in all indicators, among which the temperature prediction accuracy reached 91.8%, humidity prediction reached 83.7%, air pressure prediction reached the highest 96.5%, wind speed prediction reached 77.2%, and even the most challenging precipitation prediction reached 74.3%. The overall data shows that the model that introduces physical law constraints has an average improvement of 15-20 percentage points compared to the traditional LSTM model, fully verifying the significant effect of physical law constraints in improving the accuracy of meteorological element forecasts.

[0109] In an optional implementation, combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term includes:

[0110] Mapping the error propagation feature and the initial reconstruction feature to a feature space of the same dimension to obtain a projected error propagation feature and a projected initial reconstruction feature; dividing the projected initial reconstruction feature into a plurality of local regions, calculating a mean and a standard deviation for each of the local regions, constructing a local feature statistic based on the mean and the standard deviation, and combining the local feature statistic with the projected error propagation feature to generate a local compensation coefficient;

[0111] Calculating a correlation matrix between the projected initial reconstructed features and the projected error propagation features, calculating a global feature correlation based on the correlation matrix, and generating a global compensation coefficient according to the global feature correlation;

[0112] Performing feature splicing on the projected initial reconstructed features and the projected error propagation features, generating a weight coefficient according to the result of the feature splicing, and performing a weighted combination of the local compensation coefficient and the global compensation coefficient based on the weight coefficient to obtain a final compensation coefficient;

[0113] The final compensation coefficient is subjected to element-wise multiplication operation with the projected error propagation feature, including: taking each dimension feature of the final compensation coefficient as a weight, weighting the feature of the corresponding dimension of the projected error propagation feature, combining the weighted features of each dimension to obtain a compensation feature, fusing the compensation feature with the initial reconstruction feature to generate a corrected reconstruction feature, and obtaining a dynamic error compensation term according to the difference between the corrected reconstruction feature and the initial reconstruction feature.

[0114] The combined processing of error propagation features and initial reconstruction features first requires mapping both to a feature space of the same dimensionality. Assuming the original dimension of the error propagation features EP is 64 and the original dimension of the initial reconstruction features IR is 128, they are each mapped to a 96-dimensional feature space through a fully connected layer. In specific implementation, the error propagation features EP are linearly transformed using the parameter matrix WEP (size 64×96) to obtain the projected error propagation features EP'; similarly, the initial reconstruction features IR are linearly transformed using the parameter matrix WIR (size 128×96) to obtain the projected initial reconstruction features IR'. During the transformation process, the elements of each parameter matrix are optimized using the gradient descent method, so that the transformed features retain the key information of the original features.

[0115] The initial reconstructed feature IR' after projection needs to be divided into multiple local regions to capture the error characteristics at different locations. In practical applications, the 96-dimensional IR' feature is divided into 8 local regions, each containing 12 continuous feature dimensions. For the i-th local region IRi' (i=1,2...8), the mean μi and standard deviation σi of all eigenvalues ​​in the region are calculated. For example, for the first local region IR1' (containing the 1st to 12th dimension features of IR'), its mean μ1 is the arithmetic mean of these 12 eigenvalues, and the standard deviation σ1 is a measure of the degree of dispersion of these 12 eigenvalues.

[0116] Based on the calculated mean μi and standard deviation σi, a local feature statistic LSi is constructed. For the i-th local region, its statistic LSi is a two-tuple (μi, σi). Taking real data as an example, in a certain calculation, the mean μ1=0.32 and the standard deviation σ1=0.15 of IR1', the corresponding local feature statistic LS1=(0.32, 0.15). The local feature statistic LSi is combined with the corresponding region EPi' in the projected error propagation feature EP'. The combination is performed by comparing the mean μEPi of EPi' with μi in LSi to calculate the relative deviation; and comparing the standard deviation σEPi of EPi' with σi in LSi to calculate the fluctuation consistency. Based on these two comparisons, the local compensation coefficient LCi is generated. Typically, if μEPi is close to μi, but σEPi differs significantly from σi, it indicates that the region requires significant compensation adjustment.

[0117] Calculate the correlation matrix CM between the projected initial reconstructed features IR' and the projected error propagation features EP'. The size of the correlation matrix CM is 96×96, where the element CM[j,k] represents the correlation between the j-th dimension of IR' and the k-th dimension of EP'. The correlation is calculated by calculating the covariance of the two features on the training dataset and normalizing it. The value range is [-1, 1]. Correlations closer to 1 indicate a stronger positive correlation, closer to -1 indicates a stronger negative correlation, and closer to 0 indicates no significant correlation.

[0118] The global feature correlation (GR) is calculated based on the correlation matrix (CM). This GR is obtained by statistically analyzing elements in CM with large absolute values ​​(e.g., |CM[j,k]| > 0.5). Specifically, the proportion of elements meeting this condition relative to the total number of elements is calculated, as well as the average absolute value of these elements. In practice, if the proportion of elements meeting |CM[j,k]| > 0.5 is 30%, and the average absolute value of these elements is 0.72, then a moderate degree of global correlation exists between the two features.

[0119] A global compensation coefficient GC is generated based on the global feature correlation GR. When GR indicates a strong correlation (e.g., the proportion of correlated elements is >40% and the average absolute value is >0.6), a larger global compensation coefficient (e.g., GC = 0.8) is set; when GR indicates a moderate correlation, a medium global compensation coefficient (e.g., GC = 0.5) is set; when GR indicates a weak correlation, a smaller global compensation coefficient (e.g., GC = 0.2) is set. The global compensation coefficient GC is a 96-dimensional vector, with each dimension corresponding to a dimension in the feature space.

[0120] Concatenate the projected initial reconstructed feature IR' with the projected error propagation feature EP'. This concatenation is done by directly concatenating IR' and EP' end-to-end to form a 192-dimensional concatenated feature CF. For example, if the first three dimensions of IR' are [0.2, 0.4, -0.1] and the first three dimensions of EP' are [0.3, 0.1, -0.2], then the first six dimensions of the concatenated CF are [0.2, 0.4, -0.1, 0.3, 0.1, -0.2].

[0121] The weight coefficient WC is generated based on the feature concatenation result CF. The weight coefficient WC is calculated using a simplified neural network consisting of one hidden layer (192--32) and one output layer (32--1). The network input is CF, and the output is a scalar value in the range [0, 1]. This value represents the weight of the local compensation coefficient in the final compensation coefficient calculation, while the weight of the global compensation coefficient is 1 minus this value.

[0122] The local compensation coefficients LC and the global compensation coefficients GC are weighted together based on the weight coefficient WC to obtain the final compensation coefficient FC. The calculation method is: FC = WC × LC + (1-WC) × GC. LC is the 96-dimensional vector formed by connecting the compensation coefficients LCi of the eight local regions, and GC is the 96-dimensional global compensation coefficient calculated previously. In a real-world case, if WC = 0.6, LC = 0.75, and GC = 0.45 in a certain dimension, then FC for that dimension = 0.6 × 0.75 + 0.4 × 0.45 = 0.45 + 0.18 = 0.63.

[0123] Perform element-wise multiplication of the final compensation coefficient FC with the projected error propagation feature EP'. The specific operation is: for the j-th dimension feature EP'[j] of EP', multiply it with the j-th dimension coefficient FC[j] of FC to obtain the weighted eigenvalue WF[j]=EP'[j]×FC[j]. The weighted eigenvalues ​​of all dimensions together constitute the 96-dimensional compensation feature CF. For example, if EP'

[10] =0.42 and FC

[10] =0.63, then WF

[10] =0.42×0.63=0.2646.

[0124] The compensated features CF are fused with the initial reconstructed features IR to generate the corrected reconstructed features CR. Since the dimension of CF is 96 and the dimension of IR is 128, CF must first be mapped back to the feature space of IR. Using a fully connected layer, the parameter matrix WCF (96×128) converts CF into a 128-dimensional CF'. CF' is then added to IR: CR = IR + CF'. Addition is used here rather than multiplication to preserve the main structure of the original reconstructed features while making targeted adjustments to specific dimensions.

[0125] The dynamic error compensation term DC is calculated based on the difference between the corrected reconstructed feature CR and the initial reconstructed feature IR. This calculation is: DC = CR - IR = CF'. The dynamic error compensation term DC is also a 128-dimensional vector, representing the adjustments made to each dimension of the initial reconstructed feature.

[0126] In a specific application case, in a certain image reconstruction task, the 57th dimension eigenvalue of IR was originally 0.835, indicating the intensity of a certain texture feature in the image. After the above compensation calculation, the corresponding dynamic error compensation term DC

[57] =0.127, and the corrected reconstructed feature CR

[57] =0.962. Practical verification shows that the corrected eigenvalue more accurately reflects the texture features of the original image, reduces the error introduced by the feature extraction process, and thus improves the quality of the reconstructed image. The reconstruction error on the test dataset was reduced from an average of 0.085 to 0.042, improving the accuracy by about 50%.

[0127] The introduction of dynamic error compensation effectively mitigates information loss during feature reconstruction, particularly for high-frequency details and complex textures. By comprehensively analyzing local and global features, this method provides adaptive compensation strategies for different types of errors, significantly improving the accuracy and stability of the reconstruction results.

[0128] Figure 3 This is a bar chart comparing the accuracy of different error compensation methods in the embodiment of the present invention. The figure shows the performance comparison of three different prediction methods in meteorological element prediction scenarios. Among them, the traditional reconstruction method represents the basic model using conventional data reconstruction technology, and its performance is relatively conservative, reaching 67.0% in the temperature prediction scenario, 55.0% in the precipitation prediction scenario, and 62.0% in the wind speed prediction scenario; the local compensation method represents an algorithm that introduces a local data correction mechanism, and its performance has been significantly improved, with the temperature prediction accuracy increased to 76.0%, the precipitation prediction accuracy reached 66.0%, and the wind speed prediction accuracy increased to 72.0%; the present technical solution represents an innovative comprehensive optimization algorithm, which shows the best performance in all test scenarios, with the accuracy of the temperature prediction scenario reaching 90.0%, the precipitation prediction scenario reaching 81.0%, and the wind speed prediction scenario reaching 87.0%. From the overall data, the present technical solution has an average improvement of about 20-25 percentage points compared to the traditional reconstruction method, and a performance improvement of 10-15 percentage points compared to the local compensation method, which fully demonstrates the significant effect of the solution in improving the accuracy of meteorological element prediction. Especially in precipitation scenarios that are more difficult to predict, this technical solution still maintains a high prediction accuracy, reflecting its robustness and adaptability in complex meteorological element prediction tasks.

[0129] In an optional embodiment, the observation time series data set is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation index, including:

[0130] Performing feature reconstruction on the observation time series data set to obtain a reconstructed observation time series, wherein the reconstructed observation time series has the same temporal distribution characteristics as the original observation time series;

[0131] The timing deviation between the reconstructed observation time series and the original observation time series at corresponding moments is calculated, and the timing deviation is used as an evaluation index to evaluate the reconstruction effect.

[0132] Divide the observation time series dataset into chronological order. For example, divide the first 800 time points of a 1000-time-point observation time series dataset into a training set, and the last 200 time points into a validation set. This division method maintains the coherence and temporal characteristics of the time series, which is conducive to the model's learning of time series features.

[0133] The autoencoder model can be used to implement feature reconstruction of observation time series datasets. The autoencoder consists of two parts: an encoder and a decoder. The encoder maps the original observation time series to a latent space representation, and the decoder reconstructs the latent space representation into a sequence with the same temporal distribution characteristics as the original observation time series. Specifically, for a certain observation time series sample X, which contains data from n time points, the encoder maps X to a latent representation Z, and the decoder maps Z to a reconstructed sequence X'. In practical applications, such as in a power load forecasting scenario, the original observation time series is the power load data for each hour within a week, that is, 168 data points. After feature reconstruction through the autoencoder, the reconstructed time series should also contain 168 data points and maintain the temporal distribution characteristics of the original sequence, such as periodicity and trend.

[0134] To ensure that the reconstructed observation time series has the same temporal distribution characteristics as the original observation time series, the mean squared error (MSE) loss function can be used when training the autoencoder model, and a temporal correlation constraint can be added to the loss function. This constraint is implemented by calculating the difference in the autocorrelation coefficients between the original and reconstructed sequences, ensuring that the reconstructed sequence maintains the temporal dependencies of the original sequence. For example, in practical applications, temperature observation data from a meteorological monitoring station exhibits significant intraday fluctuations and seasonal variations. The temperature series reconstructed using these features should also reflect these temporal characteristics.

[0135] After the feature reconstruction is completed, it is necessary to calculate the time series deviation between the reconstructed observation time series and the original observation time series at the corresponding time, and use the time series deviation as an evaluation indicator. The calculation method of the time series deviation can be the absolute value or relative error of the difference between the values ​​at the corresponding time. Specifically, assuming that the original observation time series is X=(x1,x2...x n ), the reconstructed observation time series is X'=(x'1,x'2...x' n ), then the absolute deviation at each moment |x can be calculated i -x' i |or relative deviation|x i -x' i | / |x i |.

[0136] To comprehensively evaluate the reconstruction effect, various statistical indicators can be calculated, including mean absolute deviation, root mean square deviation, and maximum deviation. For example, for time series data generated by an industrial sensor, if the original sequence value at a certain time point is 20.5 and the reconstructed sequence value at the corresponding time point is 19.8, the absolute deviation is 0.7, and the relative deviation is 0.034, or 3.4%. By calculating the mean absolute deviation of the entire series, a value such as 0.65 can be obtained, indicating the average deviation of the reconstructed series from the original.

[0137] In practical applications, time series deviation evaluation metrics can be further refined into deviations at different time scales, such as short-term deviation (hourly), medium-term deviation (daily), and long-term deviation (monthly). This is particularly important for time series with multi-scale characteristics. For example, the sales data of a retail company exhibits both daily fluctuations and weekend effects and seasonal variations. By calculating deviations at different time scales, we can comprehensively evaluate the reconstruction model's ability to capture multi-scale time series characteristics.

[0138] Statistical tests, such as the Kolmogorov-Smirnov test, can be used to assess the distributional consistency of the reconstructed sequence and the original sequence. Autocorrelation functions can be used to compare the reconstructed sequence and assess whether it maintains the temporal dependency structure of the original sequence. For example, for a financial market yield time series, the autocorrelation coefficients of the original and reconstructed sequences can be calculated to compare and examine the similarity in their statistical properties.

[0139] Visualization techniques are also an effective tool in the evaluation process. By plotting the original and reconstructed sequences in the same coordinate system, one can visually observe the degree of fit and deviation between the two over different time periods. For example, for an air quality index series from an environmental monitoring station, a time series curve comparison can be used to reveal the performance of the reconstructed model in capturing outliers and identifying sudden changes.

[0140] Based on the time series deviation evaluation index, different feature reconstruction methods can be compared and optimized. For example, in an intelligent traffic flow prediction system, different autoencoder structures (such as variational autoencoders and recurrent neural network autoencoders) can be tried, and the most suitable reconstruction method can be selected based on the time series deviation index to achieve high-precision reconstruction of traffic flow time series, providing a reliable foundation for subsequent anomaly detection and prediction tasks.

[0141] In an optional embodiment, performing abnormality assessment on the data to be detected based on the physical drive model, and determining a quality level mark according to a result of the abnormality assessment includes:

[0142] Input the data to be detected into the physical drive model to generate a state feature vector, construct a state transition diagram based on the state feature vector, and generate an initial abnormal state prediction value by analyzing the matching degree between the state transition diagram and the physical constraint conditions;

[0143] A prior distribution is constructed based on the initial abnormal state prediction value, and the Markov chain Monte Carlo method is used to perform stratified repeated sampling on the prior distribution. The node weights and edge probabilities of the state transition diagram are iteratively updated according to the sampling results, and the parameter posterior distribution of the optimal parameter configuration of the state transition diagram is obtained through multiple rounds of iterative optimization;

[0144] Performing multiple stratified sampling on the parameter posterior distribution using the Markov Chain Monte Carlo method, combining the results of the multiple stratified sampling with the optimal parameter configuration to construct a new state transition diagram, and obtaining an abnormal state score by comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram;

[0145] Performing a probability density estimation on the node distribution in the new state transition graph, analyzing the fluctuation characteristics of the node weights and edge probabilities based on the Markov chain Monte Carlo method, and calculating an uncertainty measure based on the fluctuation characteristics combined with the confidence interval of the physical constraints;

[0146] A segmented mapping relationship is established according to the correlation between the numerical distribution of the abnormal state score and the topological structure of the new state transition diagram, and a threshold of the segmented mapping relationship is adaptively adjusted based on the uncertainty metric to determine a quality level mark.

[0147] like Figure 4 As shown, the method includes:

[0148] After the data to be tested is input into the physical drive model, it is extracted through an LSTM encoder to obtain a state feature vector. This state feature vector is a 64-dimensional vector, with each dimension representing a physical correlation characteristic between meteorological elements. Taking the four meteorological elements of temperature, humidity, air pressure, and wind speed as an example, the original data input to the model has a dimension of [n, 4], where n is the length of the time series. After processing by the encoder, the resulting state feature vector is a

[64] -dimensional vector.

[0149] Based on the acquired state eigenvectors, a state transition graph is constructed, consisting of multiple nodes and connecting edges. Nodes represent meteorological states, and edges represent transition probabilities between states. In implementation, the 16 most significant dimensions of the eigenvectors are selected and clustered into eight representative state nodes using the K-means clustering algorithm. The initial weight of each node is determined by the proportion of the corresponding sample size to the total sample size. Edge weights are calculated from the state transition frequency between adjacent time steps, forming the initial state transition matrix T, with dimensions [8,8].

[0150] The state transition diagram is then analyzed for compatibility with predefined physical constraints. These constraints include a set of 12 rules, including a daily temperature range of no more than 15°C, relative humidity between 0% and 100%, and a pressure change rate of no more than 3 hPa / hour. By calculating the degree of conformance between the state of each node in the state transition diagram and the physical constraints, an initial abnormal state prediction value is derived. This value is a scalar between 0 and 1, with 0 indicating complete compliance with physical laws and 1 indicating a complete anomaly.

[0151] A prior distribution is constructed based on the initial abnormal state prediction value. The prior distribution uses a Beta distribution, with parameters α and β set based on the known ratio of normal to abnormal samples in historical data. In practice, parameters α = 3 and β = 12 are empirically chosen to bias the distribution towards normal samples.

[0152] The Markov Chain Monte Carlo method is used to perform stratified repeated sampling of the prior distribution, with 500 sample points drawn each time. For each sample point, the node weight and edge probability are updated according to the topological structure of the state transition diagram. The update rule is based on the Bayesian posterior probability calculation method and is adjusted based on the consistency between the current sample and historical observation data. Specifically, if the state represented by a node is consistent with the historical meteorological patterns of the same time period and region, the node weight is increased; if the state change represented by a transition edge exceeds the physically reasonable range, the edge probability is reduced.

[0153] After 10 rounds of iterative optimization, the node weights and edge probabilities of the state transition graph gradually converged, ultimately obtaining the posterior distribution of the parameters for the optimal parameter configuration. This posterior distribution reflects the likelihood of different parameter values ​​and provides a stable and reliable parameter estimate for the state transition graph.

[0154] Using the Markov Chain Monte Carlo method, 200 stratified samplings were performed on the posterior distribution of the parameters. Each sampling step generated a set of parameter configurations based on the posterior distribution. These parameter configurations were combined with the optimal parameter configuration to construct a new state transition diagram. This new state transition diagram also contains eight nodes, but the node weights and edge probabilities have been optimized to better reflect the distribution characteristics of actual meteorological data.

[0155] The anomaly score is calculated by comparing the node weights and edge probabilities of the new state transition graph with the original state transition graph. This difference is calculated using the relative entropy method, namely the KL divergence between the two distributions. If the two graphs differ significantly, there is a significant deviation between the original and optimized graphs, indicating that the corresponding data is anomaly. If the difference is small, it indicates that the data conforms to physical laws. In practice, KL divergence values ​​below 0.05 are generally considered normal, those between 0.05 and 0.15 are slightly anomaly, and those above 0.15 are marked as significantly anomaly.

[0156] The probability density of the node distribution in the new state transition graph is estimated using kernel density estimation, a Gaussian kernel function, and a bandwidth parameter of h = 0.08. The estimated result reflects the distribution of the weights of each node, with the peak representing the weight value and the width of the distribution reflecting the degree of uncertainty.

[0157] The Markov Chain Monte Carlo method was used to analyze the fluctuation characteristics of node weights and edge probabilities. The standard deviation was calculated from 200 sampling results to form an uncertainty index. The fluctuation characteristic analysis found that the standard deviation of node weights corresponding to normal data is generally below 0.02, while abnormal data corresponds to higher volatility.

[0158] Combined with the confidence intervals for the physical constraints, a comprehensive uncertainty measure is calculated. Confidence intervals are statistically derived from historical data. For example, the 95% confidence interval for the daily temperature range is [-10°C, +8°C]. If the data falls within this interval, the certainty score is increased; otherwise, the uncertainty measure is increased. The final uncertainty measure ranges from 0 to 1, where 0 indicates complete certainty and 1 indicates extremely high uncertainty.

[0159] Based on the correlation between the numerical distribution of anomaly scores and the topological structure of the new state transition diagram, a segmented mapping relationship is established. This mapping relationship divides the anomaly scores into multiple intervals, each corresponding to a quality level. The intervals are divided as follows: scores of 0-0.05 correspond to "high-quality" data, 0.05-0.15 correspond to "good" data, 0.15-0.3 correspond to "suspicious" data, 0.3-0.5 correspond to "error" data, and 0.5 or above corresponds to "severe error" data.

[0160] The threshold for segment-by-segment mappings is adaptively adjusted based on the uncertainty metric. When the uncertainty metric is high (e.g., exceeding 0.4), the threshold is adjusted downward by 10%, making the system more inclined to mark data as lower quality. When the uncertainty metric is low (e.g., below 0.1), the threshold is adjusted upward by 5%, allowing for a more relaxed judgment criteria. This adaptive adjustment mechanism dynamically optimizes quality assessment criteria based on the degree of data certainty.

[0161] After these processing steps, the quality level of the data being tested is finally determined. This level includes both the quality level and the confidence level. For example, "Good (0.92)" indicates that the data has been rated as good with a 92% confidence level. This quality level is fed back to the LSTM-AE model for subsequent model training and optimization, forming a closed-loop optimization mechanism.

[0162] In an actual application case, a set of temperature data reported by a meteorological station showed a sudden drop from 23°C to 5°C within 1 hour. The system found through state transition diagram analysis that this change caused an abnormal state score of 0.42 and an uncertainty measure of 0.18. It was eventually marked as "Error (0.86)", indicating that this was abnormal data caused by sensor failure or data transmission error.

[0163] In an optional embodiment, the results of the multiple stratified samplings are combined with the optimal parameter configuration to construct a new state transition diagram, and the abnormal state score is obtained by comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram.

[0164] Constructing a node distribution probability matrix based on the results of the multiple stratified samplings, performing a convolution operation on the node distribution probability matrix and the optimal parameter configuration to obtain node weights, constructing a node distribution of a state transition diagram based on the node distribution probability matrix, performing Fourier transform on the node distribution to obtain frequency domain features, optimizing the node distribution probability matrix based on the frequency domain features to obtain edge probabilities, and combining the node weights and the edge probabilities to construct a new state transition diagram;

[0165] The node weights of the new state transition diagram and the original state transition diagram are subjected to the Fourier transform to obtain the node weight frequency domain features, the log likelihood ratio of the node weight frequency domain features is calculated to obtain the node weight difference value, the edge probability of the new state transition diagram and the original state transition diagram is subjected to the Fourier transform to obtain the edge probability frequency domain features, the log likelihood ratio of the edge probability frequency domain features is calculated to obtain the edge probability difference value, and the node weight difference value and the edge probability difference value are weightedly fused to obtain the abnormal state score.

[0166] A node distribution probability matrix is ​​constructed based on the results of multiple stratified samplings. Specifically, suppose a network system undergoes 10 stratified samplings, each collecting status data for five key nodes. After compiling these sampling results, a 10×5 node distribution probability matrix P is obtained, where P[i][j] represents the probability distribution value of the jth node in the i-th sampling. For example, the values ​​of P are: [[0.12, 0.23, 0.18, 0.27, 0.20], [0.15, 0.21, 0.19, 0.25, 0.20], ...].

[0167] The node weights are obtained by convolving the node distribution probability matrix with the optimal parameter configuration. Assume that the optimal parameter configuration is a weight vector W = [0.3, 0.2, 0.15, 0.25, 0.1], which represents the importance weights of the five key indicators. For each row Pi in the node distribution probability matrix P, the convolution of Pi and W is calculated to obtain the node weight NW. In specific implementation, for each node j, all sampled values ​​of the node in P are multiplied by the corresponding weight and summed, that is, NW[j] = sum(P[i][j] × W[j]), where i ranges from 1 to 10. For example, for node 1, its weight is calculated as (0.12 + 0.15 + ...) / 10 × 0.3 = 0.135.

[0168] The node distribution of the state transition graph is constructed based on the node distribution probability matrix. The node distribution probability matrix P is converted into the distribution of nodes in the state space. Specifically, for each node j, its probability distribution values ​​across all samples are collected to form a distribution vector ND[j]. For example, the distribution vector for node 1 is [0.12, 0.15, ...], representing the distribution of the node across 10 samples.

[0169] Perform a Fourier transform on the node distribution to obtain frequency domain features. Apply a discrete Fourier transform to each node's distribution vector ND[j] to convert it into frequency domain features FD[j]. The purpose of this step is to extract the periodic characteristics of the node distribution. In practice, this can be achieved using the Fast Fourier Transform algorithm. For example, performing a Fourier transform on the distribution vector [0.12, 0.15, ...] for node 1 yields the frequency domain features [0.135, 0.02-0.01i, 0.005+0.002i, ...].

[0170] The node distribution probability matrix is ​​optimized based on frequency domain features to obtain edge probabilities. The original node distribution probability matrix is ​​optimized using the frequency domain features FD to calculate the transition probabilities between nodes, namely the edge probabilities EP. Specifically, for any two nodes j and k, the correlation between their frequency domain features is calculated and combined with the original probability distribution values ​​to generate edge probabilities EP[j][k]. For example, the edge probability from node 1 to node 2 is 0.35, indicating a 35% probability of transitioning from state 1 to state 2.

[0171] The node weights and edge probabilities are combined to construct a new state transition graph. Using the calculated node weights NW and edge probabilities EP, a new state transition graph G is constructed. new In this graph, the weight of each node j is NW[j], and the weight of the edge from node j to node k is EP[j][k]. For example, G new The weight of node 1 is 0.135, and the weight of the edge from node 1 to node 2 is 0.35.

[0172] Perform Fourier transform on the node weights of the new state transition graph and the original state transition graph to obtain the frequency domain characteristics of the node weights. Assume that the original state transition graph G orig The node weights in the new graph are OW, and Fourier transforms are performed on NW and OW respectively to obtain the frequency domain features FNW and FOW. For example, for node 1, its weight 0.135 in the new graph is [0.135, 0, 0, ...] after Fourier transform, while its weight 0.15 in the original graph is [0.15, 0, 0, ...] after Fourier transform.

[0173] Calculate the log-likelihood ratio of the frequency domain features of node weights to obtain the node weight difference. For each node j, calculate the log-likelihood ratio of its frequency domain features in the old and new graphs, i.e., ND[j] = log(FNW[j] / FOW[j]). This step quantifies the degree of change in node weights in the frequency domain. For example, the weight difference of node 1 is log(0.135 / 0.15) = -0.105.

[0174] Perform Fourier transform on the edge probability of the new state transition graph and the original state transition graph to obtain the edge probability frequency domain feature. new and G orig The edge probabilities EP and OEP in the graph are Fourier transformed to obtain the frequency domain features FEP and FOEP. For example, the probability of an edge from node 1 to node 2 in the new graph is 0.35, after Fourier transform, it becomes [0.35, 0, 0, ...], while the probability of 0.4 in the original graph is [0.4, 0, 0, ...].

[0175] Calculate the log-likelihood ratio of the frequency-domain features of the edge probability to obtain the edge probability difference. For each pair of nodes (j, k), calculate the frequency-domain log-likelihood ratio of their edge probabilities, i.e., ED[j][k] = log(FEP[j][k] / FOEP[j][k]). For example, the edge probability difference between node 1 and node 2 is log(0.35 / 0.4) = -0.134.

[0176] The node weight differences and edge probability differences are weighted and combined to produce an abnormality score. The importance weights alpha and beta are set for the node weight differences and edge probability differences (alpha + beta = 1), and the overall score S is calculated as S = alpha × sum(ND[j]) + beta × sum(ED[j][k]). For example, if alpha = 0.6 and beta = 0.4, the sum of the node differences is -0.5, and the sum of the edge differences is -0.7, then the anomaly score S = 0.6 × (-0.5) + 0.4 × (-0.7) = -0.58. A larger score indicates a more significant difference between the old and new state transition graphs, and a higher probability that the system is in an abnormal state. In practical applications, a threshold T can be set. When S > T, the system is considered to be in an abnormal state and the corresponding alarm mechanism is triggered.

[0177] Through the above method, accurate monitoring of the system status is achieved, which can effectively capture abnormal changes in the state transition process and improve the security and stability of the system.

[0178] Figure 5This is a comparison chart of the accuracy of abnormal state scoring in an embodiment of the present invention, which shows the comparison curves of the evaluation performance of three different methods at different noise levels. Among them, the traditional feature comparison method represents a basic algorithm based on traditional feature extraction and comparison, with an initial performance of 90.0%. As the noise level increases, the performance gradually decreases, dropping to 84.2% at a noise level of 5%, 79.8% at 10%, 74.4% at 15%, 68.5% at 20%, 63.2% at 25%, and finally to 58.7% when the noise reaches 30%; the single parameter evaluation method represents a method that uses a single feature parameter for evaluation, with an initial performance of 93.8%, and the performance at different noise levels is 88.1 and 88.7% respectively. % (5%), 83.5% (10%), 77.3% (15%), 71.8% (20%), 67.4% (25%), and 62.6% (30%). This technical solution represents a new type of comprehensive evaluation algorithm, showing the strongest noise resistance, with an initial performance of up to 97.6%. Its performance degrades the least as the noise level increases, reaching 95.3% (5%), 91.2% (10%), 87.7% (15%), 81.3% (20%), 75.6% (25%), and 69.8% (30%), respectively. Overall, this technical solution maintains optimal performance at all noise levels and exhibits greater robustness than the other two methods. In particular, it maintains a high evaluation accuracy in high-noise environments (noise levels > 20%), fully demonstrating the superiority of this solution in complex noise environments.

[0179] The method further comprises:

[0180] During the data preprocessing phase, we collected multi-year meteorological observation data from national meteorological stations and numerical forecast data at the corresponding stations. Taking temperature observations from national stations from 2018 to 2022 as an example, we selected hourly temperature observations and simultaneous short-term temperature forecasts from 100 national stations across China. For the numerical forecast data, we extracted the forecast values ​​for the station locations from the original forecast fields using bilinear interpolation. Data normalization used mean normalization, which involves subtracting the mean from the original data and dividing it by the standard deviation. For the temperature data, the calculated mean was approximately 15.3°C and the standard deviation was approximately 10.2°C. After normalization, the temperature data were mostly distributed within the range [-3, 3].

[0181] Based on the characteristics of time series, each 24-hour observation data was organized into a sequence sample. Each sample contained the observed values ​​of four meteorological variables: temperature, air pressure, humidity, and wind speed, as well as the temperature forecast. After removing sequences with missing measurements, a total of approximately 150,000 valid sequence samples were obtained. The first 80% (approximately 120,000 samples) were used for model training in chronological order, and the remaining 20% ​​(approximately 30,000 samples) were used for model validation.

[0182] During the model architecture design phase, a dual-encoder LSTM-AE structure was constructed. The first encoder processes the multivariate observation sequence and consists of two stacked LSTM layers, each with 128 hidden units. The input is a tensor of shape [24, 4], representing the observations of four meteorological variables within a 24-hour period. The second encoder processes the temperature forecast sequence and consists of a single LSTM layer with 64 hidden units. The input is a tensor of shape [24, 1], representing the temperature forecast value within a 24-hour period. The outputs of the two encoders are combined into a 192-dimensional fused representation through a concatenation layer. The decoder uses two stacked LSTM layers, each with 128 hidden units, to reconstruct the fused representation into the original observation sequence.

[0183] Within the LSTM unit, the forget gate, input gate, and output gate each control the retention of historical information, the influence of current input, and the selection of output content through different weight matrices. Taking the first LSTM layer of the first encoder as an example, with an input feature dimension of 4 and a hidden state dimension of 128, the forget gate's weight matrix dimension is [132, 128], containing approximately 17,000 parameters. Similarly, the input and output gates each have the same number of parameters. These parameters are continuously optimized during training, enabling the model to effectively capture the temporal dependencies of meteorological data.

[0184] During model training and hyperparameter optimization, the model was trained using mini-batch gradient descent with a batch size of 64, an initial learning rate of 0.001, and the Adam optimizer. The mean squared error (MSE) loss function was used, calculating the mean squared difference between the observed and reconstructed values ​​for the target quality control variable (temperature). During training, the learning rate was halved when the loss on the validation set stopped decreasing for five consecutive epochs. The total number of training epochs was set to 100, but the model typically achieved optimal performance between 50 and 70 epochs.

[0185] A grid search was performed on the key hyperparameter, LSTM hidden unit dimension, from 32 to 256 (in powers of 2 steps). The results showed that when the first encoder hidden unit was 128 and the second encoder hidden unit was 64, the MSE on the validation set reached a minimum of 0.0042. The model has approximately 300,000 total parameters and takes about two hours to train on a standard GPU.

[0186] During the anomaly score calculation and threshold determination phase, the reconstruction error is calculated for each sample in the validation set as the anomaly score. Statistical analysis shows that the reconstruction error for normal temperature data is mostly below 0.05, while the reconstruction error for anomaly data is typically above 0.1. By calculating the percentiles of the anomaly scores in the validation set, four thresholds are obtained: the 0.5 percentile is approximately 0.025, the 0.05 percentile is approximately 0.085, the 0.005 percentile is approximately 0.135, and the 0.0005 percentile is approximately 0.175. These four thresholds correspond to different levels of anomaly detection criteria.

[0187] In the model application stage, the temperature observation of a certain station on January 5, 2023 is used as an example. The 24-hour temperature observation sequence of that day is [-5.2, -5.8, -6.3, -6.5, -6.2, -6.0, -5.8, -5.1, -4.3, -3.2, -2.1, -1.0, 0.2, 0.5, 0.3, -0.6, -1.8, -2.5, -3.1, -3.8, -4.2, -4.6, -4.9, -5.1] ℃, and the corresponding forecast series is [-5.5,-6.0,-6.2,-6.3,-6.1,-5.9,-5.5,-4.9,-4.0,-3.0,-1.8,-0.8,0.0,0.4,0.1,-0.8,-2.0,-2.8,-3.4,-4.0,-4.4,-4.8,-5.0,-5.3]℃. The normalized observation and forecast sequences were fed into the trained model, resulting in a reconstructed temperature sequence of [-5.3, -5.7, -6.2, -6.4, -6.3, -6.1, -5.7, -5.0, -4.2, -3.1, -2.0, -0.9, 0.1, 0.4, 0.2, -0.7, -1.9, -2.6, -3.2, -3.9, -4.3, -4.5, -5.0, -5.2]°C. The mean square error (MSE) between the original and reconstructed temperature observations was 0.0038, below the most relaxed threshold of 0.025. Therefore, the observation sequence was deemed normal.

[0188] If anomalies are encountered, such as a sudden temperature change caused by a sensor failure, such as an unusual jump to 10.2°C at the 13th hour of the sequence, the reconstructed temperature at that moment will remain around 0.1°C, resulting in a significant reconstruction error. The calculated MSE in this case is approximately 0.153, exceeding the 0.005 percentile threshold of 0.135 but below the 0.0005 percentile threshold of 0.175, resulting in a more stringent anomaly warning.

[0189] During the ongoing operation phase, observers can adjust the anomaly detection threshold based on actual conditions and empirical feedback. For example, in cold winter regions with relatively gradual temperature fluctuations, the threshold can be appropriately lowered to increase detection sensitivity; whereas in spring and autumn regions with large daily temperature fluctuations, the threshold can be appropriately raised to avoid false alarms. This dynamic adjustment mechanism, combined with human experience, continuously optimizes the model's quality control performance under different meteorological conditions.

[0190] Compared to traditional statistical methods, this approach can simultaneously consider the interrelationships and temporal characteristics of multiple variables and integrate numerical forecast information, significantly improving the accuracy of anomaly detection. In actual application tests, the detection rate of artificially introduced anomalous data reached over 95%, while the false alarm rate was kept below 5%, meeting the meteorological department's high data quality control requirements.

[0191] According to a second aspect of an embodiment of the present invention, an electronic device is provided, including:

[0192] processor;

[0193] a memory for storing processor-executable instructions;

[0194] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0195] According to a third aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0196] The present invention may be a method, an apparatus, a system and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for executing various aspects of the present invention.

[0197] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. The AI ​​meteorological station data quality control method integrated with numerical forecast is characterized by: include: Acquiring multivariate meteorological observation data and numerical forecast data from a meteorological station, and performing standardization processing on the multivariate meteorological observation data and the numerical forecast data to form an observation time series dataset and a numerical forecast time series dataset; Constructing a physical driving model based on the LSTM-AE structure, extracting the physical relationships and temporal evolution characteristics between meteorological elements in the observation time series dataset through a long short-term memory network encoder based on the physical driving model, parsing the observation time series dataset and the numerical forecast time series dataset to obtain observed physical characteristics and predicted physical characteristics; The observed physical features and the forecast physical features are processed based on the LSTM decoder in the physical drive model, and feature reconstruction is completed according to the physical laws and spatiotemporal evolution characteristics of meteorological elements to generate a reconstructed observation time series, including: Performing initial reconstruction of the observed physical features and the predicted physical features based on the LSTM decoder to generate initial reconstruction features, calculating a reconstruction error between the initial reconstruction features and the observed time series, and combining the reconstruction error with a hidden state to obtain an error propagation feature; Combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term; combining the dynamic error compensation term with the initial reconstructed feature to obtain a corrected reconstructed feature, calculating the meteorological physical correlation between air pressure, temperature, humidity, wind speed, and precipitation in the corrected reconstructed feature, and comparing the meteorological physical correlation with the physical correlation between corresponding meteorological elements in the observation time series to obtain a deviation value between the corrected reconstructed feature and the observation time series in terms of physical laws; Adjusting the forgetting weight, input weight, and output weight in the LSTM decoder according to the bias value, updating the cell state and hidden state using the adjusted weights, and performing a new round of reconstruction on the observed physical characteristics and the predicted physical characteristics based on the updated cell state and hidden state to generate a reconstructed observation time series that conforms to the physical laws and spatiotemporal evolution characteristics of meteorological elements; The observation time series dataset is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation indicator to complete the training and parameter optimization of the LSTM codec in the physical drive model; Based on the physical driving model, an abnormal state evaluation is performed on the data to be detected, and a quality grade label is determined according to the result of the abnormal state evaluation. The labeling result is fed back to optimize the feature extraction capability and anomaly detection performance of the LSTM-AE structure, thereby realizing quality control and continuous optimization of meteorological station data.

2. The method according to claim 1, characterized in that A physical driving model based on the LSTM-AE structure is constructed. Based on the physical driving model, the physical relationship and temporal evolution characteristics between meteorological elements in the observation time series dataset are extracted through a long short-term memory network encoder. The observation time series dataset and the numerical forecast time series dataset are parsed to obtain the observed physical characteristics and the predicted physical characteristics, including: Calculating the observation data of the observation time series data set through a long short-term memory network encoder to obtain an observation mean vector and an observation variance vector, constructing an observation feature normal distribution based on the observation mean vector and the observation variance vector, and sampling from the observation feature normal distribution to obtain a latent variable of the observation data; Introducing a priori normal distribution constraints on the network parameters of the physical driving model, and calculating the physical association weights of meteorological elements at different time scales through a Bayesian neural network based on the prior normal distribution constraints; Performing weighted reconstruction on the latent variables of the observation data according to the physical association weights to obtain a reconstructed observation feature vector, constructing a conditional probability distribution based on the reconstructed observation feature vector, and performing multiple sampling from the conditional probability distribution to obtain an observed physical feature; The numerical forecast time series data set is input into the long short-term memory network encoder to obtain a forecast mean vector and a forecast variance vector of the numerical forecast time series data set, a forecast feature normal distribution is constructed based on the forecast mean vector and the forecast variance vector, latent variables of the forecast data are sampled from the forecast feature normal distribution, and the latent variables of the forecast data are weightedly reconstructed based on the physical association weights to obtain forecast physical characteristics that characterize the physical relationship of meteorological elements.

3. The method according to claim 1, characterized in that Combining the error propagation feature and the initial reconstruction feature to calculate a compensation coefficient, and correcting the initial reconstruction feature based on the compensation coefficient to obtain a dynamic error compensation term includes: Mapping the error propagation feature and the initial reconstruction feature to a feature space of the same dimension to obtain a projected error propagation feature and a projected initial reconstruction feature; dividing the projected initial reconstruction feature into a plurality of local regions, calculating a mean and a standard deviation for each of the local regions, constructing a local feature statistic based on the mean and the standard deviation, and combining the local feature statistic with the projected error propagation feature to generate a local compensation coefficient; Calculating a correlation matrix between the projected initial reconstructed features and the projected error propagation features, calculating a global feature correlation based on the correlation matrix, and generating a global compensation coefficient according to the global feature correlation; Performing feature splicing on the projected initial reconstructed features and the projected error propagation features, generating a weight coefficient according to the result of the feature splicing, and performing a weighted combination of the local compensation coefficient and the global compensation coefficient based on the weight coefficient to obtain a final compensation coefficient; The final compensation coefficient is subjected to element-wise multiplication operation with the projected error propagation feature, including: taking each dimension feature of the final compensation coefficient as a weight, weighting the feature of the corresponding dimension of the projected error propagation feature, combining the weighted features of each dimension to obtain a compensation feature, fusing the compensation feature with the initial reconstruction feature to generate a corrected reconstruction feature, and obtaining a dynamic error compensation term according to the difference between the corrected reconstruction feature and the initial reconstruction feature.

4. The method according to claim 1, wherein The observation time series data set is divided into a training set and a validation set, and the deviation between the reconstructed observation time series and the original observation time series is used as an evaluation index, including: Performing feature reconstruction on the observation time series data set to obtain a reconstructed observation time series, wherein the reconstructed observation time series has the same temporal distribution characteristics as the original observation time series; The timing deviation between the reconstructed observation time series and the original observation time series at corresponding moments is calculated, and the timing deviation is used as an evaluation index to evaluate the reconstruction effect.

5. The method according to claim 1, wherein Performing an abnormality evaluation on the data to be detected based on the physical drive model, and determining a quality grade mark according to a result of the abnormality evaluation includes: Input the data to be detected into the physical drive model to generate a state feature vector, construct a state transition diagram based on the state feature vector, and generate an initial abnormal state prediction value by analyzing the matching degree between the state transition diagram and the physical constraint conditions; A prior distribution is constructed based on the initial abnormal state prediction value, and the Markov chain Monte Carlo method is used to perform stratified repeated sampling on the prior distribution. The node weights and edge probabilities of the state transition diagram are iteratively updated according to the sampling results, and the parameter posterior distribution of the optimal parameter configuration of the state transition diagram is obtained through multiple rounds of iterative optimization; Performing multiple stratified sampling on the parameter posterior distribution using the Markov Chain Monte Carlo method, combining the results of the multiple stratified sampling with the optimal parameter configuration to construct a new state transition diagram, and obtaining an abnormal state score by comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram; Performing a probability density estimation on the node distribution in the new state transition graph, analyzing the fluctuation characteristics of the node weights and edge probabilities based on the Markov chain Monte Carlo method, and calculating an uncertainty measure based on the fluctuation characteristics combined with the confidence interval of the physical constraints; A segmented mapping relationship is established according to the correlation between the numerical distribution of the abnormal state score and the topological structure of the new state transition diagram, and a threshold of the segmented mapping relationship is adaptively adjusted based on the uncertainty metric to determine a quality level mark.

6. The method according to claim 5, characterized in that Combining the results of the multiple stratified samplings with the optimal parameter configuration to construct a new state transition diagram, and comparing the differences in node weights and edge probabilities between the new state transition diagram and the original state transition diagram to obtain an abnormal state score includes: Constructing a node distribution probability matrix based on the results of the multiple stratified samplings, performing a convolution operation on the node distribution probability matrix and the optimal parameter configuration to obtain node weights, constructing a node distribution of a state transition diagram based on the node distribution probability matrix, performing Fourier transform on the node distribution to obtain frequency domain features, optimizing the node distribution probability matrix based on the frequency domain features to obtain edge probabilities, and combining the node weights and the edge probabilities to construct a new state transition diagram; The node weights of the new state transition diagram and the original state transition diagram are subjected to the Fourier transform to obtain the node weight frequency domain features, the log likelihood ratio of the node weight frequency domain features is calculated to obtain the node weight difference value, the edge probability of the new state transition diagram and the original state transition diagram is subjected to the Fourier transform to obtain the edge probability frequency domain features, the log likelihood ratio of the edge probability frequency domain features is calculated to obtain the edge probability difference value, and the node weight difference value and the edge probability difference value are weightedly fused to obtain the abnormal state score.

7. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 6.

8. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 6 is implemented.

Images