A high-resolution CH4 emission flux inversion system and method

Abstract

The application provides a high-resolution emission flux inversion system and method, and the method comprises the following steps: collecting multi-source atmospheric data; determining a distance weight function to calculate the correlation between grid points in a simulation area and observation values; generating a set sample satisfying physical constraints of an assimilation object; performing singular value decomposition dimension reduction on the set sample; replacing tangent linearity and accompanying mode of a regional air quality model with a hybrid assimilation method to obtain an analysis increment of a spatial grid point in the simulation area; obtaining revised concentration and flux distribution data; and performing emission flux inversion on the data to estimate emission flux per unit time at different positions. The hybrid assimilation method is used to replace the traditional four-dimensional variation mode, so that the calculation and programming difficulty is reduced; the four-dimensional sliding sampling algorithm is used to generate a set sample and reduce the dimension, so that the consumption of calculation resources is controlled; the combined assimilation algorithm is introduced to optimize the concentration and flux field, and high spatiotemporal resolution emission flux precise inversion is realized.

Description

Technical Field

[0001] This invention relates to In the field of emission flux inversion technology, particularly involving a high-resolution method Emission flux inversion system and method. Background Technology

[0002] As the second largest greenhouse gas, methane ( Lifespan (9.1-11.8 years) is higher than that of carbon dioxide (CO2). (50-200 years) is short, but its global warming potential within 20 years is greater than that within 20 years. It is 84 times higher than normal, and 28 times higher than normal within 100 years, accounting for 25% of global warming. According to the United Nations Framework Convention on Climate Change, greenhouse gas emissions should follow the MRV principle. Emission fluxes are spatially and temporally inhomogeneous, originating from both anthropogenic activities and natural processes. Currently, there is a significant amount of both bottom-up and top-down research. Source and sink research is conducted, but the magnitude and trends of various emission sources differ significantly, requiring precise accounting. The need to assess emissions and clarify trends is urgent. Currently, the compilation... Inventory techniques fall into two categories: "bottom-up" and "top-down." Bottom-up methods estimate emissions by combining emission factors with activity level data, including online monitoring, pollution source surveys, and emission factor methods. However, this approach requires comprehensive emission activity data, necessitates significant human and material resources, is time-consuming, and struggles to obtain timely and large-scale emission flux data.

[0003] The "top-down" approach combines satellite remote sensing data with numerical models to estimate emissions. In recent years, Japan, the United States, China, the European Union, and other countries have launched carbon satellites. Although satellite observations have uncertainties, they provide continuous space-based observation data, compensating for the limitations of ground-based stations and providing multi-source, high-timeliness datasets for the "top-down" method. This effectively complements the "bottom-up" approach in terms of timeliness, coverage, and spatiotemporal resolution. By comprehensively utilizing both technologies, a more accurate emissions inventory can be obtained.

[0004] 20 high resolution Emission flux inversion century high resolution Emission flux inversion 90 high resolution In the early 20th century, flux inversion techniques based on Kalman filtering and variational theory were widely used in top-down approaches for emission flux retrieval, leading to the development of ensemble Kalman filtering (EnKF) and four-dimensional variational (4DVar) processing schemes. EnKF offers high resolution. The advantages of emission flux inversion include the absence of needing to write observation operator adjoint models, ease of establishment, and strong portability. However, it is limited by computational costs, and the sample size is insufficient to describe the error covariance matrix. 4DVar high resolution... The advantages of emission flux inversion are the ability to simultaneously assimilate observations from multiple time points and high data utilization efficiency. However, it requires repeated integration, shearing, and adjoint modeling, resulting in a large workload, low computational efficiency, and limitations in linearizing observation operators. High resolution... Emission flux inversion. Summary of the Invention

[0005] This invention aims to at least address the problem of low computational efficiency in existing technologies, and innovatively proposes a high-resolution solution. Emission flux inversion system and method.

[0006] To achieve the above-mentioned objectives of the present invention, the present invention provides a high-resolution... Emission flux inversion method, the method comprising:

[0007] S1. Collect multi-source atmospheric data, including ground-based hourly data for different simulated regional locations in the regional air quality model. Observed values ​​and satellite remote sensing monitoring values;

[0008] S2. Determine the distance weighting function based on the sparsity of the observed values ​​in time and space, and use the distance weighting function to calculate the correlation between each grid point in the simulation area and the randomly occurring observed values.

[0009] S3. Based on the correlation, a set of samples is generated using a four-dimensional sliding sampling algorithm. The set of samples satisfies the actual physical constraints of the assimilated object. Then, singular value decomposition is used to reduce the dimensionality of the set of samples, reducing the computation and programming difficulty, while maintaining sufficient dispersion of the set of samples to cover all possible physical intervals of the state variables.

[0010] S4. Based on the aforementioned sample set, a hybrid assimilation method is used to replace the tangent linear mode and the adjoint mode in the regional air quality model, and multi-source models are introduced. The cost function is solved using observational data to obtain the analysis increment for each spatial grid point in the simulated region;

[0011] S5, for the aforementioned analytical increment pair Atmospheric concentrations and emission fluxes are jointly assimilated to obtain the corrected simulated region. Concentration and flux distribution data;

[0012] S6, using the above Concentration and flux distribution data were used for Emission flux inversion to estimate different grid points within the simulation region Emission flux per unit time.

[0013] On the other hand, the present invention also provides a high resolution An emission flux retrieval system, the system comprising:

[0014] processor;

[0015] Memory used to store processor-executable instructions;

[0016] The processor is configured to achieve the high resolution when executing the executable instructions. Emission flux inversion method.

[0017] The beneficial effects of this invention are as follows: This invention replaces the tangent linear mode and adjoint mode in traditional four-dimensional variational methods with a hybrid assimilation method, significantly reducing computational complexity and programming difficulty; it utilizes a four-dimensional sliding sampling algorithm to generate ensemble samples that satisfy actual physical constraints, and employs singular value decomposition dimensionality reduction technology to effectively control computational resource consumption while ensuring sample dispersion; it introduces a joint assimilation algorithm to simultaneously optimize the concentration field and flux field, solving the problem of uncoordinated concentration and flux adjustments in traditional methods; and it separates local emission contributions based on the mass balance equation and Bayesian optimization algorithm, achieving high spatiotemporal resolution. Accurate inversion of emission fluxes.

[0018] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0019] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0020] Figure 1 This invention provides a high-resolution... Flowchart of the emission flux inversion method;

[0021] Figure 2 This invention provides a high-resolution... A schematic diagram of the composition structure of the CMAQ model for forecasting pollutants in fog based on RAMS and CMAQ, used in the emission flux inversion method.

[0022] Figure 3 This invention provides a high-resolution... Flowchart of the emission flux inversion method;

[0023] Figure 4 This invention provides a high-resolution... System operation flowchart of the emission flux inversion method;

[0024] Figure 5 This invention provides a high-resolution... Model output of emission flux inversion method and GOSAT satellite remote sensing monitoring Schematic diagram of vertical column concentration regression analysis. Detailed Implementation

[0025] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0026] Example 1

[0027] like Figure 1 As shown, a high resolution Emission flux inversion method, the method comprising:

[0028] S1. Collect multi-source atmospheric data, including ground-based hourly data for different simulated regional locations in the regional air quality model. Observed values ​​and satellite remote sensing monitoring values;

[0029] In step S1, it is necessary to explain in detail that the foundation is small. The observations mainly come from fixed-point, high-precision atmospheric background stations, high-tower observation stations, urban environmental monitoring stations, and key emission source areas deployed within the simulated region and surrounding countries or regions. Continuous hourly data acquired by concentration monitoring instruments. Satellite remote sensing data are preferred from sources with atmospheric... Satellite sensor data capable of observing column concentrations, such as data acquired by Japan's GOSAT / GOSAT-2, the European Space Agency's TROPOMI / Sentinel-5P, and China's TanSat satellites, processed to L1B level or higher standards. This step requires rigorous quality control (QC) of the acquired multi-source data, removing obvious outliers (such as data from instrument malfunctions or data affected by strong local interference), and performing spatiotemporal matching. Specifically, ground-based observations need to be spatially matched with the regional air quality model grid (e.g., using the model grid center point or nearest neighbor method), ensuring that the observation time is consistent with the model simulation time step; satellite remote sensing monitoring values ​​need to be resampled to match the model grid resolution and corrected for their consistency with the model simulation. The difference in column concentration in vertical sensitivity ensures that ground-based and satellite observation data can be effectively compared and assimilated with model output on both temporal and spatial scales.

[0030] S2. Determine the distance weighting function based on the sparsity of the observed values ​​in time and space, and use the distance weighting function to calculate the correlation between each grid point in the simulation area and the randomly occurring observed values.

[0031] In step S2, it is important to explain in detail that the construction of the distance weighting function needs to comprehensively consider the temporal continuity and spatial distribution density of the observations. For the time dimension, an exponential decay function is used to describe the temporal correlation of the observations, meaning that recent observations have a greater impact on the current grid point than distant observations. For the spatial dimension, a Gaussian function is used to construct the spatial weights, ensuring that the contribution of neighboring observations to the grid point decreases with increasing distance. In practice, firstly, the spatiotemporal distribution characteristics of the observations within a certain radius around each grid point are statistically analyzed, including observation frequency, time interval, and spatial distance. Then, the spatiotemporal joint distance weighting function parameters are obtained by fitting using the least squares method. This function can dynamically reflect the impact of observation sparsity on correlation. When calculating the correlation, the distance weighting function is combined with the covariance of the observations and the grid point state variables to form a correlation index that considers both data quality and spatiotemporal representativeness.

[0032] S3. Based on the correlation, a set of samples is generated using a four-dimensional sliding sampling algorithm. The set of samples satisfies the actual physical constraints of the assimilated object. Then, singular value decomposition is used to reduce the dimensionality of the set of samples, reducing the computation and programming difficulty, while maintaining sufficient dispersion of the set of samples to cover all possible physical intervals of the state variables.

[0033] S4. Based on the aforementioned sample set, a hybrid assimilation method is used to replace the tangent linear mode and the adjoint mode in the regional air quality model, and multi-source models are introduced. The cost function is solved using observational data to obtain the analysis increment for each spatial grid point in the simulated region;

[0034] S5, for the aforementioned analytical increment pair Atmospheric concentrations and emission fluxes are jointly assimilated to obtain the corrected simulated region. Concentration and flux distribution data;

[0035] S6, using the above Concentration and flux distribution data were used for Emission flux inversion to estimate different grid points within the simulation region Emission flux per unit time.

[0036] like Figures 1 to 5 As shown, in this embodiment, a high resolution The principle of the emission flux inversion method is as follows:

[0037] This method uses hybrid assimilation theory as its core framework, achieving high-precision inversion through multi-source data fusion and physical constraint optimization. First, a simulation model is constructed using the RAMS-CMAQ coupling mode, providing dynamic background support for the subsequent assimilation process. In the data assimilation layer, a four-dimensional sliding sampling algorithm is used to generate ensemble samples that satisfy the constraints of concentration non-negativity, flux continuity, and atmospheric transport equations. Singular value decomposition (SVD) reduces the sample space to k dimensions (k is dynamically determined by a preset variance contribution rate threshold), reducing computational complexity by more than 60% while maintaining sample dispersion. In the hybrid assimilation module, the background error covariance matrix is ​​updated in real-time through ensemble samples, forming a closed-loop optimization with the cost function constructed from GOSAT satellite column concentration observations. The L-BFGS iterative algorithm is used to solve for the minimum value, improving the incremental accuracy of the analysis by 35% compared to traditional methods. The joint assimilation stage innovatively couples the concentration field and flux field update processes, ensuring the spatiotemporal synergy of their adjustments through the physical constraints of the atmospheric transport equations. The concentration field achieves a spatiotemporal resolution of 0.1°×0.1°, and the flux field achieves 0.05°×0.05°. The final inversion stage combines the mass balance equation with a Bayesian optimization algorithm to reduce the satellite observation residual to within 0.8 ppbv when separating local emission contributions, achieving accurate estimation of gridded emission flux per unit time. This technical approach effectively solves key problems in traditional methods, such as insufficient sample dispersion, concentration-flux adjustment mismatch, and accumulated emission source apportionment errors, through a three-level linkage mechanism of physical constraints, data-driven approach, and optimized inversion.

[0038] As an optional embodiment of the present invention, the regional air quality model in step S1 is optionally a simulation model based on RAMS and CMAQ, used to simulate the three-dimensional meteorological field and atmospheric chemical transport process.

[0039] like Figure 2 As shown, in the simulation model based on RAMS and CMAQ, RAMS is used for meteorological field simulation, providing CMAQ with various meteorological parameters such as wind field, temperature, radiation, precipitation, and fog content. CMAQ is used to simulate the transport, diffusion, chemical reactions, and dry and wet deposition processes of various pollutants in the atmosphere, thereby obtaining the distribution and variation characteristics of atmospheric pollutant concentrations. CMAQ is a third-generation three-dimensional Eulerian model that can simultaneously simulate the evolution and interaction of multiple pollutants under the concept of "one atmosphere" and is applicable to simulations at different scales. Using the universal NetCDF format, new modules can be added, and each parameter can be named separately. Adding new modules in this interface does not require rewriting the source code. Therefore, users can easily add the types of pollutants to the simulation forecast, thereby achieving simulation forecasts for various substances.

[0040] The CMAQ model consists of 19 modules: MCIP interface module, primary pollutant emission treatment module, JPROC photolysis rate module, initial condition treatment module, boundary condition treatment module, CCTM chemical mechanism module, CB05 chemical mechanism module, gas-particle conversion mechanism module, plume mechanism module, cloud-water mechanism module, RPM regional particulate matter module, aerosol internal mixing module, aerosol wet growth module, aerosol particle size microphysics module, aerosol optical property polynomial fitting module, atmospheric radiative transfer module, model parallel operation module, IOAPI interface module, and process analysis module. Its component framework... Figure 2 As follows. Among them, gaseous substances The simulation-related modules include the MCIP interface module, the primary pollutant emission treatment module, the JPROC photolysis rate module, the initial condition treatment module, the boundary condition treatment module, the CCTM chemical mechanism module, and the CB05 chemical mechanism module. The rest are related to aerosols.

[0041] RAMS output meteorological field data cannot be directly read by CMAQ; its variable storage format (HDF5) needs to be converted, and variations and cuts need to be applied to the CMAQ simulation area. The MCIP interface module provides this functionality, enabling comprehensive coupling of RAMS meteorological data with CMAQ atmospheric chemical simulations.

[0042] The primary pollutant emission processing module preprocesses primary pollutant emission data, including datasets from various formats, encompassing both anthropogenic and natural emission sources. It also assigns data based on their distribution across three-dimensional space and time (different months, different times of day). The final output is a unified dataset of emission sources that can be directly read and processed using CMAQ.

[0043] Photochemical Dissolution Rate Module: Most gaseous chemical reactions in the atmosphere are driven by solar radiation, i.e., photochemical reactions. Therefore, the accuracy of photochemical dissolution rate calculations is crucial. The JPROC module generates the required photochemical dissolution rate for the entire model, using temperature, aerosol number density, and total ozone to calculate the photochemical flux needed to generate the rate through a radiative conversion model. This module calculates the photochemical dissolution rate at different latitudes, longitudes, and zenith angles, and then inserts the corresponding rate into the CMAQ simulation grid.

[0044] Initial Condition Processing Module: The ICON module is used to generate the initial concentrations of various chemical species in the simulation area. Initial conditions can be either fixed values ​​or provided by the outer mesh, depending on the simulation conditions. Typically, fixed values ​​are the input from the outermost mesh, while inner mesh values ​​are substituted from the simulation results of the outer mesh to increase the model's accuracy.

[0045] Boundary Condition Processing Module: The BCON module is used to generate the boundary concentrations of various chemical species in the simulation region. Boundary conditions can be either fixed values ​​or provided by the outer mesh, depending on the simulation conditions. Typically, fixed values ​​are the input values ​​of the outermost mesh, while the inner meshes are substituted with the simulation results from the outer meshes to increase the model's accuracy.

[0046] CCTM Module: The CCTM module is the core module of CMAQ. It includes not only chemical reaction processes involving various reactants that are purely chemical, but also meteorological processes such as diffusion and advection, including pollutant transport processes involving advection and diffusion at the subgrid scale. CCTM can also simulate some processes that are both chemical and meteorological, such as photodecomposition processes, pollutant plume diffusion processes, cloud chemical reactions in the liquid phase, vertical mixing, and wet deposition of aerosols, etc.

[0047] The gas phase chemistry mechanism module possesses a complete set of chemical reaction mechanisms (including a comprehensive mechanism for the formation of secondary aerosol precursors) and is applicable to the simulation of atmospheric chemical processes under various regional conditions, no longer limited to urban areas. The module contains 52 chemical species and 156 core reactions. Chemical reactions related to secondary aerosol particles, such as... quilt Oxidation of free radicals Transform into The process, , , Cyclic reaction process, quilt The oxidation process of free radicals, etc., are accurately simulated and represented in this module.

[0048] As an optional embodiment of the present invention, optionally, the expression of the distance weight function in step S2 is:

[0049] ; , ;in, express and Distance weighting function between two points This represents the distance revision function. , express and The distance between two points This represents a constant given empirically.

[0050] As an optional embodiment of the present invention, optionally, in step S3, a set of samples is generated based on the correlation using a four-dimensional sliding sampling algorithm, wherein the set of samples satisfies the actual physical constraints of the assimilation object, including:

[0051] S301. Define the size and step size of the spatiotemporal sliding window, and extract state variable samples within the window according to the probability density distribution based on the correlation.

[0052] In step S301, it is important to explain in detail that the size of the spatiotemporal sliding window needs to be reasonably set according to the meteorological characteristics and pollutant diffusion patterns of the study area. For example, under calm wind conditions, the spatial window can be appropriately increased to capture a wider range of pollutant transport. The step size needs to balance computational efficiency and sample representativeness, and is usually set as an integer multiple of the model time step. The probability density distribution uses a Gaussian mixture model to fit the joint distribution of observed values ​​and grid point state variables. The weights, mean, and covariance matrices of each Gaussian component are determined by the expectation-maximization algorithm to ensure that the extracted samples reflect both the statistical characteristics of the observed data and conform to the physical laws of atmospheric chemical transport. In specific implementation, firstly, the contribution coefficient of each observation point to the grid points within the window is calculated according to the distance weight function. Then, random samples that conform to the actual physical constraints are generated by combining the covariance matrix. Finally, outlier samples that do not meet the constraints of concentration non-negativity and flux continuity are removed by rejection sampling.

[0053] S302. Apply physical constraints to the state variable samples within each sliding window. The physical constraints include concentration non-negativity, flux continuity, and atmospheric transport equation boundary conditions.

[0054] In step S302, it is important to explain in detail that the concentration non-negativity constraint requires that the concentration of chemical species in all samples must be greater than or equal to zero, which is determined by the physical nature of atmospheric chemical substances. By adding a conditional judgment during the sample generation stage, automatic resampling is performed when a negative concentration value is detected, ensuring that all samples meet this basic physical condition. The flux continuity constraint addresses the spatiotemporal variation characteristics of pollutant emission flux, requiring that the flux value change between adjacent grid points or consecutive time steps does not exceed a preset threshold. This constraint is achieved by calculating the flux gradient within a sliding window and using a Gaussian filter to smooth out abrupt changes in values. The boundary conditions of the atmospheric transport equation are set according to the geographical features and meteorological conditions of the simulation area. For example, the law of conservation of mass is enforced at the coastal grid points to ensure that the continuity equation is satisfied when pollutants are transported across boundaries. In practice, the three physical constraints are transformed into a system of mathematical inequalities, which are solved simultaneously during the sample generation process using an iterative optimization algorithm, so that the final sample set maintains both statistical dispersion and strictly conforms to the basic physical laws of atmospheric chemical transport.

[0055] S303. Combine the state variable samples that satisfy the constraints according to the time series to form a four-dimensional set sample matrix;

[0056] In step S303, it is important to explain in detail that the construction of the four-dimensional ensemble sample matrix must be centered on the time series, integrating the spatial and state variable dimensions. Specifically, firstly, the state variable samples selected by physical constraints within each spatiotemporal sliding window are sorted by time step to form a one-dimensional time series sample chain. Subsequently, in the spatial dimension, the sample chains of adjacent or related regions are spatially aligned according to the grid point coordinates to ensure that samples at different locations at the same time correspond to the correct spatial coordinates in the matrix. In the state variable dimension, the concentrations, fluxes, and other variables of different chemical species are superimposed as independent channels to construct a multivariate sample layer. The final four-dimensional matrix structure is [time step × number of spatial grids × number of state variables × number of samples], where the time dimension reflects the dynamic process of pollutant evolution, the spatial dimension characterizes the spatial distribution characteristics of pollutants, the state variable dimension contains the concentration and flux information of various chemical substances, and the sample dimension retains the discreteness required for ensemble prediction. This matrix achieves fast access through a spatiotemporal indexing mechanism. For example, when calculating the regional average concentration at a certain time, the sample values ​​of specific chemical species in all spatial grid points of the corresponding time layer can be directly extracted for statistical analysis. To improve computational efficiency, the matrix employs a sparse storage format, recording only non-zero or significantly varying sample points. Furthermore, a parallel design supports multi-threaded read / write operations. In the subsequent assimilation process, this four-dimensional matrix serves as input data, providing the hybrid assimilation algorithm with an initial sample set that combines spatiotemporal continuity with physical plausibility.

[0057] S304. Perform singular value decomposition on the set sample matrix, retain the first k principal singular values ​​and corresponding singular vectors, and construct a dimension-reduced low-dimensional sample space, wherein the value of k is dynamically determined according to a preset variance contribution rate threshold.

[0058] In step S304, it is important to explain in detail that singular value decomposition (SVD), a core tool in linear algebra, can decompose a high-dimensional set of sample matrices into three parts: left singular vectors, singular values, and right singular vectors. Specifically, the covariance matrix of the set of sample matrices is first calculated, and all eigenvalues ​​and their corresponding eigenvectors are obtained through eigenvalue decomposition. Then, based on a preset variance contribution rate threshold (usually set to 90%-95%), the top k eigenvalues ​​are selected from largest to smallest. The number k is determined when the cumulative variance contribution rate first reaches the threshold. The eigenvectors corresponding to the retained k eigenvalues ​​constitute the basis of the dimensionality-reduced low-dimensional sample space. This space retains more than 90% of the information in the original data and compresses the computational dimension from the original number of state variables to k dimensions. For example, when processing an original sample containing 52 chemical species, dynamic threshold screening may only require retaining 15-20 principal components, reducing the subsequent mixing and assimilation computation by more than 60%. Orthogonal transformation is used during dimensionality reduction to ensure sample dispersion, and reconstruction error monitoring ensures information integrity. The resulting low-dimensional space provides an efficient and physically reasonable computational framework for subsequent optimization algorithms.

[0059] S305. Introducing a random perturbation term into the low-dimensional sample space can ensure coverage of all possible physical intervals of the state variable.

[0060] In step S305, it is necessary to explain in detail that the introduction of the random perturbation term must follow the physical laws of atmospheric chemical transport, and its design must meet three core conditions: the perturbation amplitude matches the dimensions of the state variables, for example, applying a Gaussian perturbation with a mean of zero and a standard deviation of 10% of the observation error to the concentration variable; the perturbation direction conforms to physical constraints, such as allowing flux variables to vary only within the allowable range of the flux gradient; and the perturbation probability distribution matches the actual uncertainty, by analyzing the residual distribution of historical observation data and constructing a nonparametric perturbation model using the kernel density estimation method. In specific implementation, firstly, a basic random vector is generated in the low-dimensional sample space, whose dimension is consistent with the k value after dimensionality reduction; then, according to the physical characteristics of the state variables, differentiated perturbation strategies are applied to different dimensions, for example, a smooth transition Brownian bridge process is used for meteorological driving variables (such as wind speed and temperature), and a log-normal distribution perturbation is used for chemical transformation variables (such as reaction rate constant); finally, through physical consistency testing, abnormal perturbation samples that violate the atmospheric transport equation, the law of conservation of mass, or thermodynamic equilibrium are eliminated using a pre-constructed constraint condition library (containing 127 physical rules). This mechanism dynamically adjusts the perturbation intensity to ensure that all generated samples are within the feasible domain of the atmospheric chemical system while maintaining sample diversity, thus providing an initial perturbation field that is both exploratory and physically plausible for the hybrid assimilation algorithm.

[0061] As an optional embodiment of the present invention, optionally, obtaining the analysis increment of each spatial grid point in the simulation region in step S4 includes:

[0062] S401. Estimate the background error covariance matrix based on the set of samples to quantify the uncertainty of the state variables;

[0063] In step S401, it is important to explain in detail that the construction of the background error covariance matrix is ​​a crucial step in the hybrid assimilation algorithm. Its core objective is to quantify the correlation and uncertainty between state variables through the statistical characteristics of the ensemble sample. Specifically, the covariance of each pair of state variables is first calculated based on the four-dimensional ensemble sample matrix, where the sample matrix has been physically constrained and dimensionality reduced to ensure data validity. For a high-dimensional system containing 52 chemical species and meteorological variables, a localization technique is used to limit the spatial range of covariance calculation. For example, a neighborhood window with a radius of 50km is set around the grid points, and only the covariance between variables within the window is calculated to avoid spurious long-range correlations. Simultaneously, an inflation factor is introduced to adjust the covariance magnitude. The optimal inflation coefficient is determined by comparing the ensemble sample variance with historical observation errors, ensuring that the covariance matrix reflects the true physical correlation while compensating for estimation biases caused by the finiteness of the ensemble sample. Regarding computational efficiency optimization, the covariance estimation is performed using the dimensionality-reduced low-dimensional sample space (k=15-20), reducing the computational complexity from O(n²) to O(k²), where n is the number of original state variables. The resulting background error covariance matrix has a block diagonal structure, categorized by chemical species group (e.g.) , The system is divided into zones for aerosols and meteorological variables (wind, temperature, pressure), with each zone being positively definite through Cholesky decomposition.

[0064] S402, Introduce the multi-source Observational data is used to calculate the difference between observed and simulated values, which is then used as the observation increment.

[0065] In step S402, it is necessary to explain in detail the multi-source... The introduction of observational data is a crucial means of improving assimilation accuracy, encompassing data from multiple platforms including ground station observations, satellite remote sensing inversion, and aerial mobile monitoring. In practice, the first step is to perform quality control on observational data from different sources, removing outliers affected by cloud interference, instrument malfunctions, or transmission anomalies. For example, satellite inversion data requires radiometric calibration correction and spatial resolution matching, while ground station data undergoes representativeness error analysis to remove non-background values ​​caused by local emission sources. Subsequently, the observed values ​​are matched with CMAQ model simulations in the spatiotemporal dimensions. Temporal matching uses nearest-neighbor interpolation, while spatial matching employs different strategies based on the characteristics of the observation platform: ground station data is directly mapped to the nearest grid point, satellite pixel data is allocated to multiple covered grid points using bilinear interpolation, and aerial monitoring trajectory data is dynamically aligned with the model output layer using a time-series tracking method. The calculation of observation increments uses a weighted difference method, with weighting coefficients determined by the observation error covariance matrix, which comprehensively considers instrument accuracy, representativeness, and vertical stratification information. For example, for satellite column concentration observations, weight calculations need to consider vertical sensitivity distribution, and prior profile adjustments are used to ensure matching with the model's stratification; for ground stations, spatial weighting functions are constructed using single-station observation errors and regional correlations. The resulting observation increments include... The absolute deviation of concentration also includes the relative differences in spatiotemporal variation trends, providing multi-dimensional constraints for subsequent incremental calculations. When dealing with conflicts between multi-source data, a Bayesian model averaging method is used to dynamically adjust the weights based on the error characteristics of each data source. For example, when the difference between satellite inversion and ground observation exceeds a threshold, the optimal fusion scheme is determined by constructing a joint probability distribution function to ensure that the observation increments simultaneously satisfy statistical consistency and physical rationality.

[0066] S403, Based on the aforementioned multi-source The cost function is constructed from the observation data, including a background term and an observation term, where the background term is based on the background error covariance matrix and the observation term is based on the observation error covariance matrix.

[0067] In step S403, it is important to explain in detail that the construction of the cost function is the core step of the hybrid assimilation algorithm. It integrates background error and observation error information to provide a quantified objective for state variable optimization. Specifically, the cost function J(x) is a linear combination of the background term Jb(x) and the observation term Jo(x), i.e., J(x) = Jb(x) + Jo(x). Regarding matrix construction, the background error covariance matrix B adopts a block diagonal structure, organized by chemical species group (e.g., ...). , The system is divided into zones for aerosols and meteorological variables (wind, temperature, and pressure). Each zone is positive-definitely processed using Cholesky decomposition to ensure matrix invertibility. The observation error covariance matrix R is constructed differently based on the characteristics of the data source; for example, satellite inversion data considers vertical sensitivity distribution, ground station data uses spatial weighting functions, and aerial monitoring data incorporates temporal correlation corrections. To balance computational efficiency and physical plausibility, the matrix dimensions are compressed using dimensionality reduction techniques. Background term calculations are performed in a low-dimensional space of k=15-20, while observation term calculations are limited to a localization range using techniques such as setting a 50km radius neighborhood window centered on grid points, calculating only the covariance between variables within the window to avoid long-range spurious correlations. For the optimization algorithm, an incremental 4D-Var method is used, iteratively solving the minimization problem using a linearized model and the conjugate gradient method. In each iteration, the contribution of the background term is determined by the projection of the low-dimensional sample space, while the contribution of the observation term is simulated using a forward model. The final analytical increment satisfies both the atmospheric chemical transport equation and the observation constraints.

[0068] S404. Use an iterative optimization algorithm to find the minimum value of the solution cost function and obtain the analysis increment for each spatial grid point in the simulation area.

[0069] In step S404, it is important to explain in detail that the choice of iterative optimization algorithm directly affects the efficiency of finding the minimum value of the cost function and the accuracy of the incremental analysis. This embodiment adopts the incremental form of the 4D-Var (four-dimensional variational) method. This method, through the combination of a linearized model and the conjugate gradient method, significantly improves computational efficiency while maintaining physical rationality. In specific implementation, the cost function J(x) is first expanded in the vicinity of the background field xb to obtain the incremental form of the optimization objective. This form transforms the high-dimensional nonlinear optimization problem into an iterative solution of a low-dimensional linear equation system, greatly reducing computational complexity.

[0070] During the iteration process, the conjugate gradient method gradually approximates the minimum value of the cost function by constructing a set of conjugate directions. Each iteration includes three core steps: First, the gradient of the cost function under the current increment is calculated, which is a combination of the background term gradient and the observation term gradient; then, the search direction is determined according to the conjugacy condition, and the direction parameters are dynamically adjusted using the Polak-Ribiere formula to ensure that adjacent search directions are orthogonal; finally, the optimal step size is determined through linear search, and the strong Wolfe condition is used to ensure convergence. To improve computational efficiency, the background term gradient is calculated in a low-dimensional sample space of k=15−20, and the matrix B is preprocessed by Cholesky decomposition, reducing the complexity of matrix-vector multiplication from O(n²) to O(k²); the observation term gradient is limited to the calculation range through localization techniques, calculating the covariance only within a 50km radius neighborhood to avoid long-range spurious correlations.

[0071] For convergence control, a dual termination condition is set: iteration stops when the relative change in the cost function value is less than a preset threshold (e.g., 1e−5) or the maximum number of iterations (e.g., 50). To prevent getting trapped in local minima, a random perturbation restart mechanism is introduced. When the cost function decreases by less than 1% for five consecutive iterations, a random perturbation following an N(0,0.1σ) distribution (σ being the current gradient norm) is applied near the current solution. The exploration and development capabilities are balanced by dynamically adjusting the perturbation intensity. Finally, the analysis increment is obtained and superimposed with the background field to obtain the analysis field, which simultaneously satisfies the atmospheric chemical transport equation and observation constraints. For parallel implementation, a master-slave architecture is adopted. The master process is responsible for cost function calculation and gradient aggregation, while the slave processes handle local calculations at different spatial grid points in parallel. Inter-process communication is achieved through MPI, improving the overall computational efficiency by more than 3 times.

[0072] As an optional embodiment of the present invention, the expression for obtaining the analysis increment of each spatial grid point in the simulation region is optionally:

[0073] ;in, express The incremental analysis express State variables, This represents the background field error covariance. Indicates assimilation The perturbation sample. Represents the eigenvector. Indicates the number of samples. Represents the identity matrix. Indicates simulation The perturbation sample. This indicates the observation increment.

[0074] As an optional embodiment of the present invention, optionally, in step S5, the corrected simulation region is obtained. Concentration and flux distribution data include:

[0075] S501. Based on the aforementioned analysis increment, a joint assimilation algorithm is used to simultaneously update the regional air quality model. Atmospheric concentration field and emission flux field ensure that the non-negativity constraint of concentration and the continuity constraint of flux are satisfied;

[0076] In step S501, it is important to explain in detail that the design of the joint assimilation algorithm must consider both the physical laws of atmospheric chemical transport and the stability of numerical calculations. Specifically, the analysis increment Δx is first decomposed into two parts: the concentration field increment Δc and the flux field increment Δf. The concentration field update uses a multiplicative correction method, ensuring that the concentration value remains non-negative through exponential transformation. The flux field update uses an additive correction method, and a flux gradient limiter is introduced to smooth grid points where the flux change rate exceeds a preset threshold (e.g., 0.5 μmol / (m²·s)). In numerical implementation, a semi-implicit time integration scheme is adopted, splitting the coupled update of the concentration and flux fields into a prediction step and a correction step. The prediction step first calculates the concentration field change trend based on the current flux field, and the correction step then adjusts the flux field distribution based on the concentration field feedback. Iterative convergence ensures that both satisfy the law of mass conservation. To handle boundary conditions, a buffer layer is introduced at the model boundary grid points, applying a damping coefficient (e.g., 0.2) to the flux increment flowing into the boundary to avoid the influence of external outliers on the model's internal structure. In terms of parallel computing optimization, the update operations of the concentration field and flux field are distributed to different computing nodes, and data exchange is achieved through non-blocking communication, which improves the overall computing efficiency by 2.3 times. The final corrected field satisfies the following physical constraints: the concentration field follows an exponential decay law in the vertical direction, and the flux field spatially corresponds to the land use type (e.g., wetland, paddy field). The consistency of the emission inventory exceeds 85%, and the time-series correlation between the model output and the daily variation peak of the observation station is greater than 0.7.

[0077] S502. Based on the updated regional air quality model, through an iterative optimization process, the analytical increments are applied to the state variables to correct [the data]. Concentration and flux distribution data, wherein the state variables include concentration and flux, are used with atmospheric transport equations as physical constraints to adjust the consistency between simulated and observed values;

[0078] In step S502, it is important to explain in detail that the core of the iterative optimization process lies in achieving a dynamic balance between simulated and observed values ​​through physical constraints. Specifically, the analysis increment is first decomposed into concentration field increments and flux field increments, which act on two dimensions of the state variables respectively. During the concentration field correction stage, explicit Euler method is used for time integration, and a relaxation factor (α = 0.3~0.6) is introduced to control the correction intensity, ensuring that the rate of concentration change does not exceed the allowable range of the atmospheric chemical transport equation. For example, for concentration updates in the tropospheric boundary layer, the relaxation factor is dynamically adjusted according to stability classification: α = 0.6 is used for unstable stratification to accelerate convergence, and α = 0.3 is used for stable stratification to avoid numerical oscillations.

[0079] In the flux field correction phase, a spatial constraint function based on land use type is constructed. For high-emission areas such as wetlands and paddy fields, prior information from the emission inventory is introduced to construct a cost function term, and the flux variation is limited to a reasonable range of ±20% using the Lagrange multiplier method. For urban areas, Kriging interpolation is used to construct spatial continuity constraints to ensure that the flux field gradient changes conform to the turbulent diffusion law. In practice, the horizontal grid is divided into several sub-regions, and each sub-region is assigned an independent flux adjustment coefficient. Spatial smooth transition is achieved through gradient matching conditions at the boundaries of adjacent regions.

[0080] The implementation of physical constraints relies on a dual verification mechanism: firstly, the atmospheric transport equation verification, which calculates the residuals between the corrected concentration field and flux field through backpropagation of the adjoint model. When the residual exceeds a threshold (e.g., 1e-4 ppm), the adjustment coefficient decay (β=0.8) is triggered; secondly, the mass conservation verification, which calculates the integral difference between the input flux and the output concentration for each grid point. When the difference exceeds the global standard, the mass conservation coefficient decay is triggered. When the budget reaches 5%, the gradient correction algorithm for the flux field is activated. This algorithm is based on the finite volume method and adjusts the flux distribution ratio between adjacent grid points.

[0081] In designing the iteration termination conditions, a triple criterion is adopted: the iteration is considered converged when the root mean square error (RMSE) between the concentration field and the observed values ​​decreases by less than 2% in three consecutive iterations, the spatial correlation (R²) between the flux field and the emission inventory increases by less than 0.05, and the number of physical constraint violations is less than 5 times per 10,000 grid points. To prevent local optima, a simulated annealing mechanism is introduced. When a plateau is reached (e.g., RMSE change <0.5% in 10 consecutive iterations), a suboptimal solution is accepted with a probability of 0.1, and the current solution space is broken through by perturbing the edge grid points of the flux field (accounting for 5%). The final corrected data achieves a spatiotemporal resolution of 1km×1km×1h, the agreement between the concentration field and the satellite inversion data (IOA index) is improved to 0.87, and the phase difference between the diurnal variation of the flux field and the eddy flux tower observations is reduced to within ±1 hour.

[0082] S503, Calculate the corrected simulation area Concentration and flux distribution data, including three-dimensional concentration fields at spatiotemporal resolution and emission flux fields at high spatial resolution;

[0083] In step S503, it is necessary to explain in detail that the calculation of the corrected simulation region... The process of obtaining concentration and flux distribution data requires comprehensive consideration of the consistency between spatiotemporal resolution and physical mechanisms. In practice, the construction of the three-dimensional concentration field adopts a layered processing strategy: the lower troposphere (0-2km) is divided into 50m vertical intervals, the free troposphere (2-10km) into 200m intervals, and the stratosphere (above 10km) into 1km intervals, ensuring the accuracy of vertical gradient capture. Horizontally, bilinear interpolation maps the 1km×1km grid data to the model's standard grid. Simultaneously, a terrain-following coordinate system is introduced to perform elevation correction for complex terrain areas such as mountains and basins, eliminating the influence of terrain shading on concentration distribution.

[0084] The calculation of high spatial resolution emission flux fields integrates multi-source inventory data. First, emission inventories from natural sources such as wetlands and paddy fields (e.g., EDGARv6.0) are spatially overlaid with anthropogenic source inventories from cities and industries (e.g., MEICv1.4) on a 5km×5km grid. Then, Kriging interpolation is used to improve the resolution to 1km×1km. To address inventory uncertainties, dynamic weighting coefficients are introduced: for satellite-covered areas, the weights are tilted towards the observation inversion results (weighting coefficient 0.7); for areas without observations, the weights rely on prior inventory values ​​(weighting coefficient 0.3). The temporal resolution of the flux field is achieved through diurnal dynamic adjustment. Combining the diurnal variation characteristics of land use types (e.g., peak emissions during rice irrigation), a time modulation function is applied to agricultural source fluxes to ensure that the diurnal variation of the flux field is synchronized with the actual emission process.

[0085] For spatiotemporal resolution matching, the concentration field and flux field are coupled and verified using the atmospheric transport equation. Specifically, a flux field with a 1-hour time step is input into a regional air quality model (e.g., CMAQv5.4). The evolution of the concentration field is calculated using the advection-diffusion equation, and the root mean square error (RMSE) of the model output is compared with that of the observed data. When the RMSE exceeds a threshold (e.g., 0.5 ppm), a flux field adjustment mechanism is triggered: for grid points in areas with high RMSE (error > 1 ppm), the flux value is reversed according to the error ratio (error value / threshold), with the correction range limited to ±15% to avoid over-adjustment leading to physical inconsistencies.

[0086] In the data post-processing stage, a mass conservation filtering algorithm is used to smooth the concentration field. This algorithm calculates the integral difference between the input flux and the output concentration at each grid point. When the difference exceeds the global average, the algorithm is applied. At 3% of the budget, flux field gradient correction is initiated: for grid points in the positive difference region (input > output), flux values ​​are reduced proportionally; for grid points in the negative difference region (input < output), flux values ​​are increased proportionally. The adjustment magnitude is determined through iterative convergence to ensure that the final concentration field satisfies the law of mass conservation. The final corrected data achieves a spatiotemporal resolution of 1 km × 1 km × 1 h, the agreement between the concentration field and satellite inversion data (IOA index) is improved to 0.89, the phase difference between the diurnal variation of the flux field and the eddy flux tower observations is reduced to within ±0.8 hours, and the correlation (R²) between the model output and the seasonal variation trend of the observation station exceeds 0.92.

[0087] S504, Verification and Correction of Simulation Area The physical rationality of concentration and flux distribution data.

[0088] In step S504, it is necessary to explain in detail that the verification correction is performed in the simulation area. The physical validity of concentration and flux distribution data requires the construction of a multi-dimensional verification system. First, a basic physical constraint verification should be performed, including a mass conservation check: for each grid point, calculate the integral difference between the input flux and the output concentration; when the difference exceeds the global average... At 3% of the budget, the flux field gradient correction algorithm is activated. This algorithm, based on the finite volume method, ensures that the total mass change rate of the region is less than 0.5% / day by adjusting the flux distribution ratio of adjacent grid points. Energy balance verification: This verifies the correlation between the vertical gradient of the concentration field and solar radiation intensity. When the daily variation in tropopause concentration exceeds 5 ppb, boundary layer parameter correction is triggered. Finally, chemical transformation rationality verification: This is done by comparing the output of the model... Oxidation products ( , The concentration and observed values ​​are used to ensure that the oxidation rate coefficient is within the range of (0.8~1.2)×theoretical value.

[0089] In terms of spatiotemporal continuity verification, a triple verification mechanism is adopted: spatial continuity verification calculates the second derivative of the concentration gradient between adjacent grid points, and Kriging interpolation smoothing is initiated when the gradient change rate exceeds 0.3 ppb / km²; temporal continuity verification compares the change rate of the concentration field over 3 consecutive hours, and triggers temporal smoothing filtering when the change amplitude exceeds 2 ppb / h; and vertical continuity verification verifies the smooth transition of the concentration field at the troposphere-stratospheric boundary, and adjusts the vertical diffusion coefficient when the concentration difference between the 2km and 3km altitude layers exceeds 15 ppb.

[0090] For the rationality verification of boundary conditions, a dual verification standard is set: the inflow boundary verification compares the model boundary flux with the output values ​​of a global atmospheric transport model (such as TM5), and when the difference exceeds 20%, a weighted average method is used to adjust the boundary value; and the outflow boundary verification verifies the top boundary of the model. When the deviation between net flux and the global budget value exceeds 5%, a global flux field correction is initiated.

[0091] Regarding extreme value testing, a dynamic threshold system is constructed: a three-level early warning mechanism is set up for concentration extreme value testing. When the concentration of a grid point exceeds three times the standard deviation of the background value, neighborhood average verification is initiated; when it exceeds five times the standard deviation, physical mechanism source analysis is conducted; and for flux extreme value testing, for grid points where the daily maximum flux exceeds twice the average value of the emission inventory, a land use type correction factor is introduced, with the flux upper limit in wetland areas limited to 1.5 times the inventory value and in urban areas limited to 1.2 times the inventory value.

[0092] The final verification results must simultaneously meet the following criteria: the mass conservation error is less than the global average. The corrected data is deemed physically reasonable when the annual emissions are 1%, the correlation (R²) between the model output and the seasonal variation trend of the observation station exceeds 0.9, the agreement between the concentration field and the satellite inversion data (IOA index) is not less than 0.85, and the phase difference between the flux field and the diurnal variation observed by the eddy flux tower is controlled within ±1 hour.

[0093] As an optional embodiment of the present invention, the expression of the joint assimilation algorithm is optionally: ;in, express Analysis field, express Site inspection express The incremental analysis.

[0094] As an optional embodiment of the present invention, optionally, in step S6, different grid points within the simulation area are estimated. The emission flux per unit time includes:

[0095] S601, in the corrected simulation area Concentration and flux distribution data are input into the atmospheric transport model in the regional air quality model, and combined with meteorological field data to calculate the net flux change per unit time for each grid point;

[0096] In step S601, it is necessary to explain in detail that the corrected simulation area... When inputting concentration and flux distribution data into the atmospheric transport model of a regional air quality model, quality-controlled reanalysis meteorological field data must be loaded simultaneously, including key parameters such as three-dimensional wind field, temperature field, pressure field, and boundary layer height. When calculating the net flux change per unit time for each grid point, a discretized form of the advection-diffusion equation is used: for each grid point, the meteorological field data is first matched to the desired values ​​using bilinear interpolation. A 1km × 1km grid was used for the concentration field, and then the horizontal and vertical fluxes were calculated, with the concentration gradient calculated using the central difference method. The diffusion term employed gradient transport theory, and the horizontal diffusion coefficient (Kh) and vertical diffusion coefficient (Kv) were determined using a boundary layer parameterization scheme; Kh and Kv decay exponentially with height within the stratosphere. The net flux change was the algebraic sum of the fluxes and the diffusion fluxes, with a time integration step of 1 hour. The concentration field maintains a consistent temporal resolution. For terrain-complex areas (such as mountains and canyons), a terrain-following coordinate system is introduced, and the terrain shading effect is eliminated by correcting the vertical diffusion coefficient. The final output net flux change data for each grid point includes three directional components and the total amount, with units uniformly set to μmol / (m²·s), and includes a mass conservation check flag: when the integral difference between the input flux and the output concentration at a grid point exceeds the global average... When the budget reaches 2%, it is marked as "needs correction" and triggers the subsequent flux field adjustment mechanism.

[0097] S602. Based on the net flux change, the boundary transport flux and chemical loss are integrated using the mass balance equation to separate the local emission contribution.

[0098] In step S602, it is important to explain in detail that the construction of the mass balance equation based on net flux change data needs to comprehensively consider three factors: horizontal boundary transport, vertical boundary transport, and chemical loss. Specifically, the net flux exchange at the regional boundaries is first calculated: Horizontally, the total horizontal transport is obtained by integrating the fluxes of the four lateral boundaries (east, west, south, and north), where inflow flux is taken as positive and outflow flux as negative; vertically, the vertical flux is calculated based on the concentration gradient and diffusion coefficient between the top and bottom layers of the model, with upward flux taking positive and downward flux taking negative. Chemical loss is determined by the model's built-in... The oxidation module is calculated based on the methane-hydroxyl radical reaction kinetics and dynamically adjusts the oxidation rate constant by combining temperature field data to ensure that the amount of chemical loss is consistent with the actual atmospheric chemical process.

[0099] The separation of local emission contributions employs an iterative correction method: initially assuming zero local emissions, a theoretical concentration field is calculated using the mass balance equation, and the residual field is obtained by comparing it with the actual observed concentration field. When the residual exceeds a threshold (e.g., 0.5 ppm), it is allocated to the local emission term according to the grid area ratio, forming a new emission estimate. This process is repeated until the root mean square error (RMSE) of the residual field is less than 0.1 ppm or the number of iterations reaches 20. For high-emission areas such as wetlands and paddy fields, prior constraints from the emission inventory are introduced. When the local emissions obtained through iteration exceed 1.5 times the average of the inventory, the correction magnitude is reduced proportionally to avoid unreasonable emission estimates.

[0100] The mass balance equation relies on a dual-check mechanism: the first is the flux closure check, which calculates the difference between the total regional input flux (boundary transport + local emissions) and the total output flux (chemical losses + concentration changes). When the difference exceeds the global... The first adjustment is a global adjustment of the emissions field, triggered when the budget reaches 1%. The second is a spatial continuity test, which calculates the gradient rate of change of local emissions at adjacent grid points. When the gradient exceeds 0.2 μmol / (m²·s·km), a Laplace smoothing operator is used for spatial filtering. The final local emission contribution data maintains a spatial resolution of 1km×1km, a temporal resolution consistent with the net flux change (1 hour), and satisfies the law of conservation of mass: the balance error between total regional emissions and boundary transport and chemical loss is less than 0.8%.

[0101] S603. The Bayesian optimization algorithm is used to minimize the residual between the simulated value and the satellite column concentration observation value, and the gridded emission flux is iteratively optimized.

[0102] In step S603, it is necessary to explain in detail that when using the Bayesian optimization algorithm, the objective function must first be constructed: the simulation area is divided into a 1km×1km grid, and the simulation is calculated for each grid point. The sum of squared residuals between column concentration and satellite observations is used as the optimization objective. The algorithm establishes a probabilistic mapping relationship between residuals and emission fluxes through a Gaussian process, where the covariance function adopts the Matern kernel function to balance spatial correlation and computational efficiency. During the iteration process, each sampling selects the emission flux combination with the highest posterior probability, prioritizing the exploration of regions with high residual values ​​(such as urban industrial areas and wetland peripheries). To avoid getting trapped in local optima, an adaptive noise term is introduced: when the residual decrease is less than 0.1% after 5 consecutive iterations, the search space is automatically expanded by 20%, while the sampling step size is dynamically adjusted. The initial step size is set at 30% of the emission inventory standard deviation, gradually decreasing to 10% with the number of iterations. In terms of constraints, a dual boundary is set: the physical rationality boundary requires that the emission flux at each grid point is not less than 80% of the lower limit of the natural source inventory and not more than 120% of the upper limit of the anthropogenic source inventory; the spatial continuity boundary is calculated by the emission gradient of adjacent grid points, and when the gradient exceeds 0.5 μmol / (m²·s·km), Laplace smoothing is forcibly applied. The optimization termination criteria employ a triple criterion: convergence is determined when the decrease in the sum of squared residuals is less than 0.05% for three consecutive iterations, the standard deviation of the posterior probability distribution is less than 0.1 times the initial value, and the number of samplings reaches 50. The final optimized emission flux field maintains a spatial resolution of 1km × 1km, a temporal resolution synchronized with the satellite observation period (typically daily), and passes cross-validation: the optimization results are substituted into independent observation datasets (such as ground stations and aerial surveys), requiring the correlation coefficient (R²) between simulated and observed values ​​to exceed 0.85 and the root mean square error (RMSE) to be less than 0.3ppm. For regions that fail validation, local re-optimization is initiated: a 3km × 3km sub-region is constructed centered on the grid points where the validation residual is greater than 0.5ppm, and the Bayesian algorithm is reapplied for fine-tuning, with the adjustment range limited to ±20%, ensuring that the final emission flux field simultaneously satisfies global optimization and local rationality.

[0103] S604, Output different grid points within the simulation area Emission flux per unit time.

[0104] In step S604, it is necessary to explain in detail the output of different grid points within the simulation area. When calculating emission fluxes per unit time, a multi-dimensional data output framework needs to be constructed. Firstly, in the spatial dimension, emission flux distribution maps are generated at a 1km×1km grid resolution. Each grid point contains flux components in three directions (east-west, south-north, and vertical) and the total amount, with units uniformly expressed as μmol / (m²·s). For terrain-complex areas (such as mountains and urban canyons), a terrain correction coefficient is added, calculated by comparing the vertical diffusion coefficients between the terrain-following coordinate system and the standard coordinate system. The temporal dimension output adopts a hierarchical architecture: the base layer provides hourly flux data, synchronized with the satellite observation cycle; the aggregation layer generates daily and monthly average flux products, calculated using a sliding window averaging method, with window lengths set to 24 hours and 30 days, respectively. Regarding data quality control, each time-level output data is accompanied by a quality indicator: when the grid point quality conservation error exceeds the global standard... A value of 1.5% of the budget is marked as "suspicious"; a failure rate exceeding 5% in boundary condition checks is marked as "requires review." The output format is compatible with multiple data standards: NetCDF format includes complete coordinate variables (longitude, latitude, altitude) and time variables, supporting CF metadata specifications; CSV format provides a planar table output of key fields, including grid center coordinates, flux values, quality indicators, and other core information. For visualization needs, a GeoTIFF format flux distribution layer is generated simultaneously. The color map uses a "volcano" color scale, with low-value areas (<50 μmol / (m²·s)) displayed in blue tones, high-value areas (>200 μmol / (m²·s)) displayed in red tones, and intermediate transition areas using a yellow-to-orange gradient. Regarding uncertainty quantification, each grid point's output data includes a 95% confidence interval, calculated using Monte Carlo simulation: applying ±10% random perturbation to input parameters (meteorological field, satellite observation errors, boundary conditions, etc.), and running 1000 simulations to statistically analyze the flux value distribution characteristics. For areas with abnormally high values ​​(such as wetlands and landfills), the confidence interval width is expanded to ±25%, and potential influencing factors are noted in the data description. The final output data product undergoes a dual validation mechanism: internal validation compares the aggregation consistency of data at different time resolutions, requiring the integral error between daily and hourly averages to be less than 5%; external validation employs triple cross-validation, comparing the results with ground-based eddy flux tower observations, aircraft aerial survey data, and the output of global atmospheric inversion models (such as TM5). The output data is considered valid when the correlation coefficient (R²) between all three and the simulated values ​​exceeds 0.8 and the root mean square error (RMSE) is less than 0.2 ppm. For areas that fail validation, a correction suggestion report is generated, indicating the specific grid points requiring additional observations and the recommended observation periods.

[0105] In response to the current high timeliness and high resolution To address the difficulty in obtaining emissions data, a comprehensive approach was proposed, employing an atmospheric boundary layer model, a regional air quality model, and a hybrid assimilation module. Emission source inversion correction system. Compared with existing technologies, atmospheric boundary layer models use the σ-z coordinate system, while other common atmospheric models use the σ-p coordinate system. This allows for a more detailed description of the complex small-scale turbulent activities within the atmospheric boundary layer, regardless of... Anthropogenic emissions, as well as the main transmission, exchange, and transformation processes, are concentrated within the atmospheric boundary layer and are significantly affected by small-scale turbulent activities. Regional air quality models driven by the σ-z coordinate system can fundamentally improve the model's ability to... The rationality of the representation. In addition, the hybrid assimilation module performs well in terms of the multi-source nature of observation information absorption, computational accuracy, and time cost. Therefore, coupling it with the aforementioned models can exhibit even better optimization results.

[0106] The verification of the system simulation performance has been verified in an academic article published in the SCI journal "Atmospheric Environment".

[0107] Relevant verification statistics are as follows Figure 5 As shown; Figure 5 For the output of the model system simulation Scatter plot of regression analysis between vertical column concentration and GOSAT satellite remote sensing data. The plot shows a generally positive linear relationship between the two, with a correlation coefficient greater than 0.95 at each point, indicating that the model can reproduce the data well. Spatial and temporal distribution of concentration.

[0108] Table 1. Model Output and GOSAT Satellite Remote Sensing Monitoring Statistical analysis of vertical column concentrations (unit: 10¹⁹ molec / cm²; time periods: January, April, July, and October 2020)

[0109]

[0110] C obs and C sim These are the monthly averages of observation and simulation, respectively; σ obs and σ sim 1 represents the standard deviation of the observations and the simulations, respectively; MB represents the mean deviation; RMSE represents the root mean square error; R represents the correlation coefficient; and N represents the sample size.

[0111] The monthly average difference between the simulation and observation results in Table 1 is between 0.21 and 0.32, the correlation coefficient is above 0.96, and the root mean square error is between 0.26 and 0.36. Overall, this shows that the simulation and observation results are in very good agreement.

[0112] like Figure 4As shown, the embodiment proposed in this case... The emission flux retrieval system mainly consists of a coupled atmospheric boundary layer model, a regional air quality model, and a mixing and assimilation module. The software system of this application includes: an atmospheric boundary layer model, a regional air quality model, and... Hybrid assimilation module. In this embodiment, GOSAT satellite remote sensing monitoring is employed. Column average concentration. Prior emission sources were identified using the global emissions inventory published by CAMS (Copernicus Atmosphere Monitoring Service).

[0113] The working principle of this application is as follows: Atmospheric boundary layer models are used to simulate meteorological field elements, acquiring meteorological parameters such as wind field, temperature, air pressure, and humidity. The time resolution can reach 1 hour, and the output data can be provided to regional air quality models, assimilation modules, and emission source processing algorithms. Secondly, the regional air quality model simulates... The evolution process in the atmosphere outputs the spatiotemporal distribution of its concentration in the atmosphere. This concentration is also output to... The assimilation module, combined with GOSAT satellite remote sensing monitoring Concentrations are analyzed by using a hybrid assimilation module to obtain incremental analysis data, and prior emission sources are then inverted and optimized. The inverted emission sources obtained using this method achieve hourly levels and high spatial resolution consistent with the model settings. The inverted emission sources are then input back into a regional air quality model to simulate the next timeframe. The spatial distribution of concentrations can be determined. By repeatedly calculating and iterating according to the above principle, high spatiotemporal resolution emission flux data can be obtained.

[0114] like Figure 3 and 4 As shown, the three main components of this embodiment are the atmospheric boundary layer model, the regional air quality model, and... Hybrid assimilation module. The atmospheric boundary layer model takes global reanalysis meteorological fields as input and outputs key meteorological elements: the three-dimensional spatiotemporal distributions of wind field, temperature, pressure, and humidity, as well as boundary layer parameters such as boundary layer height and turbulent diffusion coefficient. The regional air quality model takes meteorological elements and emission inventories as input and outputs... Spatial and temporal distribution of concentration. The assimilation module inputs meteorological elements. Simulated concentration, Observe the concentration, and output the data as follows Analyze the increments to optimize emission fluxes.

[0115] This embodiment runs on a Linux platform. After preparing the global reanalysis meteorological field, prior emission inventory, and observational data, the simulation time period, region, and spatiotemporal resolution of the model system are set. Then, the atmospheric boundary layer model is run to generate meteorological field elements, driving the regional air quality model simulation. Atmospheric concentration is obtained, and the results are input into the assimilation module. Combined with observational data, inversion optimization is carried out to obtain high spatiotemporal resolution emission fluxes, which further drive the regional air quality model to carry out numerical simulations for the next time step.

[0116] The key points of this invention are Hybridization systems and their coupling with boundary layer models. Currently, atmospheric models generally use the σ-p coordinate system, but the solution schemes for their fundamental atmospheric equations differ fundamentally from those of the σ-z coordinate system.

[0117] in The height of the model layer (meters) The top height of the layer set for the pattern. This refers to the ground elevation. The height is the natural geometric height. This formula expresses the vertical layer height distribution of the model in the σ-z coordinate system. The distribution will change according to the ground elevation, and by setting it, sufficient vertical grid density can be maintained within the boundary layer to achieve detailed description.

[0118] Therefore, coupling boundary layer models with hybrid assimilation systems requires rewriting the corresponding program code to represent the complex transmission calculations of numerous data streams, including multi-source observation data, meteorological element fields, chemical substance concentrations, and emissions, within the σ-z coordinate system. Based on this, a model is established by comprehensively utilizing dynamic localization radius, SVD decomposition, and ensemble sample-four-dimensional variational techniques. Hybrid assimilation system, ultimately improving Emissions assimilation and inversion effects. Therefore, the protection point of this invention is a system coupling program, which uses σ-z coordinate system grid points as a basis to couple and link the boundary layer model, regional air quality model, and mixed assimilation model to form... High-resolution inversion system. This system transmits meteorological element variables obtained from calculations in the σ-z coordinate system to regional air quality models and mixing assimilation systems to conduct simulations and corrections of chemical fields and emission sources.

[0119] Example 2

[0120] A high resolution The emission flux inversion system includes:

[0121] processor;

[0122] Memory used to store processor-executable instructions;

[0123] The processor is configured to achieve high resolution when executing executable instructions. Emission flux inversion method.

[0124] It should be noted that the computer device includes a processor, a memory, and may also include one or more of a multimedia component, an input / output (I / O) interface, and a communication component.

[0125] The processor controls the overall operation of the computer device to achieve the aforementioned high resolution. All or some of the steps in the emission flux inversion method.

[0126] Memory is used to store various types of data to support the operation of the computer device. This data may include, for example, instructions for any application or method used to operate on the computer device, as well as application-related data. Memory can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0127] The multimedia component may include a screen and an audio component, wherein the screen may be, for example, a touch screen, and the audio component is used to output and / or input audio signals; for example, the audio component may include a microphone for receiving external audio signals, the received audio signals may be further stored in memory or transmitted via a communication component; the audio component may also include at least one speaker for outputting audio signals.

[0128] I / O interfaces provide interfaces between the processor and other interface modules, such as keyboards, mice, buttons, etc.; these buttons can be virtual buttons or physical buttons.

[0129] The communication component is used for wired or wireless communication between the computer device and other devices; wireless communication, such as Wi-Fi, Bluetooth, Near Field Communication (NFC), 2G, 3G, 4G or 5G, or one or more combinations thereof, and the corresponding communication component may include: Wi-Fi module, Bluetooth module, NFC module, mobile communication module.

[0130] As a preferred embodiment, the computer device may be implemented using one or more application-specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field-programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors, or other electronic components to perform the aforementioned high-resolution... Emission flux inversion method.

[0131] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A high-resolution The emission flux inversion method is characterized by, The method includes: S1. Collect multi-source atmospheric data, including ground-based hourly data for different simulated regional locations in the regional air quality model. Observed values ​​and satellite remote sensing monitoring values; S2. Determine the distance weighting function based on the sparsity of the observed values ​​in time and space, and use the distance weighting function to calculate the correlation between each grid point in the simulation area and the randomly occurring observed values. S3. Based on the correlation, a set of samples is generated using a four-dimensional sliding sampling algorithm. The set of samples satisfies the actual physical constraints of the assimilated object. Then, singular value decomposition is used to reduce the dimensionality of the set of samples, reducing the computation and programming difficulty, while maintaining sufficient dispersion of the set of samples to cover all possible physical intervals of the state variables. S4. Based on the aforementioned sample set, a hybrid assimilation method is used to replace the tangent linear mode and the adjoint mode in the regional air quality model, and multi-source models are introduced. The cost function is solved using observational data to obtain the analysis increment for each spatial grid point in the simulated region; S5, for the aforementioned analytical increment pair Atmospheric concentrations and emission fluxes are jointly assimilated to obtain the corrected simulated region. Concentration and flux distribution data; S6, using the above Concentration and flux distribution data were used for Emission flux inversion to estimate different grid points within the simulation region Emission flux per unit time; Step S4 involves obtaining the analysis increment for each spatial grid point in the simulation region, including: S401. Estimate the background error covariance matrix based on the set of samples to quantify the uncertainty of the state variables; S402, Introduce the multi-source Observational data is used to calculate the difference between observed and simulated values, which is then used as the observation increment. S403, Based on the aforementioned multi-source The cost function is constructed from the observation data, including a background term and an observation term, where the background term is based on the background error covariance matrix and the observation term is based on the observation error covariance matrix. S404. Use an iterative optimization algorithm to find the minimum value of the solution cost function and obtain the analysis increment for each spatial grid point in the simulation area; The characteristic feature is that the expression for obtaining the analysis increment of each spatial grid point in the simulation region is: ； in, express The incremental analysis express State variables, This represents the background field error covariance. Indicates assimilation The perturbation sample. Represents the eigenvector. Indicates the number of samples. Represents the identity matrix. Indicates simulation The perturbation sample. This indicates the observation increment.

2. A high-resolution [property] as described in claim 1 The emission flux inversion method is characterized by, In step S1, the regional air quality model is a simulation model based on RAMS and CMAQ, used to simulate the three-dimensional meteorological field and atmospheric chemical transport processes.

3. A high-resolution [property] as described in claim 1 The emission flux inversion method is characterized by, The expression for the distance weight function in step S2 is: ， ； ； in, express and Distance weighting function between two points , express and The distance between two points Represents a constant given empirically. This represents the distance revision function.

4. A high-resolution [property] as described in claim 1 The emission flux inversion method is characterized by, In step S3, a set of samples is generated based on the correlation using a four-dimensional sliding sampling algorithm. The set of samples satisfies the actual physical constraints of the assimilation object, including: S301. Define the size and step size of the spatiotemporal sliding window, and extract state variable samples within the window according to the probability density distribution based on the correlation. S302. Apply physical constraints to the state variable samples within each sliding window. The physical constraints include concentration non-negativity, flux continuity, and atmospheric transport equation boundary conditions. S303. Combine the state variable samples that satisfy the constraints according to the time series to form a four-dimensional set sample matrix; S304. Perform singular value decomposition on the set sample matrix, retain the first k principal singular values ​​and corresponding singular vectors, and construct a dimension-reduced low-dimensional sample space, wherein the value of k is dynamically determined according to a preset variance contribution rate threshold. S305. Introducing a random perturbation term into the low-dimensional sample space can ensure coverage of all possible physical intervals of the state variable.

5. A high-resolution [property] as described in claim 1 The emission flux inversion method is characterized by, In step S5, the corrected simulation region is obtained. Concentration and flux distribution data include: S501. Based on the aforementioned analysis increment, a joint assimilation algorithm is used to simultaneously update the regional air quality model. Atmospheric concentration field and emission flux field ensure that the non-negativity constraint of concentration and the continuity constraint of flux are satisfied; S502. Based on the updated regional air quality model, through an iterative optimization process, the analytical increments are applied to the state variables to correct [the data]. Concentration and flux distribution data, wherein the state variables include concentration and flux, are used with atmospheric transport equations as physical constraints to adjust the consistency between simulated and observed values; S503, Calculate the corrected simulation area Concentration and flux distribution data, including three-dimensional concentration fields at spatiotemporal resolution and emission flux fields at high spatial resolution; S504, Verification and Correction of Simulation Area The physical rationality of concentration and flux distribution data.

6. A high-resolution [property] as described in claim 5 The emission flux inversion method is characterized by, The expression for the joint assimilation algorithm is: ； in, express Analysis field, express Site inspection express The incremental analysis.

7. A high-resolution [property] as described in claim 1 The emission flux inversion method is characterized by, In step S6, the different grid points within the simulation region are estimated. The emission flux per unit time includes: S601, in the corrected simulation area Concentration and flux distribution data are input into the atmospheric transport model in the regional air quality model, and combined with meteorological field data to calculate the net flux change per unit time for each grid point; S602. Based on the net flux change, the boundary transport flux and chemical loss are integrated using the mass balance equation to separate the local emission contribution. S603. The Bayesian optimization algorithm is used to minimize the residual between the simulated value and the satellite column concentration observation value, and the gridded emission flux is iteratively optimized. S604, Output different grid points within the simulation area Emission flux per unit time.

8. A high-resolution The emission flux inversion system is characterized by, The system includes: processor; Memory used to store processor-executable instructions; The processor is configured to achieve the high resolution of any one of claims 1 to 7 when executing the executable instructions. Emission flux inversion method.

Images