Atmospheric environment monitoring fixed station pollutant monitoring data backtracking method

Abstract

The invention belongs to the field of atmospheric environment monitoring, and relates to an atmospheric environment monitoring fixed station pollutant monitoring data backtracking method. According to the method, a BPCA-LSTM model combining Bayesian Principal Component Analysis (BPCA) and a Long Short-Term Memory (LSTM) is constructed, an available complete data set is selected and formed from pollutant concentration monitoring historical data of all N stations according to the range of pollutant concentration monitoring data needing to be backtracked by M stations, and the available complete data set is selected from the pollutant concentration monitoring historical data of all N stations according to the range of the pollutant concentration monitoring data needing to be backtracked by M stations. And performing data cleaning and data dimension reduction on the selected data set, and learning a concentration time sequence data mapping relationship between N-M sites and M sites. And through the model, data of N-M sites corresponding to the data missing period of the M sites are used for backtracking missing data of the M sites. The method is simple and universal, provides a reliable data backtracking scheme for atmospheric environment monitoring, and has important practical significance and wide application prospects.

Description

Technical Field

[0001] The invention belongs to the field of atmospheric environment monitoring and relates to a pollutant monitoring data backtracking method for atmospheric environment monitoring fixed sites. Background Art

[0002] With the continuous advancement of urbanization and the transformation of the economy to the post-industrial stage, the importance of green travel in modern urban governance has become increasingly prominent, and it is an important way to achieve urban ecological civilization construction. The promotion and implementation of green travel depends on the assessment and forecast of air quality. Complete pollutant monitoring data from fixed atmospheric environmental monitoring stations is an important data basis for air quality assessment and forecasting, as well as for capturing the causes, evolution and dissipation laws of concentration. The lack of such data not only affects the comprehensive assessment of the urban atmospheric environment, but also restricts the accurate location of pollution sources and the formulation of atmospheric pollution control measures.

[0003] In the actual monitoring process, due to equipment failure, maintenance, additions and reductions and other reasons, data loss at monitoring sites is inevitable. There are two main solutions to the missing data of atmospheric pollutant points. One is to use pollution diffusion models such as AERMOD, CALINE 4, and ADMS to simulate pollution diffusion. Such models rely on a complete pollution emission inventory, which is time-consuming to calculate and cannot guarantee the accuracy of the simulation results. The other is to use spatiotemporal interpolation methods, Bayesian methods, and deep learning methods to complete the data from a data perspective. The research mainly focuses on the block scale, and the pollution situation in the area is analyzed based on relatively dense data sampling. At present, there is still a lack of a method that can effectively trace data at the urban scale, which makes the comprehensive assessment of pollution at the mesoscale challenging. Summary of the invention

[0004] In order to solve the above problems, the present invention provides a method for backtracking pollutant monitoring data at fixed sites for atmospheric environment monitoring, which can learn the concentration time series data mapping relationship between the data of NM sites other than the M monitoring sites with missing data and the data of the M monitoring sites based on a method combining statistics and deep learning, and then use the data of the NM sites corresponding to the period when the data of the M sites were missing to backtrack and complete the missing data of the M sites.

[0005] Technical solution of the present invention:

[0006] A method for tracing back pollutant monitoring data of a fixed site for atmospheric environment monitoring comprises the following steps:

[0007] S1. Construct a BPCA-LSTM network model, including a Bayesian principal component analysis (BPCA) part and a long short-term memory (LSTM) network part; at the same time, according to the range of pollutant concentration monitoring data that needs to be traced back at M stations, select and form a usable complete data set from the historical data of pollutant concentration monitoring at all N stations;

[0008] S2. Perform data cleaning on the selected data set, including pollutant concentration measurement standardization, outlier removal and missing point filling, and apply the Bayesian principal component analysis part to the cleaned data for data dimensionality reduction;

[0009] S3. Set the initial data sliding window and BPCA-LSTM network model hyperparameters, and use the long short-term memory network part to construct the time series pollutant concentration data mapping relationship between the NM stations and each of the M stations;

[0010] S4. Optimize and adjust the sliding window period length and the hyperparameters of the BPCA-LSTM network model;

[0011] S5. Taking the historical data of NM stations corresponding to the period when the data of M stations were missing as input, the missing pollution concentration monitoring data of M stations were backtracked through the constructed BPCA-LSTM network model.

[0012] The present invention relies on Bayesian principal component analysis (BPCA) and long short-term memory network (LSTM) technology to realize the backtracking of pollutant monitoring data of fixed sites for atmospheric environment monitoring. Through the BPCA-LSTM model, the concentration time series data mapping relationship between the data-missing monitoring site and other data-complete sites is learned to complete the missing data. This method effectively solves the problem of missing monitoring data caused by equipment failure, maintenance, increase and decrease, fills the gap in the backtracking method of urban scale monitoring data, provides a reliable solution for the government and environmental protection departments, supports more scientific air quality assessment and pollution control, promotes the development of green travel and the construction of ecological civilization, and has important practical significance and broad application prospects.

[0013] Beneficial effects of the present invention:

[0014] Achieve accurate data backtracking: By combining statistics and deep learning methods, the temporal relationship patterns of atmospheric pollutant concentrations can be effectively captured and missing data can be accurately traced back.

[0015] Filling the technical gap: It fills the technical gap in the backtracking of atmospheric environment monitoring data at the urban scale, which has important practical significance and application prospects.

[0016] Wide applicability: It is suitable for backtracking of different types of pollutant monitoring data, has strong versatility, and can be used in a variety of urban environmental monitoring application scenarios.

[0017] Improve work efficiency: Use data cleaning and dimensionality reduction techniques (such as BPCA) to pre-process data, reduce data dimensions and noise, and optimize sliding windows and network model hyperparameters through proxy models to reduce computing time, ensure data backtracking is completed within a limited time, and improve efficiency.

[0018] Support scientific decision-making: Provide reliable pollutant monitoring data traceback results, provide the government and environmental protection departments with a basis for air quality assessment and air pollution control, support scientific decision-making, promote the development of green travel, and help build urban ecological civilization. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 This is the main process of a method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring of the present invention;

[0020] Figure 2 A schematic diagram of a weight matrix for data dimensionality reduction using BPCA in step 2 of the present invention;

[0021] Figure 3 The present invention uses complete period data to construct a mapping relationship and traces back the overall process of missing period data;

[0022] Figure 4 Schematic diagram of the bidirectional LSTM structure used in step 3 of the present invention;

[0023] Figure 5 It is a scatter plot of true value-predicted value on the training set of the optimal model obtained by the present invention in the embodiment;

[0024] Figure 6 It is the site data backtracking result obtained by relying on the present invention in the embodiment. DETAILED DESCRIPTION

[0025] The technical solution provided by the present application will be further described below in conjunction with specific embodiments and accompanying drawings. The advantages and features of the present application will become more apparent with the following description.

[0026] like Figure 1 As shown, a method for tracing back pollutant monitoring data of a fixed site for atmospheric environment monitoring comprises the following steps:

[0027] S1. Construct a BPCA-LSTM network model, including a Bayesian principal component analysis (BPCA) part and a long short-term memory (LSTM) network part; at the same time, according to the range of pollutant concentration monitoring data that needs to be traced back at M stations, select and form a usable complete data set from the historical data of pollutant concentration monitoring at all N stations;

[0028] S2. Perform data cleaning on the selected data set, including pollutant concentration measurement standardization, outlier removal and missing point filling, and apply the Bayesian principal component analysis part to the cleaned data for data dimensionality reduction;

[0029] S3. Set the initial data sliding window and BPCA-LSTM network model hyperparameters, and use the long short-term memory network part to construct the time series pollutant concentration data mapping relationship between the NM stations and each of the M stations;

[0030] S4. Optimize and adjust the sliding window period length and the hyperparameters of the BPCA-LSTM network model;

[0031] S5. Taking the historical data of NM stations corresponding to the period when the data of M stations were missing as input, the missing pollution concentration monitoring data of M stations were backtracked through the constructed BPCA-LSTM network model.

[0032] Furthermore, the data selected from the historical data of all N stations in step S1 are generally unified as the pollutant concentration monitoring data of N stations in the next 2 to 3 years adjacent to the pollutant concentration data missing interval of the M stations that need to trace back data, so as to ensure that the training and test data in step S3 are sufficient, and the relationship pattern of the concentration time series data between the N stations covered by the selected data can effectively characterize the concentration time series data mapping relationship between the NM stations and the M stations in the data missing interval of the M stations. The above data needs to discard the station data with data integrity less than 90% in the NM stations, and the data integrity is defined by the following formula:

[0033]

[0034] Among them, Ω represents the completeness of the data, Q represents the total amount of data that should be monitored in the data set, and Q a Indicates the total amount of data that has not been monitored in the dataset.

[0035] Furthermore, in step S2, the selected data set is cleaned and dimensionally reduced. Data cleaning includes standardization of pollutant concentration, removal of outliers and filling in missing points, as follows:

[0036] (a3.1) For pollutant concentration measurement standardization, if the stored data uses the concentration value under standard conditions (273.15K, 101.325kPa), it is converted into the concentration value under actual conditions. Considering the small difference in atmospheric pressure, the simplified conversion formula is as follows:

[0037]

[0038] Among them, C r Indicates the concentration value under actual conditions, C s Indicates the concentration value under standard conditions, T r Indicates the ambient temperature under actual conditions, which uses the Kelvin temperature scale;

[0039] (a3.2) For outlier removal, abnormal detection records (for example, data of 9999 and -9999) and data outliers are removed. Outliers are defined as values ​​less than the lower limit or greater than the upper limit. The upper and lower limit calculation formulas are:

[0040]

[0041] Among them, U and L represent the upper and lower limits respectively, Q3 and Q1 represent the upper quartile and lower quartile respectively, and IQR represents the interquartile range, that is, the distance between Q3 and Q1;

[0042] (a3.3) For missing point filling, the missing points of the original data are filled with the average concentration of pollutants within the selected data range. For the data removed due to outliers, if they were originally abnormal detection records, they are filled with the above-mentioned average concentration value. If they are outliers, those less than the lower limit are filled with the minimum value of the data after removing the outliers, and those greater than the upper limit are filled with the maximum value of the data after removing the outliers.

[0043] Next, Bayesian principal component analysis is used to reduce the dimension of the cleaned data. The dimension reduction process is completed using the expectation-maximization algorithm (EM algorithm). For the actual concentration vector at a certain moment D is the feature dimension, where x is the vector of actual concentrations of all stations except the station with missing data at that moment. The goal of dimensionality reduction is to find a low-dimensional representation and the corresponding weight matrix (like Figure 2 right) and a bias vector So that:

[0044] x=Wz+μ+ε (4)

[0045] Among them, ε is the noise term. There are the following prior distribution assumptions about x, z, and W, where represents Gaussian distribution, I is the identity matrix, σ 2 is the observation noise variance:

[0046]

[0047] The parameter set that needs to be estimated under the current assumption is θ = {W, μ, σ 2 ,α}, where μ is estimated using the following formula:

[0048]

[0049] Parameters W, σ 2 , α is iteratively updated through the EM algorithm. The EM algorithm estimates the model parameters by alternating between the E step (expectation step) and the M step (maximization step). The following are the specific steps of the EM algorithm:

[0050] (b3.1) Initialize parameters W, σ 2 , α;

[0051] (b3.2) In step E, the expectation and covariance of the posterior distribution of the latent variable z are calculated as follows:

[0052] Where M = W T W+σ 2 I (10)

[0053]

[0054] (b3.3) In the M step, the parameters W and σ are updated according to the expectation maximization. 2 , α, and its update formula is as follows:

[0055]

[0056] Where Tr(·) represents the trace of the matrix;

[0057] (b3.4) Repeat steps E and M until the model parameters converge, that is, the iterative change value of the parameters is less than the set threshold.

[0058] Furthermore, step S3 constructs a mapping relationship between the time series pollutant concentration data of each of the NM stations and the M stations, so as to use the data of the remaining NM stations to trace back the data of the M stations when the data of the M stations are missing. Specifically, for the current time T, a time window [Tn, T+n] is selected, where n is a hyperparameter for adjusting the window size, and the pollutant concentration data within the time window of the NM stations is used as input (Z T-n ,…,ZT ,…,Z T+n ), the concentration value of the backtracking site required at the current time T is output (Y T ), and use the long short-term memory network (LSTM) to construct the mapping relationship between the two under the selected network hyperparameters (such as Figure 3 As shown in Figure 1). This mapping relationship takes into account the lag of the impact of pollution diffusion, and uses the concentration information of other stations around the unknown data station from the past to the future relative to the retrospective time to infer the concentration value of the unknown data station at the retrospective time. The network model adopts a bidirectional structure, that is, a bidirectional LSTM (BidirectionalLSTM, Bi-LSTM network) is constructed, as shown in Figure 1. Figure 4 As shown, it is used to capture and memorize the concentration relationship pattern of time series data in the time window in two directions: past to future (forward time series) and future to past (reverse time series). The training and validation sets are divided into 5:1.

[0059] Furthermore, step S4 optimizes the hyperparameters, adopts a differentiated optimization strategy, and adjusts the network hyperparameters including the time window size, learning rate, regularization parameter, batch size, network layer, number of neurons, and loss function. In this method, mean square error is used as the loss function, the Adadelta optimizer is used to dynamically adjust the learning rate, and the L2 regularization parameter is used. The remaining hyperparameters include the time window size, batch size, network layer, and number of neurons. The hyperparameter space is defined as x, and the Gaussian Process (Gaussian Process, ) is used as a proxy model for hyperparameter optimization. The specific steps are as follows:

[0060] (a4.1) Initialization phase: define the hyperparameter space x, and randomly sample the initial hyperparameter point set A = {a1, a2, …, a n}, based on the initial hyperparameter point set, the Bi-LSTM network model is trained, and the validation set loss y={y1,y2,…,y n}, define the Gaussian process f(a):

[0061]

[0062] Among them, μ(a) represents the mean function, that is, k(a,a′) is a kernel function that characterizes the correlation between a and a′ in the hyperparameter space. Here, the radial basis function (RBF) is used, and the initial sampling set is used. Train the Gaussian process model f(a);

[0063] (a4.2) Calculate the predicted distribution. For candidate point a, the predicted mean μ(a) and the predicted variance σ 2The calculation formula of (a) is as follows:

[0064] μ(a)=k(a,A)(K(A,A)+σ 2 I) -1 y (16)

[0065] σ 2 (a) = k(a,a) - k(a,A)(K(A,A) + σ 2 I) -1 k(A,a) (17)

[0066] Where A is the matrix of evaluated points, y is the corresponding objective function value vector, K(A,A) is the covariance matrix of the kernel function on the evaluated points, k(a,A) is the covariance vector between the candidate point and the evaluated points, and k(a,a) is the covariance between the candidate point and itself;

[0067] (a4.3) Select and optimize the acquisition function and use the expected improvement (EI) as the acquisition function. The calculation formula is as follows:

[0068]

[0069] Among them, a + is the current optimal point, ξ is a non-negative exploration parameter used to control the trade-off between exploration and exploitation of the algorithm, Φ(·) and φ(·) are the cumulative distribution function and probability of the standard normal distribution, respectively. Find the acquisition function α on the surrogate model. EI (a) The largest point a * , that is:

[0070]

[0071] (a4.4) Update the model, for the hyperparameter a * , train the model to get the validation set loss y * , add the new sampling point to the evaluated point set to form a set {(a1,y1),(a2,y2),…,(a n ,y n ),(a * ,y * )}, retrain the Gaussian process model using the updated dataset;

[0072] (a4.5) Check the stopping condition to see whether the stopping condition (maximum number of iterations, time limit or expected objective function value) is met. If so, stop. Otherwise, repeat (a4.2) to (a4.4).

[0073] Furthermore, in step S5, the BPCA-LSTM model is constructed and compiled with the help of the PyTorch deep learning framework, and the model is trained under the hyperparameter optimization framework to obtain the iteratively optimal hyperparameter settings and the corresponding network model under the preset time limit. For the period when the data of M stations are missing, the model is applied to the pollutant concentration time series window data of the remaining NM stations as input to trace the pollutant concentration data of the M stations.

[0074] Example

[0075] Taking the 2018 NO2 missing data backtracking of Jing'an West Yan'an Road Traffic Station in Shanghai as an example, relying on the data backtracking method proposed in the present invention, using the complete full-station data from 2019 to 2021, the time series mapping relationship between the NO2 concentration data of 11 national control stations and Jing'an West Yan'an Road Traffic Station is constructed by training the network model, and then the data of the 11 national control stations in 2018 is used as input, and the missing data backtracking of Jing'an West Yan'an Road Traffic Station in 2018 is completed using the trained network model. The specific implementation steps are as follows:

[0076] S1. Select data from the historical data for 2 to 3 years after the data missing period to form training and validation data sets to ensure that the training and test data in step S3 are sufficient, and the relationship pattern of the concentration time series data between all stations covered by the selected data can effectively characterize the NO2 concentration data mapping relationship between the 11 national control stations and the Jing'an Yan'an West Road Traffic Station in the data missing interval of the Jing'an Yan'an West Road Traffic Station, including the following data (as shown in Table 1):

[0077] (1) Training set: NO2 standard condition data from 11 national monitoring stations and Jing'an Yan'an West Road Traffic Station in the second half of 2019, 2020, and 2021;

[0078] (2) Validation set: NO2 standard condition data from 11 national monitoring stations and Jing'an Yan'an West Road Traffic Station in the first half of 2019;

[0079] The above data should exclude the site data with data integrity less than 90%. The data integrity is defined by the following formula:

[0080]

[0081] Among them, Ω represents the completeness of the data, Q represents the total amount of data that should be monitored in the data set, and Q a Indicates the total amount of data in the dataset that has not been monitored.

[0082] Table 1 Training and validation data sets (only part of them)

[0083]

[0084]

[0085] S2. Perform data cleaning and data dimension reduction on the selected data set. Data cleaning includes pollutant concentration measurement standardization, outlier removal and missing point filling, as follows:

[0086] (a3.1) For pollutant concentration measurement standardization, if the stored data uses the concentration value under standard conditions (273.15K, 101.325kPa), it is converted into the concentration value under actual conditions. Considering the small difference in atmospheric pressure, the simplified conversion formula is as follows:

[0087]

[0088] Among them, C r Indicates the concentration value under actual conditions, C s Indicates the concentration value under standard conditions, T r Indicates the ambient temperature under actual conditions, which uses the Kelvin temperature scale. The conversion example is shown in Table 2 below;

[0089] Table 2 Pollutant concentration measurement standardization

[0090]

[0091]

[0092] (a3.2) For outlier removal, in addition to removing abnormal detection records (for example, data of 9999 and -9999), data outliers are also removed. Outliers are defined as values ​​less than the lower limit or greater than the upper limit. The upper and lower limits are calculated as follows:

[0093]

[0094] Among them, U and L represent the upper and lower limits respectively, Q3 and Q1 represent the upper quartile and lower quartile respectively, and IQR represents the interquartile range, that is, the distance between Q3 and Q1;

[0095] (a3.3) For missing point filling, the missing points of the original data are filled with the average concentration of pollutants within the selected data range. For the data removed due to outliers, if they were originally abnormal detection records, they are filled with the above-mentioned average concentration value. If they are outliers, those less than the lower limit are filled with the minimum value of the data after removing the outliers, and those greater than the upper limit are filled with the maximum value of the data after removing the outliers.

[0096] Next, Bayesian principal component analysis is used to reduce the dimension of the cleaned data. The dimension reduction process is completed using the expectation-maximization algorithm (EM algorithm). For the actual concentration vector at a certain moment D is the feature dimension, where x is the vector of actual concentrations of all stations except the station with missing data at that moment. The goal of dimensionality reduction is to find a low-dimensional representation and the corresponding weight matrix and a bias vector So that:

[0097] x=Wz+μ+ε (4)

[0098] Among them, ε is the noise term. There are the following prior distribution assumptions about x, z, and W, where represents Gaussian distribution, I is the identity matrix, σ 2 is the observation noise variance:

[0099]

[0100]

[0101] The parameter set that needs to be estimated under the current assumption is θ = {W, μ, σ 2 ,α}, where μ is estimated using the following formula:

[0102]

[0103] Parameters W, σ 2 , α is iteratively updated through the EM algorithm. The EM algorithm estimates the model parameters by alternating between the E step (expectation step) and the M step (maximization step). The following are the specific steps of the EM algorithm:

[0104] (b3.1) Initialize parameters Q, σ 2 , α;

[0105] (b3.2) In step E, the expectation and covariance of the posterior distribution of the latent variable z are calculated as follows:

[0106] Where M = W T W+σ 2 I (10)

[0107]

[0108] (b3.3) In the M step, the parameters Q and σ are updated according to the expectation maximization. 2 , α, and its update formula is as follows:

[0109]

[0110] Where Tr(·) represents the trace of the matrix;

[0111] (b3.4) Repeat steps E and M until the model parameters converge, that is, the iterative change value of the parameters is less than the set threshold.

[0112] S3. Construct a mapping relationship between the NO2 concentration time series data of the 11 national control stations and the Jing'an Yan'an West Road Traffic Station, so as to use the data of the 11 national control stations to trace back the data of the traffic stations when the data of the traffic stations are missing. Select the initial time window [T-2, T+2], and train a long short-term memory network (Bi-LSTM) with a bidirectional structure under the initial network hyperparameters to construct a mapping relationship between the two. This mapping relationship takes into account the lag of the impact of pollution diffusion, and uses the concentration information of other stations around the unknown data station from the past to the future relative to the traceback time to infer the concentration value of the unknown data station at the traceback time.

[0113] S4. Optimize the hyperparameters and adopt differentiated optimization strategies. The adjusted network hyperparameters include the time window size, learning rate, regularization parameter, batch size, network layer, number of neurons, and loss function. This method uses mean square error as the loss function, uses the Adadelta optimizer to dynamically adjust the learning rate, and uses the L2 regularization parameter. The remaining hyperparameters include the time window size, batch size, network layer, number of neurons, and define its hyperparameter space as Using Gaussian Process ) is used as a proxy model for hyperparameter optimization. The specific steps are as follows:

[0114] (a4.1) Initialization phase, defining the hyperparameter space In which the initial hyperparameter point set A={a1,a2,…,a n}, based on the initial hyperparameter point set, the Bi-LSTM network model is trained, and the validation set loss y={y1,y2,…,y n}, define the Gaussian process f(a):

[0115]

[0116] Among them, μ(a) represents the mean function, that is, E[f(a)], k(a,a′) is the kernel function, which characterizes the correlation between a and a′ in the hyperparameter space. Here, the radial basis function (RBF) is used, and the initial sampling set is used. Train the Gaussian process model f(a);

[0117] (a4.2) Calculate the predicted distribution. For candidate point a, the predicted mean μ(a) and the predicted variance σ 2 The calculation formula of (a) is as follows:

[0118] μ(a)= k(a,A)(K(A,A)+σ 2 I) -1 y (16)

[0119] σ 2 (a) = k(a,a) - k(a,A)(K(A,A) + σ 2 I) -1 k(A,a) (17)

[0120] Where A is the matrix of evaluated points, y is the corresponding objective function value vector, K(A,A) is the covariance matrix of the kernel function on the evaluated points, k(a,A) is the covariance vector between the candidate point and the evaluated points, and k(a,a) is the covariance between the candidate point and itself;

[0121] (a4.3) Select and optimize the acquisition function and use the expected improvement (EI) as the acquisition function. The calculation formula is as follows:

[0122]

[0123] Among them, a + is the current optimal point, ξ is a non-negative exploration parameter used to control the trade-off between exploration and exploitation of the algorithm, Φ(·) and φ(·) are the cumulative distribution function and probability of the standard normal distribution, respectively. Find the acquisition function α on the surrogate model. EI (a) The largest point a * , that is:

[0124]

[0125] (a4.4) Update the model, for the hyperparameter a * , train the model to get the validation set loss y * , add the new sampling point to the evaluated point set to form a set {(a1,y1),(a2,y2),…,(a n ,y n ),(a * ,y * )}, retrain the Gaussian process model using the updated dataset;

[0126] (a4.5) Check the stopping condition to see whether the stopping condition (maximum number of iterations, time limit or expected objective function value) is met. If so, stop. Otherwise, repeat (a4.2) to (a4.4).

[0127] Finally, the time window in this example is determined to be [T-6, T+6], the batch size is 32, the maximum number of iterations is set to 200, and the network parameters are shown in Table 3 below:

[0128] Table 3 Network parameters

[0129]

[0130]

[0131] S5. Using the PyTorch deep learning framework, the BPCA-LSTM model is constructed and compiled. The model is trained under the hyperparameter optimization framework to obtain the optimal hyperparameter settings and the corresponding network model under the preset time limit. The True vs Predicted Plot of the optimal model on the validation set is shown in the figure below. Figure 5 As shown in the figure, for the missing NO2 data of Jing'an West Yan'an Road Traffic Station in 2018, the model was applied with the NO2 pollutant concentration time series window data of 11 national control stations as input, and the NO2 hourly data series of Jing'an West Yan'an Road Traffic Station in 2018 was traced back. The results are as follows: Figure 6 shown.

[0132] The above description is only a description of the preferred embodiments of the present application, and is not intended to limit the scope of the present application. Any changes or modifications made by any person skilled in the art based on the above disclosed technical contents shall be deemed as equivalent effective embodiments and shall fall within the scope of protection of the technical solution of the present application.

Claims

1. A method for tracing back pollutant monitoring data at fixed sites of atmospheric environment monitoring, characterized in that: The following steps are involved: S1. Construct a BPCA-LSTM network model, including a Bayesian principal component analysis (BPCA) part and a long short-term memory (LSTM) network part; at the same time, according to the range of pollutant concentration monitoring data that needs to be traced back at M stations, select and form a usable complete data set from the historical data of pollutant concentration monitoring at all N stations; S2. Perform data cleaning on the selected data set, including pollutant concentration measurement standardization, outlier removal and missing point filling, and apply the Bayesian principal component analysis part to the cleaned data for data dimensionality reduction; S3. Set the initial data sliding window and BPCA-LSTM network model hyperparameters, and use the long short-term memory network part to construct the time series pollutant concentration data mapping relationship between the NM stations and each of the M stations; S4. Optimize and adjust the sliding window period length and the hyperparameters of the BPCA-LSTM network model; S5. Taking the historical data of NM stations corresponding to the period when the data of M stations were missing as input, the missing pollution concentration monitoring data of M stations were backtracked through the constructed BPCA-LSTM network model.

2. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 1, characterized in that: The data selected from the historical data of all N stations in step S1 are generally unified as the pollutant concentration monitoring data of N stations adjacent to the pollutant concentration data missing interval of the M stations that need to trace back data in the next 2 to 3 years. The relationship pattern of the concentration time series data between the N stations covered by the selected data effectively represents the concentration time series data mapping relationship between the NM stations and the M stations in the data missing interval of the M stations; Data integrity is defined by the following formula: Among them, Ω represents the completeness of the data, Q represents the total amount of data that should be monitored in the data set, and Q a Indicates the total amount of data in the dataset that has not been monitored.

3. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 1, characterized in that: In step S2, data cleaning includes pollutant concentration measurement standardization, outlier removal and missing point filling, as follows: (a3.1) For pollutant concentration measurement standardization, if the stored data uses the concentration value under standard conditions (273.15K, 101.325kPa), it is converted to the concentration value under actual conditions. The conversion formula is as follows: Among them, C r Indicates the concentration value under actual conditions, C s Indicates the concentration value under standard conditions, T r Indicates the ambient temperature under actual conditions, which uses the Kelvin temperature scale; (a3.2) For outlier removal, that is, to remove abnormal detection records and data outliers, outliers are defined as values ​​less than the lower limit or greater than the upper limit. The upper and lower limit calculation formulas are: Among them, u and L represent the upper and lower limits respectively, Q3 and Q1 represent the upper quartile and lower quartile respectively, and IQR represents the interquartile range, that is, the distance between Q3 and Q1; (a3.3) For missing point filling, the missing points of the original data are filled with the average concentration of pollutants within the selected data range. For the data removed due to outliers, if they were originally abnormal detection records, they are filled with the above-mentioned average concentration value. If they are outliers, those less than the lower limit are filled with the minimum value of the data after removing the outliers, and those greater than the upper limit are filled with the maximum value of the data after removing the outliers.

4. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 1, characterized in that: In step S2, Bayesian principal component analysis is used to reduce the dimension of the cleaned data. The dimension reduction process is completed using the expectation-maximization algorithm (EM algorithm). For the actual concentration vector at a certain moment D is the feature dimension, where x is the vector of actual concentrations of all stations except the station with missing data at that moment. The goal of dimensionality reduction is to find a low-dimensional representation and the corresponding weight matrix and a bias vector So that: x=Wz+μ+ε (4) Among them, ε is the noise term; there are the following prior distribution assumptions about x, z, and W, where represents Gaussian distribution, I is the identity matrix, σ 2 is the observation noise variance: The parameter set that needs to be estimated under the current assumption is θ = {W, μ, σ 2 ,α}, where μ is estimated using the following formula: Parameters W, σ 2 , α is iteratively updated through the EM algorithm, and the EM algorithm estimates the model parameters by alternately updating the E step, i.e., the expectation step, and the M step, i.e., the maximization step.

5. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 4, characterized in that: The EM algorithm: Specific steps: (b3.1) Initialize parameters W, σ 2 , α; (b3.2) In step E, the expectation and covariance of the posterior distribution of the latent variable z are calculated as follows: where M = W T W + σ 2 I (10) (b3.3) In the M step, the parameters W and σ are updated according to the expectation maximization. 2 , α, and its update formula is as follows: Where Tr(·) represents the trace of the matrix; (b3.4) Repeat steps E and M until the model parameters converge, that is, the iterative change value of the parameters is less than the set threshold.

6. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 1, characterized in that: In step S3: A mapping relationship between the time series pollutant concentration data of the NM stations and each of the M stations was constructed to trace back the data of the M station using the data of the remaining NM stations when the data of the M station is missing.

7. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 6, characterized in that: In step S3: Specifically, For the current moment T, a time window [Tn, T+n] is selected, where n is a hyperparameter for adjusting the window size. The pollutant concentration data within the time window of NM stations is used as input, and the concentration value of the backtracked station required at the current moment T is used as output. The mapping relationship between the two is constructed using a long short-term memory network (LSTM) under the selected network hyperparameters. The concentration information of a time window from the past to the future relative to the backtracking moment of other stations around the unknown data station is used to infer the concentration value of the unknown data station at the backtracking moment.

8. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 7, characterized in that: In step S3: Long short-term memory network (LSTM), this network model adopts a bidirectional structure, that is, a bidirectional LSTM (Bi-LSTM) is constructed to capture and memorize the concentration relationship patterns of time series data in the time window in both the past to future and future to past directions.

9. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 1, characterized in that: Step S4 adopts a differentiated hyperparameter optimization strategy. The adjusted network hyperparameters include the time window size, learning rate, regularization parameter, batch size, network layer, number of neurons, and loss function. In this method, mean square error is used as the loss function, the Adadelta optimizer is used to dynamically adjust the learning rate, and the L2 regularization parameter is used. The remaining hyperparameters include the time window size, batch size, network layer, and number of neurons. The hyperparameter space is defined as Using Gaussian Process ) is used as a proxy model for hyperparameter optimization. The specific steps are as follows: (a4.1) Initialization phase, defining the hyperparameter space In which the initial hyperparameter point set A={a1,a2,…,a n }, based on the initial hyperparameter point set, the Bi-LSTM network model is trained, and the validation set loss y={y1,y2,…,y n }, define the Gaussian process f(a): Among them, μ(a) represents the mean function, that is, k(a,a′) is a kernel function that characterizes the correlation between a and a′ in the hyperparameter space. Here, the radial basis function (RBF) is used, and the initial sampling set is used. Train the Gaussian process model f(a); (a4.2) Calculate the predicted distribution. For candidate point a, the predicted mean μ(a) and the predicted variance σ 2 The calculation formula of (a) is as follows: μ(a)=k(a,A)(K(A,A)+σ 2 I) -1 y (16)p 2 (a)=k(a,a)-k(a,A)(K(A,A)+σ 2 I) -1 k(A,a)(17) Where A is the matrix of evaluated points, y is the corresponding objective function value vector, K(A,A) is the covariance matrix of the kernel function on the evaluated points, k(a,A) is the covariance vector between the candidate point and the evaluated points, and k(a,a) is the covariance between the candidate point and itself; (a4.3) Select and optimize the acquisition function and use the expected improvement (EI) as the acquisition function. The calculation formula is as follows: Among them, a + is the current optimal point, ξ is a non-negative exploration parameter used to control the trade-off between exploration and exploitation of the algorithm, Φ(·) and φ(·) are the cumulative distribution function and probability of the standard normal distribution, respectively. Find the acquisition function α on the surrogate model. EI (a) The largest point a * , that is: (a4.4) Update the model, for the hyperparameter a * , train the model to get the validation set loss y * , add the new sampling point to the evaluated point set to form a set {(a1,y1),(a2,y2),…,(a n ,y n ),(a * ,y * )}, retrain the Gaussian process model using the updated dataset; (a4.5) Check the stopping condition and see if it is met. If so, stop. If not, repeat (a4.2) to (a4.4).

10. A method for tracing back pollutant monitoring data at a fixed site for atmospheric environment monitoring as claimed in claim 2, characterized in that: In step S1, the selected data: The data of NM stations with data completeness less than 90% were discarded.

Images