An industrial source emission inventory spatial downscaling method

Abstract

The application discloses an industrial source emission inventory spatial downscaling method, and belongs to the technical field of atmospheric pollution source inventory compilation and satellite remote sensing application. The method comprises the following steps: obtaining a coarse resolution industrial source emission inventory of a target region, thermal infrared remote sensing data and satellite atmospheric column concentration observation data; based on the thermal infrared remote sensing data, a thermal infrared radiation flux is inversed, and a physical lower boundary constraint of emission intensity is constructed according to the principle of thermodynamic energy conservation; based on the analytical solution of the Gaussian plume model, an inversion equation is constructed, and the physical lower boundary constraint is taken as an inequality constraint, so that the spatial distribution of the emission intensity of the target region is inversed and solved; and a high-resolution industrial source emission intensity distribution is output. Through the application, the problem that the thermal infrared data is only used as a spatial position weight and cannot utilize the thermodynamic quantitative relationship between the thermal infrared data and the emission intensity in the prior art is solved, the role of the thermal infrared data is changed from a statistical agent to a physical constraint, and the physical accuracy of the downscaling result is improved.

Description

Technical Field

[0001] This invention relates to a spatial downscaling method for industrial source emission inventories, belonging to the technical field of atmospheric pollution source inventory compilation and satellite remote sensing applications. Background Technology

[0002] Spatial downscaling of industrial emission inventories is a key technology for allocating low-resolution total emissions onto a high-resolution grid, and is an important foundation for carbon emission accounting and air quality simulation. Current mainstream downscaling techniques typically utilize socioeconomic proxy factors such as nighttime light pollution, population density, and land use, establishing the relationship between these proxy factors and emission intensity through statistical regression or machine learning models. In recent years, thermal infrared remote sensing data has been introduced into this field due to its ability to identify industrial heat sources, and is used to assist in identifying emission areas.

[0003] However, current technologies typically utilize thermal infrared data only as an indicator of spatial location, i.e., as an independent variable in statistical regression or a soft prior weight in an inversion model, to tell the model "where emissions might occur." This approach fails to uncover the intrinsic thermodynamic quantitative relationship between thermal infrared radiation flux and emission intensity. Because waste heat emissions from industrial heat sources are physically coupled with the fuel combustion process, the radiative heat dissipation component observed by thermal infrared remote sensing is actually a physical lower bound indicator of emission intensity. Current technologies ignore this physical mechanism, resulting in a lack of physical constraints in the downscaling process, weak physical consistency of the results, and difficulty in accurately quantifying emission intensity.

[0004] Therefore, a new solution is needed to address this problem. Summary of the Invention

[0005] The technical problem to be solved by this invention is to provide a spatial downscaling method for industrial source emission inventories. This method solves the problem in the prior art that thermal infrared data is only used as a spatial location weight and fails to utilize its thermodynamic quantitative relationship with emission intensity. It realizes the transformation of thermal infrared data from a statistical proxy to a physical constraint, thereby improving the physical accuracy of the downscaling results.

[0006] The technical problem to be solved by this invention is achieved by the following technical solution:

[0007] A spatial downscaling method for industrial source emission inventories includes the following steps:

[0008] Obtain a coarse-resolution industrial source emission inventory of the target area, and acquire simultaneous thermal infrared remote sensing data and satellite atmospheric column concentration observation data;

[0009] Based on the aforementioned thermal infrared remote sensing data, thermal infrared radiation flux is retrieved, and a physical lower bound constraint for emission intensity is constructed according to the principle of thermodynamic energy conservation. The physical lower bound constraint is constructed as follows: thermal efficiency, radiative heat dissipation ratio, and emission factor are obtained. The radiative heat dissipation ratio is the proportion of radiative heat dissipation in waste heat. The waste heat ratio is obtained by subtracting the thermal efficiency from 1. The product of the waste heat ratio and the radiative heat dissipation ratio is calculated as a first product. The thermal infrared radiation flux is divided by the first product to obtain an estimated total fuel combustion energy value. The estimated total fuel combustion energy value is multiplied by the emission factor to obtain a lower bound value for emission intensity, and the emission intensity is not less than the lower bound value for emission intensity.

[0010] An inversion equation is constructed based on the analytical solution of the Gaussian plume model. The satellite atmospheric column concentration observation data is used as the input of the inversion equation, and the physical lower bound constraint is used as the inequality constraint to invert and solve the spatial distribution of emission intensity in the target area.

[0011] Based on the spatial distribution of emission intensity obtained from the inversion solution, the coarse-resolution industrial source emission inventory is spatially downscaled to output a high-resolution industrial source emission intensity distribution.

[0012] Optionally, the step of retrieving thermal infrared radiative flux based on the thermal infrared remote sensing data specifically includes:

[0013] Radiometric calibration and atmospheric correction were performed on the thermal infrared remote sensing data to obtain brightness and temperature data;

[0014] According to Planck's law, the brightness temperature data is integrated within a preset thermal infrared band to obtain the surface radiation flux, and the surface radiation flux is determined as the thermal infrared radiation flux.

[0015] Alternatively, when the thermal infrared remote sensing data is a VIIRS fire radiation power product, the radiation power value of the VIIRS fire radiation power product can be directly determined as the thermal infrared radiation flux.

[0016] Optionally, before applying the inequality constraints to the inversion equations, the method further includes a step of determining the set of meshes where the physical constraints are effective, specifically including:

[0017] The initial solution of the emission intensity is obtained by solving the inversion equation without applying the aforementioned inequality constraints.

[0018] Calculate the thermal infrared signal-to-noise ratio for each grid and calculate the constraint confidence index, which is the ratio of the initial solution to the lower limit of the emission intensity.

[0019] When the signal-to-noise ratio of the thermal infrared signal is greater than the preset signal-to-noise ratio threshold and the constraint confidence index is less than or equal to 1, the corresponding mesh is included in the set of meshes in effect for physical constraints.

[0020] In the step of using the physical lower bound constraint as an inequality constraint, the inequality constraint is applied only to the meshes in the set of meshes where the physical constraint is effective.

[0021] Optionally, the step of inverting and solving the spatial distribution of emission intensity in the target area specifically includes:

[0022] The study area was discretized into high-resolution grid cells, and the Gaussian plume transfer matrix was calculated based on meteorological data.

[0023] A regularization objective function with inequality constraints is constructed, wherein the regularization objective function includes a data fitting term, a spatial smoothing regularization term, and a thermal infrared prior regularization term;

[0024] The emission intensity vector is obtained by solving the convex quadratic programming problem with inequality constraints using the effective set method or the interior point method.

[0025] Optionally, the inversion solution process may also include a parameter iterative calibration step:

[0026] Input the spatial distribution of emission intensity obtained from the current iteration into the Gaussian plume forward model to simulate the satellite column concentration;

[0027] Calculate the relative deviation between the simulated satellite column concentration and the observed satellite atmospheric column concentration;

[0028] When the relative deviation is greater than the preset convergence threshold, the thermal efficiency, the emission factor and the proportion of radiative heat dissipation are adjusted within a physically reasonable range, and the lower limit of emission intensity is recalculated. The process returns to the step of performing the inversion solution of the spatial distribution of emission intensity in the target area and performs the next iteration solution.

[0029] In a single iteration, the adjustment range of each parameter does not exceed the preset first adjustment threshold, and the cumulative adjustment range does not exceed the preset second adjustment threshold.

[0030] Optionally, the parameter iterative calibration process also includes steps for constraint conflict detection and anomaly source identification:

[0031] Monitor whether, after parameter adjustment, the current iterative solution is less than the recalculated lower bound of the emission intensity in the grid.

[0032] If it exists, then calculate the constraint conflict degree, which is the ratio of the difference between the lower bound of the emission intensity and the current iterative solution to the lower bound of the emission intensity.

[0033] When the constraint conflict degree is greater than the preset conflict threshold, the corresponding mesh is marked as a constraint abnormal mesh;

[0034] For the constrained abnormal grid, the data quality of its thermal infrared remote sensing data is detected. If the data quality is lower than the preset quality standard, the thermal infrared data is determined to be unreliable. If the data quality meets the preset quality standard, it is determined to be a candidate region of non-thermal radiation dominant emission source, and the grid is removed from the set of grids with effective physical constraints in subsequent iterations.

[0035] Based on the constraint conflict degree, the data quality of the thermal infrared remote sensing data, the quantity and quality of satellite observations, and the iterative convergence status, the confidence level of each grid is comprehensively evaluated, and the confidence level is correlated with the spatial distribution of emission intensity and output; the confidence level includes at least high confidence, medium confidence and low confidence.

[0036] Optionally, the thermal efficiency, the proportion of radiative heat dissipation, and the emission factor are determined in the following ways:

[0037] The industry type of the target grid is determined based on thermal infrared temperature characteristics, land use data, and / or point of interest data;

[0038] Based on the industry type, the corresponding range of thermal efficiency, range of radiative heat dissipation ratio, and range of emission factor are searched from the preset parameter library;

[0039] In this process, a conservative strategy is adopted when selecting parameters, with the radiative heat dissipation ratio taking the upper limit of its preset range and the thermal efficiency taking the lower limit of its preset range.

[0040] During the parameter iterative calibration process, when it is necessary to adjust parameters first, the proportion of radiative heat dissipation should be adjusted first to keep the thermal efficiency and the emission factor relatively stable.

[0041] Optionally, the industry types include iron and steel smelting, thermal power, cement, and chemical industries;

[0042] The preset parameter library stores the mapping relationship between the industry type and the range of thermal efficiency and the range of radiative heat dissipation ratio; the range of thermal efficiency and the range of radiative heat dissipation ratio are determined according to the typical process temperature and heat dissipation characteristics of the industry.

[0043] Optionally, the thermal infrared remote sensing data includes remote sensing image data of at least one thermal infrared band;

[0044] The satellite atmospheric column concentration observation data includes TROPOMI NO2 tropospheric column concentration products and OCO-2 / 3XCO2 products;

[0045] The inversion solution adopts a multi-species joint inversion strategy, using the TROPOMI NO2 tropospheric column concentration product to constrain the spatial distribution structure of emissions, and using the OCO-2 / 3 XCO2 products to constrain the total emission level.

[0046] Optionally, the preset convergence threshold includes an absolute convergence threshold and a relative convergence threshold;

[0047] When the relative deviation is less than the absolute convergence threshold, it is determined to be converged;

[0048] Alternatively, convergence is determined when the relative deviation change rate between two consecutive iterations is less than the relative convergence threshold.

[0049] If convergence is not achieved after reaching the preset maximum number of iterations, output the result with the smallest relative deviation in the history of iterations.

[0050] The beneficial effects of this invention are: by using thermal infrared radiation flux to back-calculate the minimum value of total fuel combustion energy, a physical lower bound for emission intensity is established, and this lower bound is embedded as an inequality constraint in the inversion solution process. This approach upgrades the application of thermal infrared data from a simple statistical correlation to a physical constraint relationship, effectively avoiding unreasonable low values ​​in the inversion results in active heat source regions, and significantly improving the physical consistency and inversion accuracy of the emission inventory spatial downscaling. Attached Figure Description

[0051] Figure 1 This is a flowchart of the main method of the present invention.

[0052] Figure 2 This is a flowchart of the physical lower bound constraint construction process of the present invention.

[0053] Figure 3 This is a flowchart of the parameter iterative calibration and conflict detection process of the present invention. Detailed Implementation

[0054] To facilitate a clear understanding of the technical means, creative features, objectives, and effects of this invention, the invention will be further described below in conjunction with specific illustrations.

[0055] To address the technical problem in existing technologies where thermal infrared data is only used for spatial location weights in emission retrieval and fails to utilize the thermodynamic quantitative relationship between thermal infrared radiation flux and emission intensity, this application provides a spatial downscaling method for industrial source emission inventories, referencing... Figure 1 This includes the following steps.

[0056] S1, Data Acquisition and Preprocessing.

[0057] Obtain a coarse-resolution industrial emission inventory for the target area, which can be a MEIC 27 km resolution industrial emission inventory, and extract gridded total emission data for the target area from it. Acquire concurrent thermal infrared remote sensing data and satellite atmospheric column concentration observation data.

[0058] Thermal infrared remote sensing data can be nighttime thermal infrared imagery from the SDGSAT-1 Thermal Infrared Spectrometer (TIS), comprising three thermal infrared bands with ranges of 8.0–10.5 μm, 10.3–11.3 μm, and 11.5–12.5 μm, respectively, and a spatial resolution of 30 meters. Nighttime observations effectively avoid interference from solar radiation. Fire radiation power data from the VIIRS 375-meter resolution thermal anomaly product can also be used as a supplement. Radiometric calibration and atmospheric correction are performed on the thermal infrared imagery to retrieve brightness and temperature data for each band.

[0059] Satellite atmospheric column concentration observation data include the TROPOMI NO2 tropospheric column concentration product with a spatial resolution of 3.5 x 5 km, and the OCO-2 or OCO-3 XCO2 product with a spatial resolution of approximately 2 km.

[0060] It also acquires reanalysis meteorological data, including wind speed, wind direction, boundary layer height, and atmospheric stability, as well as auxiliary information such as land use data and points of interest data.

[0061] S2, Thermal infrared radiation flux inversion and physical lower bound constraint construction.

[0062] refer to Figure 2 In conventional emission spatial allocation methods, thermal infrared data is only used as a soft prior weight for spatial location, and the thermodynamic quantitative relationship between thermal infrared data and emission intensity is not utilized. This step utilizes the principle of thermodynamic energy conservation to construct thermal infrared data as a physical lower bound constraint for emission intensity.

[0063] S21, thermal infrared radiation flux inversion.

[0064] For the i-th grid, the surface temperature T_s,i is retrieved from the brightness temperature data. The radiative flux is calculated using Planck's blackbody radiation formula, combined with the surface emissivity ε_i and the retrieved surface temperature T_s,i, within three preset thermal infrared bands. After performing integration operations on each component and summing them, the surface radiant flux is obtained and determined as the thermal infrared radiant flux F_IR,i. When the thermal infrared remote sensing data is the VIIRS fire radiant power product, the radiant power value of the product is directly determined as the thermal infrared radiant flux F_IR,i, without the need for brightness-temperature inversion and band integration operations.

[0065] S22, Construction of physical lower bound constraints.

[0066] Of the total chemical energy released by fuel combustion, a portion is converted into useful work, while the remainder is dissipated as waste heat through various pathways, including radiation, convection, and conduction. Thermal infrared remote sensing captures the radiative heat dissipation component of this waste heat. Let the total fuel combustion energy of the i-th grid be E_fuel,i, and the thermal efficiency η_th,i be the proportion of useful work to the total fuel combustion energy. Then the waste heat proportion is 1-η_th,i. Let β_i be the proportion of radiative heat dissipation in the waste heat, then the radiative heat dissipation flux F_IR,i is expressed as:

[0067]

[0068] This leads to the total energy released during fuel combustion:

[0069]

[0070] Let α_i be the emission factor per unit fuel energy, then the emission intensity Q_i of the i-th grid is:

[0071]

[0072] This relationship is used as a physical lower bound constraint, meaning that the emission intensity is not less than the lower bound value of the emission intensity:

[0073]

[0074] Where Q_i^{lower} is the lower bound of the emission intensity of the i-th grid, F_IR,i is the thermal infrared radiation flux of the i-th grid, α_i is the emission factor corresponding to unit fuel energy, η_th,i is the thermal efficiency, and β_i is the proportion of radiative heat dissipation in waste heat. The physical meaning of this constraint is: for a grid with a significant thermal infrared radiation signal, its emission intensity cannot be lower than the emission corresponding to the minimum fuel consumption implied by that thermal infrared radiation flux.

[0075] A conservative strategy is adopted in parameter selection to ensure the strictness of the physical lower bound constraint. β_i is taken as the upper bound of its preset range, assuming that the proportion of radiative heat dissipation to waste heat is as large as possible; η_th,i is taken as the lower bound of its preset range, assuming that the thermal efficiency is as low as possible and the proportion of waste heat is as large as possible. This makes the calculated value of Q_i^{lower} smaller, thereby ensuring that the actual emission intensity satisfies Q_i ≥ Q_i^{lower} with a high probability.

[0076] S23, Industry type identification and parameter determination.

[0077] Thermal efficiency, radiative heat dissipation ratio, and emission factor are related to industry type, and the corresponding parameter range needs to be found from the preset parameter library according to the industry type of the target grid.

[0078] The industry type of the target grid is determined based on thermal infrared temperature characteristics and land use data or point-of-interest (POI) data. Specifically, preliminary temperature classification is performed based on the absolute value of single-band brightness temperature, combined with thermal infrared spatial morphology characteristics for auxiliary determination, and cross-validation is conducted using land use data or POI data to comprehensively determine the industry type. Alternatively, the industry type can be determined directly based on land use data or POI data.

[0079] The industry types include steel smelting, thermal power, cement, and chemical industries. The preset parameter library stores the mapping relationship between industry types and thermal efficiency ranges and radiative heat dissipation ratio ranges. The thermal efficiency range and radiative heat dissipation ratio range are determined based on the typical process temperatures and heat dissipation characteristics of each industry.

[0080] For the steel smelting industry, the typical temperature characteristic is high, exceeding 800 degrees Celsius, with a thermal efficiency ranging from 0.50 to 0.65 and a radiative heat dissipation ratio ranging from 0.30 to 0.50. For the thermal power industry, the typical temperature characteristic is medium-high, between 400 and 800 degrees Celsius, with a thermal efficiency ranging from 0.35 to 0.45 and a radiative heat dissipation ratio ranging from 0.15 to 0.25. For the cement industry, the typical temperature characteristic is high, exceeding 800 degrees Celsius, with a thermal efficiency ranging from 0.40 to 0.55 and a radiative heat dissipation ratio ranging from 0.25 to 0.40. For the chemical industry, the typical temperature characteristic is medium, between 200 and 600 degrees Celsius, with a thermal efficiency ranging from 0.45 to 0.60 and a radiative heat dissipation ratio ranging from 0.15 to 0.30.

[0081] When selecting parameters, the proportion of radiative heat dissipation is taken as the upper limit of its preset range, and the thermal efficiency is taken as the lower limit of its preset range, to ensure the strictness of the physical lower limit constraint.

[0082] S3, Unconstrained pre-solution and constraint-effective mesh determination.

[0083] Before applying the physical lower bound constraint as an inequality constraint to the inversion equation, an unconstrained pre-solution is first performed to obtain an initial solution for determining the constraint confidence.

[0084] The inversion equation is solved without applying the aforementioned inequality constraints to obtain an initial solution for emission intensity. The objective function of this unconstrained inversion problem contains only data fitting and spatial smoothing terms, and does not include thermal infrared constraint terms.

[0085] Calculate the signal-to-noise ratio (SNR) of the thermal infrared signal for each grid cell. This SNR is equal to the thermal infrared radiant flux divided by the noise level of the thermal infrared data. Simultaneously, calculate the constraint confidence index, which is the ratio of the initial solution to the lower bound of the emission intensity.

[0086] When the signal-to-noise ratio (SNR) of the thermal infrared signal is greater than a preset SNR threshold, and the constraint confidence index is less than or equal to 1, the corresponding mesh is included in the set of meshes where the physical constraints are effective. A typical value for the preset SNR threshold is 3.

[0087] When the initial solution already satisfies the physical lower bound (i.e., the constraint confidence index is greater than 1), the constraints of the grid are inactive and have no effect. When the initial solution is below the physical lower bound (i.e., the constraint confidence index is less than or equal to 1), the constraints of the grid are active, and the solution needs to be raised above the lower bound.

[0088] In subsequent steps where the physical lower bound constraint is used as the inequality constraint, the inequality constraint is applied only to the meshes in the set of meshes where the physical constraint is effective.

[0089] S4, Gaussian inversion with physical constraints.

[0090] Based on the analytical solution of the Gaussian plume model, an inversion equation is constructed. Satellite atmospheric column concentration observation data is used as the input to the inversion equation, and physical lower bound constraints are used as inequality constraints to solve for the spatial distribution of emission intensity in the target area.

[0091] S41, Discretization of the Gaussian plume forward model.

[0092] Based on the steady-state analytical solution of the Gaussian plume model, the concentration contribution of each grid cell to the satellite pixel is calculated. The study area is discretized into N high-resolution grid cells, and the emission intensity of each grid cell is Q_j. For the k-th satellite pixel, its column concentration observation C_k is represented as a linear superposition of the contributions from all grid cells:

[0093]

[0094] Where G_{kj} is the element in the k-th row and j-th column of the Gaussian plume transfer matrix, representing the contribution of the emission unit of the j-th grid cell to the column concentration of the k-th satellite pixel. This transfer matrix is ​​calculated based on the steady-state analytical solution of the Gaussian plume model, meteorological data, and grid spatial location; ε_k is the observation noise. Written in matrix form:

[0095]

[0096] Where C_obs is the satellite plume concentration observation vector, Q is the emission intensity vector to be inverted, and G is the Gaussian plume transfer matrix.

[0097] S42, Construction of regularized objective function with inequality constraints.

[0098] Construct a regularization objective function with inequality constraints:

[0099]

[0100] The constraints are:

[0101]

[0102] The first term is the data fitting term, which requires the inversion results to be consistent with satellite observations. This means that the emission distribution obtained from the inversion, after being forward modeled using a Gaussian plume, should approximate the column concentration actually observed by the satellite as closely as possible. The second term is the spatial smoothing regularization term, where L is the spatial smoothing matrix and λ_1 is the smoothing regularization parameter. This term constrains the continuous variation of emission intensity between adjacent grids, avoiding unreasonable and drastic fluctuations. The third term is the thermal infrared soft prior regularization term, where P_hs is the thermal infrared spatial prior and λ_2 is the prior constraint regularization parameter. This term guides the inversion results to be consistent with the thermal infrared spatial distribution when the confidence level of the thermal infrared data is high. The inequality constraint is a physical lower bound constraint, which only applies to grids in the set Ω_active where the physical constraints are effective, ensuring that the emission intensity of these grids is not lower than their thermodynamic lower bound value.

[0103] The thermal infrared spatial prior, P_hs, is specifically a spatial weight vector constructed based on the thermal infrared radiation flux. This vector is obtained by normalizing the thermal infrared radiation flux of each grid: P_hs,i = F_IR,i / max(F_IR), where max(F_IR) is the maximum value of the thermal infrared radiation flux of all grids within the study area. For grids where the thermal infrared signal is below the noise level, P_hs,i is set to 0. This construction method allows regions with high thermal infrared radiation flux to receive higher weights in the soft prior constraints, guiding the inversion results to maintain spatial consistency with the distribution pattern of thermal infrared radiation flux. The regularization parameters λ_1 and λ_2 are determined using the L-curve method or the generalized cross-validation method.

[0104] S43, Solving convex quadratic programming problems and multi-species joint inversion.

[0105] This optimization problem is a convex quadratic programming problem with linear inequality constraints, which can be solved using the effective set method or the interior point method. The solution process is as follows: Initialize Q to the initial solution obtained in step S3; check whether the constraint Q_i ≥ Q_i^{lower} is satisfied; for meshes that do not satisfy the constraint, add the constraint as the equality constraint Q_i = Q_i^{lower} to the effective set; solve the equality-constrained quadratic programming problem on the effective set; check whether there are any new meshes that violate the constraints, and if so, add them to the effective set and solve again; repeat until all constraints are satisfied.

[0106] The inversion solution employs a multi-species joint inversion strategy. The TROPOMI NO2 tropospheric column concentration product and either OCO-2 or OCO-3 XCO2 product are simultaneously used as observation inputs to construct a joint objective function. The joint objective function is a weighted sum of data fitting terms:

[0107]

[0108] Here, γ is a weighting coefficient, determined based on the signal-to-noise ratio of NO2 and XCO2 observation data. The spatial distribution structure of emissions is constrained by the Tropomi NO2 tropospheric column concentration product because NO2 has high spatial resolution and a strong correlation with industrial combustion emissions; the total emission level is constrained by OCO-2 or OCO-3 XCO2 products because CO2 is a major component of industrial emissions and column concentration observations have high accuracy. Q_NOx and Q_CO2 are correlated through the industry emission factor α, i.e., Q_CO2 = α·Q_NOx, allowing both sets of observation data to jointly constrain the solution of the same emission field Q.

[0109] S5, Parameter Iteration Calibration and Convergence Determination.

[0110] refer to Figure 3 The inversion solution process also includes a parameter iteration calibration step to match satellite observations and improve inversion accuracy.

[0111] S51, forward modeling and verification.

[0112] Inputting the spatial distribution of emission intensity obtained from the current iteration into the Gaussian plume forward model, we can simulate the satellite column concentration values:

[0113]

[0114] in, This is the simulated satellite column concentration vector.

[0115] S52, Deviation Calculation and Convergence Judgment.

[0116] Calculate the relative deviation between the simulated satellite column concentration and the observed satellite atmospheric column concentration data. This relative deviation is defined as the norm of the difference between the simulated and observed values ​​divided by the norm of the observed values.

[0117]

[0118] Where δ is the relative deviation, representing the degree of deviation between the simulated value and the observed value.

[0119] The preset convergence thresholds include an absolute convergence threshold and a relative convergence threshold. Convergence is determined when the relative deviation is less than the absolute convergence threshold. A typical value for the absolute convergence threshold is 5%. Alternatively, convergence is determined when the rate of change of the relative deviation between two consecutive iterations is less than the relative convergence threshold. A typical value for the relative convergence threshold is 1%.

[0120] S53, parameter adjustment and iterative loop.

[0121] When the relative deviation exceeds the preset convergence threshold, the thermal efficiency, emission factor, and radiative heat dissipation ratio are adjusted within a physically reasonable range, and the lower limit of emission intensity is recalculated. The adjustment range of each parameter in a single iteration does not exceed the preset first adjustment threshold, and the cumulative adjustment range does not exceed the preset second adjustment threshold. A typical value for the first adjustment threshold is 10%, and a typical value for the second adjustment threshold is 20%.

[0122] During parameter iterative calibration, when it is necessary to adjust parameters first, the proportion of radiative heat dissipation should be adjusted first to keep thermal efficiency and emission factor relatively stable, because the proportion of radiative heat dissipation is the parameter with the greatest uncertainty among the three.

[0123] If convergence is not achieved after reaching the preset maximum number of iterations, output the result with the smallest relative deviation in the historical iterations. A typical value for the maximum number of iterations is 10.

[0124] S54, Constraint Conflict Detection and Anomaly Source Identification.

[0125] During the parameter iterative calibration process, the lower bound of emission intensity is recalculated after each parameter update, and it is monitored whether there are grids where the current iterative solution is less than the recalculated lower bound of emission intensity.

[0126] If a conflict exists, the grid is marked as a conflict grid, and the conflict degree is calculated. The conflict degree is defined as the ratio of the difference between the lower bound of emission intensity and the current iterative solution to the lower bound of emission intensity.

[0127]

[0128] in, Let be the constraint conflict degree of the i-th grid, representing the degree to which the lower bound of the emission intensity exceeds the current iterative solution. When the constraint conflict degree is greater than a preset conflict threshold, the corresponding grid is marked as a constraint anomalous grid. A typical value for the preset conflict threshold is 0.1.

[0129] For constrained anomalous grids, the data quality of their thermal infrared remote sensing data is assessed. If the data quality is below a preset quality standard, such as high cloud cover or low signal-to-noise ratio, the thermal infrared data is deemed unreliable. If the data quality meets the preset quality standard, it is identified as a candidate region for non-thermal radiation-dominated emission sources. In subsequent iterations, this grid is removed from the set of physically constrained grids, degenerating it into a soft prior constraint mode.

[0130] The confidence level of each grid is comprehensively evaluated based on the degree of constraint conflict, the data quality of thermal infrared remote sensing data, the quantity and quality of satellite observations, and the iterative convergence status. The confidence level includes at least high confidence, medium confidence, and low confidence. The confidence level is then correlated with the spatial distribution of emission intensity, and the confidence level is labeled for each grid.

[0131] S6, spatial downscaling output.

[0132] Based on the spatial distribution of emission intensity obtained from the inversion solution, the coarse-resolution industrial source emission inventory is spatially downscaled. Specifically, using the spatial distribution of emission intensity obtained from the inversion solution as spatial allocation weights, the total amount of each grid in the coarse-resolution inventory is allocated to the high-resolution grid according to the weights, outputting a high-resolution industrial source emission intensity distribution. The spatial resolution of this output can reach 100 meters.

[0133] Example 1: Taking a concentrated steel industry area as an example.

[0134] This embodiment uses a concentrated area of ​​the steel industry as an application scenario to verify the effectiveness of the present invention.

[0135] During the data acquisition and preprocessing phase, the MEIC 2020 industrial CO2 emission inventory for the region at a resolution of 27 km was obtained, and regional grid data was extracted. SDGSAT-1 TIS nighttime thermal infrared imagery from July to September 2020 was acquired, including six cloud-free images. Radiometric calibration and atmospheric correction were performed, and brightness and temperature in the three thermal infrared bands were retrieved. Tropospheric column concentration products (TROPOMINO2) and OCO-2 XCO2 products were acquired concurrently. ERA5 reanalysis meteorological data and point-of-interest (POI) data for the region were also acquired.

[0136] In the stage of thermal infrared radiation flux inversion and physical lower bound construction, based on the brightness temperature of the three bands of SDGSAT-1 TIS, the thermal infrared radiation flux of each grid is calculated by integrating over the three bands according to Planck's law. Combining the data of points of interest and thermal infrared temperature characteristics, the main industry type in this area is identified as iron and steel smelting. Iron and steel industry parameters are retrieved from the parameter library: thermal efficiency is set to 0.55, emission factor is set to 1.1 x 10^-7 kg / joule, and radiative heat dissipation ratio is set to 0.40. A conservative strategy is adopted, with the radiative heat dissipation ratio taking the upper bound of the parameter range and the thermal efficiency taking the lower bound of the parameter range. The physical lower bound of emission intensity for each grid is calculated.

[0137]

[0138] In the unconstrained pre-solution stage, a simplified problem containing only data fitting terms and spatial smoothing terms is solved to obtain the initial solution.

[0139] In the constraint confidence determination stage, the thermal infrared signal-to-noise ratio (SNR) of each grid is calculated, and weak signal grids with an SNR less than 3 are removed. The constraint confidence index is calculated, and grids with a constraint confidence index less than or equal to 1 are determined to be constraint-effective grids. In this embodiment, constraint-effective grids account for approximately 25% to 35% of the total number of grids.

[0140] In the physically constrained Gaussian inversion stage, the region is discretized into a 100m x 100m grid. The Gaussian plume transfer matrix is ​​calculated based on ERA5 meteorological data. An objective function with inequality constraints is constructed and solved using the effective set method.

[0141] During the iterative calibration and convergence phase, the inversion results are input into the Gaussian forward model to simulate column concentration, with an initial relative deviation of 9.6%. Thermal efficiency, emission factor, and radiative heat dissipation ratio are fine-tuned within physically reasonable limits. If parameter adjustments cause the lower bound of emission intensity for a certain grid to increase beyond the current iterative solution, the information for that grid is passed to the constraint conflict detection step. After four iterations, the deviation decreases to 4.3%, reaching the absolute convergence threshold.

[0142] In the constraint conflict detection phase, step S54 receives constraint conflict meshes caused by parameter changes. In this embodiment, approximately 1.5% of the constraint-effective meshes exhibit minor conflicts during the iteration process, with a constraint conflict degree of less than 0.05. Upon inspection, the conflicting meshes are mainly located at the edges of cloud-covered areas, and are determined to be unreliable thermal infrared data. They are removed from the set of physically constraint-effective meshes, and iteration continues until convergence is achieved.

[0143] In the comparative verification phase, a control group was set up to verify the effectiveness of the invention. The data was compared with the emission data of key polluting enterprises published by the local ecological and environmental authorities in 2020. The correlation coefficient of the nighttime light method at a 100-meter grid scale was 0.58, the normalized mean deviation was 35%, and there was no physical consistency. The correlation coefficient of the soft prior regularization method (i.e., the method without hard constraints) was 0.72, the normalized mean deviation was 18%, and the physical consistency was weak. The correlation coefficient of the method of the present invention was 0.84, the normalized mean deviation was 10%, and the physical consistency was strong. The results show that the method of the present invention has a significant improvement compared with the soft prior method without hard constraints.

[0144] Example 2: Taking a mixed industrial area with multiple industries as an example.

[0145] The main difference between this embodiment and Embodiment 1 is that the industry types involved in this embodiment are diverse, requiring the industry type to be identified grid by grid and the parameters to be dynamically determined.

[0146] During the data acquisition and preprocessing phase, the region's MEIC emission inventory, SDGSAT-1 TIS thermal infrared imagery, TROPOMI NO2 and OCO-2 or OCO-3 XCO2 products, as well as meteorological data and point-of-interest (POI) data are acquired.

[0147] In the multi-industry parameterization stage, industry type is determined grid-by-grid by combining point-of-interest data and thermal infrared temperature characteristics. For grids identified as steel smelting, a conservative lower bound is set for thermal efficiency (0.50-0.65) and a conservative upper bound is set for radiative heat dissipation (0.30-0.50). For grids identified as thermal power, a conservative lower bound is set for thermal efficiency (0.35-0.45) and a conservative upper bound is set for radiative heat dissipation (0.15-0.25). For grids identified as cement, a conservative lower bound is set for thermal efficiency (0.40-0.55) and a conservative upper bound is set for radiative heat dissipation (0.25-0.40). For grids identified as chemical industry, a conservative lower bound is set for thermal efficiency (0.45-0.60) and a conservative upper bound is set for radiative heat dissipation (0.15-0.30). The physical lower bound of emission intensity for each grid is dynamically determined.

[0148] In the multi-species joint inversion stage, the spatial distribution structure of emissions is constrained by the TROPOMI NO2 tropospheric column concentration product, and the total emission level is constrained by the OCO-2 or OCO-3 XCO2 product. A joint objective function is constructed and solved.

[0149] During the constraint anomaly source identification phase, approximately 1.8% of the constraint-effective grids were found to have constraint conflicts due to parameter changes, most of which were located in areas with a high concentration of chemical enterprises. The thermal infrared data of these grids were of good quality, and they were identified as candidate regions for non-thermal radiation-dominated emission sources, which is consistent with the characteristics of the chemical industry, which is mainly characterized by process emissions.

[0150] The final regional total deviation was 4.5%, and the correlation coefficient between the 100-meter grid scale and the independent survey data was 0.82.

[0151] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments, and various changes and modifications can be made without departing from the spirit and scope of the invention, all of which fall within the scope of protection claimed by the present invention. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A spatial downscaling method for industrial source emission inventories, characterized in that, Includes the following steps: Obtain a coarse-resolution industrial emission inventory of the target area, and acquire simultaneous thermal infrared remote sensing data and satellite atmospheric column concentration observation data; Based on the aforementioned thermal infrared remote sensing data, thermal infrared radiation flux is retrieved, and a physical lower bound constraint for emission intensity is constructed according to the principle of thermodynamic energy conservation. The physical lower bound constraint is constructed as follows: thermal efficiency, radiative heat dissipation ratio, and emission factor are obtained. The radiative heat dissipation ratio is the proportion of radiative heat dissipation in waste heat. The waste heat ratio is obtained by subtracting the thermal efficiency from 1. The product of the waste heat ratio and the radiative heat dissipation ratio is calculated as a first product. The thermal infrared radiation flux is divided by the first product to obtain an estimated total fuel combustion energy value. The estimated total fuel combustion energy value is multiplied by the emission factor to obtain a lower bound value for emission intensity, and the emission intensity is not less than the lower bound value for emission intensity. An inversion equation is constructed based on the analytical solution of the Gaussian plume model. The satellite atmospheric column concentration observation data is used as the input of the inversion equation, and the physical lower bound constraint is used as the inequality constraint to solve for the spatial distribution of emission intensity in the target area. Based on the spatial distribution of emission intensity obtained from the inversion solution, the coarse-resolution industrial source emission inventory is spatially downscaled to output a high-resolution industrial source emission intensity distribution.

2. The method according to claim 1, characterized in that, The step of retrieving thermal infrared radiant flux based on the thermal infrared remote sensing data specifically includes: Radiometric calibration and atmospheric correction were performed on the thermal infrared remote sensing data to obtain brightness and temperature data; According to Planck's law, the brightness temperature data is integrated within a preset thermal infrared band to obtain the surface radiation flux, and the surface radiation flux is determined as the thermal infrared radiation flux. Alternatively, when the thermal infrared remote sensing data is a VIIRS fire radiation power product, the radiation power value of the VIIRS fire radiation power product can be directly determined as the thermal infrared radiation flux.

3. The method according to claim 1, characterized in that, Before applying the inequality constraints to the inversion equations, the method further includes a step of determining the set of meshes where the physical constraints are effective, specifically including: The initial solution of the emission intensity is obtained by solving the inversion equation without applying the aforementioned inequality constraints. Calculate the thermal infrared signal-to-noise ratio for each grid and calculate the constraint confidence index, which is the ratio of the initial solution to the lower bound of the emission intensity. When the signal-to-noise ratio of the thermal infrared signal is greater than the preset signal-to-noise ratio threshold and the constraint confidence index is less than or equal to 1, the corresponding mesh is included in the set of meshes where the physical constraint is effective. In the step of using the physical lower bound constraint as an inequality constraint, the inequality constraint is applied only to the meshes in the set of meshes where the physical constraint is effective.

4. The method according to claim 3, characterized in that, The steps for inverting and solving the spatial distribution of emission intensity in the target region specifically include: The study area was discretized into high-resolution grid cells, and the Gaussian plume transfer matrix was calculated based on meteorological data. A regularization objective function with inequality constraints is constructed, wherein the regularization objective function includes a data fitting term, a spatial smoothing regularization term, and a thermal infrared prior regularization term; The emission intensity vector is obtained by solving the convex quadratic programming problem with inequality constraints using the effective set method or the interior point method.

5. The method according to claim 4, characterized in that, The inversion solution process also includes a parameter iterative calibration step: Input the spatial distribution of emission intensity obtained from the current iteration into the Gaussian plume forward model to simulate the satellite column concentration; Calculate the relative deviation between the simulated satellite column concentration and the observed satellite atmospheric column concentration; When the relative deviation is greater than the preset convergence threshold, the thermal efficiency, the emission factor and the proportion of radiative heat dissipation are adjusted within a physically reasonable range, and the lower limit of emission intensity is recalculated. The process returns to the step of performing the inversion solution of the spatial distribution of emission intensity in the target area and performs the next iteration solution. In a single iteration, the adjustment range of each parameter does not exceed the preset first adjustment threshold, and the cumulative adjustment range does not exceed the preset second adjustment threshold.

6. The method according to claim 5, characterized in that, The parameter iterative calibration process also includes steps for constraint conflict detection and anomaly source identification: Monitor whether, after parameter adjustment, a grid appears where the current iterative solution is less than the recalculated lower bound of the emission intensity; If it exists, then calculate the constraint conflict degree, which is the ratio of the difference between the lower bound of the emission intensity and the current iterative solution to the lower bound of the emission intensity. When the constraint conflict degree is greater than the preset conflict threshold, the corresponding mesh is marked as a constraint abnormal mesh; For the constrained abnormal grid, the data quality of its thermal infrared remote sensing data is detected. If the data quality is lower than the preset quality standard, the thermal infrared data is determined to be unreliable. If the data quality meets the preset quality standard, it is determined to be a candidate region of non-thermal radiation dominant emission source, and the grid is removed from the set of grids with effective physical constraints in subsequent iterations. Based on the constraint conflict degree, the data quality of the thermal infrared remote sensing data, the quantity and quality of satellite observations, and the iterative convergence status, the confidence level of each grid is comprehensively evaluated, and the confidence level is correlated with the spatial distribution of emission intensity and output; the confidence level includes at least high confidence, medium confidence and low confidence.

7. The method according to claim 1, characterized in that, The thermal efficiency, the proportion of radiative heat dissipation, and the emission factor are determined in the following ways: The industry type of the target grid is determined based on thermal infrared temperature characteristics, land use data, and / or point of interest data; Based on the industry type, the corresponding range of thermal efficiency, range of radiative heat dissipation ratio, and range of emission factor are searched from the preset parameter library; In this process, a conservative strategy is adopted when selecting parameters, with the radiative heat dissipation ratio taking the upper limit of its preset range and the thermal efficiency taking the lower limit of its preset range. During the parameter iterative calibration process, when it is necessary to adjust parameters first, the proportion of radiative heat dissipation should be adjusted first to keep the thermal efficiency and the emission factor relatively stable.

8. The method according to claim 7, characterized in that, The industry types mentioned include iron and steel smelting, thermal power, cement, and chemical industries; The preset parameter library stores the mapping relationship between the industry type and the range of thermal efficiency and the range of radiative heat dissipation ratio; the range of thermal efficiency and the range of radiative heat dissipation ratio are determined according to the typical process temperature and heat dissipation characteristics of the industry.

9. The method according to claim 1, characterized in that, The thermal infrared remote sensing data includes remote sensing image data in at least one thermal infrared band; The satellite atmospheric column concentration observation data includes TROPOMI NO2 tropospheric column concentration products and OCO-2 / 3 XCO2 products; The inversion solution adopts a multi-species joint inversion strategy, using the TROPOMI NO2 tropospheric column concentration product to constrain the spatial distribution structure of emissions, and using the OCO-2 / 3 XCO2 products to constrain the total emission level.

10. The method according to claim 5, characterized in that, The preset convergence threshold includes an absolute convergence threshold and a relative convergence threshold; When the relative deviation is less than the absolute convergence threshold, it is determined to be converged; Alternatively, convergence is determined when the relative deviation change rate between two consecutive iterations is less than the relative convergence threshold. If convergence is not achieved after reaching the preset maximum number of iterations, output the result with the smallest relative deviation in the history of iterations.

Images