A PM2.5 and O3 collaborative pollution precision tracing method based on multi-source heterogeneous data fusion

Abstract

This invention relates to the field of environmental monitoring technology and discloses a method for accurate source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion. The method includes: resampling multi-source heterogeneous atmospheric environmental data to obtain a multi-source fusion dataset; using this dataset as a constraint filter to update the emission inventory sample set to obtain a dynamic emission inventory; normalizing the regulatory intensity of meteorological factors in a specific meteorological process to obtain a dynamic meteorological constraint matrix; using this matrix to perform meteorological correction on the concentration field and emission inventory to obtain an emission characteristic concentration field and a net source inventory; constructing a spatiotemporal correlation matrix based on the two and performing convolution processing to obtain gridded spatiotemporal convolution response values; resolving the source of the response values ​​from the gridded points to the receiver to obtain the synergistic pollution source tracing results; and superimposing the source tracing results with the emission inventory to obtain the source tracing operational results. This invention can improve the efficiency of accurate source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring technology, and in particular to a method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion. Background Technology

[0002] Current mainstream air pollution source tracing technologies have significant limitations when facing combined PM2.5 and O3 pollution scenarios, failing to meet practical prevention and control needs in terms of accuracy and timeliness. Existing technologies largely rely on monitoring data from single ground stations, failing to integrate multi-source, three-dimensional monitoring methods such as UAV vertical detection, high-resolution component monitoring, and satellite remote sensing inversion. This results in the incomplete capture of the spatiotemporal distribution characteristics of key precursors, fundamentally incomplete data foundations for source tracing analysis. Furthermore, traditional source tracing models simply correlate conventional meteorological parameters, failing to establish differentiated coupling mechanisms for special meteorological scenarios such as stable weather, drought, and severe convection. This makes it impossible to distinguish between natural fluctuations in meteorological factors and the true contribution of anthropogenic emissions, leading to severe biases in the source location of emission sources in the source tracing results. In addition, existing technologies employ a static analysis model with fixed time periods, lacking responsiveness to dynamic processes such as intermittent emissions from industrial parks, seasonal releases from ecosystems, and cross-regional pollution transport. Source tracing results lag significantly behind changes in the pollution situation, essentially rendering them ineffective for the time-sensitive requirements of emergency control during heavy pollution weather.

[0003] Existing source tracing technologies also suffer from serious shortcomings in targeted analysis and operational implementation of synergistic pollution, failing to support the actual control needs of synergistic emission reduction for PM2.5 and O3. Most existing technologies focus on independent source tracing of single pollutants, neglecting the close coupling reaction mechanism between PM2.5 secondary formation and O3 photochemical pollution. They cannot accurately identify common precursors that simultaneously exert key influences on both types of pollution, nor can they quantify the respective contribution weights of different industries, natural sources, and anthropogenic sources in synergistic pollution. Therefore, synergistic emission reduction cannot be implemented due to the lack of a unified and accurate data foundation. More significantly, the outputs of traditional technologies are mostly theoretical analysis data or conclusions at the academic paper level, lacking a standardized, dynamically updated grid-based screening system. They also lack hotspot lists that can directly guide environmental enforcement personnel in conducting on-site control and case libraries adapted to different emission reduction scenarios. Many technological achievements remain at the research level because they cannot be transformed into operational control tools, creating a huge gap between them and the intelligent supervision and precise pollution control capabilities actually required by local governments and environmental protection departments. Summary of the Invention

[0004] This invention provides a precise source tracing method for PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, this invention provides a precise source tracing method for the synergistic pollution of PM2.5 and O3 based on multi-source heterogeneous data fusion, comprising: S1. Spatiotemporally resample the multi-source heterogeneous atmospheric environment data of the target process to obtain the multi-source fusion dataset of the target process; S2. Using the multi-source fusion dataset as an observation constraint, the emission inventory sample set of the target process is iteratively updated through integrated filtering to obtain the dynamic emission inventory of the target process; S3. Normalize the control intensity of the meteorological factors in the multi-source fusion dataset during a specific meteorological process to obtain the dynamic meteorological constraint matrix of the target process. S4. Based on the dynamic meteorological constraint matrix, perform meteorological correction on the PM2.5 and O3 pollutant concentration fields of the multi-source fusion dataset and the dynamic emission inventory to obtain the emission characteristic concentration field and net source inventory of the target process. S5. Construct a spatiotemporal correlation matrix based on the dynamic emission inventory and the emission characteristic concentration field, perform spatiotemporal convolution processing on the spatiotemporal correlation matrix to obtain the grid spatiotemporal convolution response value of the target process, and perform receptor source analysis on the source grid data corresponding to the grid spatiotemporal convolution response value to obtain the collaborative pollution source tracing result of the target process. S6. Overlay the collaborative pollution source tracing results with the dynamic emission inventory in a grid format to obtain the source tracing business results of the target process.

[0006] In a preferred embodiment, the step of performing spatiotemporal resampling on the multi-source heterogeneous atmospheric environment data of the target process to obtain the multi-source fused dataset of the target process includes: Acquire atmospheric multi-source monitoring data for the target process, including fixed-site monitoring data, mass spectrometer component monitoring data, UAV vertical profile data, meteorological observation data, satellite remote sensing inversion data, and enterprise online emission time series data; By statistical outlier detection, abnormal equipment values ​​and extreme meteorological interference values ​​are removed from the atmospheric multi-source monitoring data to obtain the cleaned multi-source data stream of the target process; Spatial interpolation is used to fill in the missing data points in the multi-source data stream being cleaned, thereby obtaining the gridded data field of the target process; Align the timestamps of the data layers in the gridded data field to the same time base, and project the spatial coordinates of the data layers to the same spatial reference system to obtain the multi-source fusion dataset of the target process.

[0007] In a preferred embodiment, the step of iteratively updating the emission inventory sample set of the target process using the multi-source fusion dataset as an observation constraint and through ensemble filtering to obtain the dynamic emission inventory of the target process includes: Obtain an initial emission inventory sample set for the target process, the initial emission inventory sample set containing emission inventory samples generated based on the regional pollution source ledger; The ground component monitoring data, UAV vertical profile data and satellite remote sensing data of the multi-source fusion dataset are combined into the observation constraint vector of the target process; The residual between the simulated concentration values ​​of the samples in the initial emission inventory sample set on a unified spatial grid and the observation constraint vector is calculated to obtain the deviation sequence of the target process; The Kalman gain matrix of the initial emission inventory sample set is calculated based on the deviation sequence, and the emission amount of the sample is weighted and corrected using the Kalman gain matrix to obtain the assimilation sample set of the target process. The average value of the assimilated sample set is calculated along the sample dimension to obtain the dynamic emission inventory of the target process.

[0008] In a preferred embodiment, the step of normalizing the regulation intensity of the meteorological factors in the multi-source fusion dataset during a specific meteorological process to obtain the dynamic meteorological constraint matrix of the target process includes: The wind field vector, temperature stratification parameters, atmospheric boundary layer height, turbulence intensity time series observation sequence, and precipitation intensity time series observation sequence of the multi-source fusion dataset are used as the meteorological factor time series of the target process. The meteorological process type is identified by performing meteorological process type identification on the time series of the meteorological factors to obtain the specific meteorological process type of the target process. The specific meteorological process type includes stable weather process, drought and high temperature process and severe convective weather process. Extract time-period subsequences corresponding to the specific meteorological process type from the time-series sequence of the meteorological factors, and use the correlation coefficients between the meteorological factors in the time-period subsequences and the concentrations of PM2.5 and O3 pollutants in the corresponding time period as the regulation intensity value of the target process; Based on the total control intensity of the sub-meteorological factors, the control intensity value is transformed to obtain the allocation coefficient of the target process, and the allocation coefficient is mapped to the grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process.

[0009] In a preferred embodiment, the step of converting the control intensity value based on the total control intensity of the sub-meteorological factors to obtain the allocation coefficient of the target process, and mapping the allocation coefficient to grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process, includes: Obtain the sequence of regulation intensity values ​​of the sub-meteorological factors in the specific meteorological process, and use the ratio of the regulation intensity value in the sequence of regulation intensity values ​​to the sum of the regulation intensity values ​​of the sub-meteorological factors as the allocation coefficient of the target process; The comprehensive meteorological constraint value of the target process is calculated based on the measured values ​​of the sub-meteorological factors and the allocation coefficients; The comprehensive meteorological constraint values ​​are combined according to the spatial arrangement and time step order of the standard grid to obtain the dynamic meteorological constraint matrix of the target process.

[0010] In a preferred embodiment, the formula for calculating the comprehensive meteorological constraint value is: in, The comprehensive meteorological constraint value is... For static and stable cumulative weighting coefficients, For atmospheric stability parameters, For temperature stratification parameters, For humidity parameters, Ventilation coefficient, The temperature and humidity saturation feedback coefficient is... For precipitation removal weighting coefficient, Given the current rainfall intensity, The attenuation coefficient for ventilation wet removal is... This is a reference value for the ventilation coefficient. The lagging precipitation coupling weighting coefficient, For lag Precipitation intensity The number of time steps is the lag time. This is the humidity suppression feedback coefficient.

[0011] In a preferred embodiment, the step of performing meteorological correction on the PM2.5 and O3 pollutant concentration fields of the multi-source fusion dataset and the dynamic emission inventory based on the dynamic meteorological constraint matrix to obtain the emission characteristic concentration field and net source inventory of the target process includes: Based on the multi-source fusion dataset, the concentration fields of PM2.5 and O3 pollutants on the baseline grid in the target process are extracted to obtain the original pollutant concentration field of the target process; The dynamic meteorological constraint matrix is ​​divided into transport constraint sub-matrices, diffusion constraint sub-matrices, accumulation constraint sub-matrices, and wet deposition constraint sub-matrices according to the mechanism of action of meteorological factors. Based on the transmission constraint submatrix, the diffusion constraint submatrix, the cumulative constraint submatrix, and the wet settlement constraint submatrix, the transmission contribution, diffusion contribution, cumulative contribution, and wet settlement reduction of the grid points on the reference grid are calculated respectively. The emission contribution concentration value of the target process is obtained by subtracting the weighted sum of the transport contribution, the diffusion contribution, the cumulative contribution, and the wet deposition reduction from the grid point concentration values ​​of the original pollutant concentration field. The emission contribution concentration value is then spatially reorganized to obtain the emission characteristic concentration field of the target process.

[0012] In a preferred embodiment, the step of constructing a spatiotemporal correlation matrix based on the dynamic emission inventory and the emission characteristic concentration field, performing spatiotemporal convolution processing on the spatiotemporal correlation matrix to obtain the grid-point spatiotemporal convolution response value of the target process, and performing receptor-source analysis on the source grid data corresponding to the grid-point spatiotemporal convolution response value to obtain the collaborative pollution source tracing result of the target process includes: The gridded emissions of pollutant species in the dynamic emission inventory are expanded along the spatial dimension to obtain the emission space vector of the target process. The gridded concentration values ​​of PM2.5 and O3 in the emission characteristic concentration field are expanded along the spatial dimension to obtain the concentration space vector of the target process. The emission space vector and the concentration space vector are stacked in time series to obtain the emission time series space matrix and the concentration time series space matrix of the target process. Calculate the spatial cross-correlation matrix between the emission time series spatial matrix and the concentration time series spatial matrix, and extend the spatial cross-correlation matrix along the time step direction in the time dimension to obtain the spatiotemporal correlation matrix of the target process; The neighborhood association strength of the grid points in the spatiotemporal correlation matrix is ​​causally convolved in the time step direction to obtain the grid spatiotemporal convolution response value of the target process. Based on the distribution of the spatiotemporal convolution response values ​​of the grid points, the grid points whose response values ​​exceed the target and a set verification threshold are identified as pollution source grid points. Based on the location information of the pollution source grid points, the concentration time series of the pollutant species are extracted from the emission characteristic concentration field. The concentration time series is subjected to non-negative matrix decomposition to obtain the source contribution matrix and source spectrum matrix of the target process. Based on the source contribution matrix and the source spectrum matrix, the key precursor categories and pollution contribution of the target process are determined. The PM2.5 to O3 concentration ratio, temperature parameters, and radiation parameters of the pollution source grid points are obtained. Based on the temperature parameters and radiation parameters, the control type threshold of the target process is determined. Based on the control type threshold, the concentration ratio is judged to obtain the sensitivity type of the target process. The location and response value intensity of the pollution source grid points, the category of the key precursor, the pollution contribution, and the sensitivity type are combined to form the collaborative pollution source tracing result of the target process.

[0013] In a preferred embodiment, the step of overlaying the collaborative pollution source tracing results with the dynamic emission inventory in a grid-like manner to obtain the source tracing business results of the target process includes: The response values ​​of the source tracing grid points in the collaborative pollution source tracing results are matched with the pollutant species emission amounts at the same grid point in the dynamic emission inventory to obtain the pollution source matching table for the target process; According to the pollution source matching table, the pollutant proportion set of the source tracing grid points is extracted from the dynamic emission inventory, and the precursor category corresponding to the pollutant proportion set is used as the key precursor identifier of the target process; The location information of the source tracing grid points is spatially overlaid with the emission source category information in the dynamic emission inventory to identify the emission source category to which the source tracing grid points belong. The proportion of the response value corresponding to the emission source category in the total response value of the source tracing grid is taken as the pollution contribution of the target process; The location information of the source tracing grid points, the identification of the key precursors, the emission source category, and the pollution contribution are combined to form the source tracing business results of the target process.

[0014] In a preferred embodiment, the step of using the proportion of the response value corresponding to the emission source category in the total response values ​​of the source tracing grid as the pollution contribution of the target process includes: The response values ​​of the source grid points are accumulated over time steps to obtain the cumulative response intensity of the target process; The grid points on the reference grid in the target process are sorted according to the cumulative response intensity to obtain the cumulative response intensity sequence of the target process, and the grid points in the cumulative response intensity sequence whose response intensity reaches the cumulative contribution threshold are used as the hot spot screening grid point set of the target process. Spatially match the location of the hotspot screening grid points with the emission source categories marked in the dynamic emission inventory to obtain the actual emission source identifier of the target process; Extract the component emissions corresponding to the actual emission source identification from the dynamic emission inventory, and calculate the contribution value of the hot spot screening grid points to the common precursors of PM2.5 secondary formation and O3 photochemical pollution based on the component emissions; The proportion of the contribution value corresponding to the emission source category in the contribution value of the common precursors is used as the pollution contribution of the target process.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention constructs an integrated air-ground multi-source three-dimensional monitoring data acquisition system, which integrates six types of heterogeneous data—fixed-site monitoring, high-resolution mass spectrometry component monitoring, UAV vertical profile scanning, refined meteorological observation, satellite remote sensing inversion, and enterprise online emissions—through spatiotemporal unified registration and standardized fusion. This transforms the data foundation for source tracing analysis from a single-dimensional ground concentration monitoring to a multi-dimensional three-dimensional information network covering the vertical range from the ground to the boundary layer and including chemical components, meteorological conditions, and emission dynamics. This complete data fusion architecture enables the comprehensive capture and accurate characterization of the spatiotemporal distribution characteristics of key precursors, providing complete data input for subsequent assimilation inversion and source tracing modeling. The localized ENKF ensemble Kalman filter assimilation and inversion framework built on this basis uses a multi-source fusion dataset as an observation constraint to iteratively optimize and update the initial emission inventory. This significantly improves both the spatial and temporal resolution of the emission inventory and enables continuous iteration to correct emission biases at fixed periods. This fundamentally ensures the accurate representation of the actual emission status of the region by the dynamic emission inventory. At the same time, by splitting the dynamic meteorological constraint matrix into four independent constraint sub-matrices based on the mechanism of meteorological factors—transmission, diffusion, accumulation, and deposition—and quantifying the degree of meteorological interference at each grid point, the emission contribution concentration values ​​obtained after removing meteorological contributions from the original concentration field truly reflect the actual contribution level of anthropogenic emission sources, completely eliminating the interference of meteorological fluctuations on the source tracing results.

[0016] 2. This invention also achieves significant technological advancements in collaborative pollution analysis and operational output. By constructing a spatiotemporal graph convolutional neural network source tracing model, it uses dynamic emission inventories and emission characteristic concentration fields as core inputs to build a spatiotemporal correlation matrix and performs causal convolution operations. This allows the source tracing results to simultaneously acquire spatial location information, temporal evolution patterns, and temporal causal relationships of pollution sources. This enables precise hourly and grid-by-grid capture of dynamic pollution scenarios such as intermittent emissions, seasonal releases, and cross-regional transmission in industrial parks. Based on this, by integrating the dual analysis mechanism of the PMF source apportionment model and the OBM box model, it can simultaneously identify common key precursors of PM2.5 secondary generation and O3 photochemical pollution and quantify the collaborative contribution share of different emission source categories to the two types of pollution. The key precursor identifiers output directly point to the specific control direction of VOCs. Ultimately, by constructing a standardized grid-based hotspot screening system and emission reduction case library, the source tracing results are transformed into complete source tracing business outcomes, including precise location of pollution sources, identification list of key precursors, quantified data on the dynamic contribution of multiple pollution sources, and priority control ranking list. These outcomes record in detail, by grid point, the type of emission source belonging to each key area, the types of pollutants that should be prioritized for control, the contribution ratio of this type of source in the overall total, and the corresponding specific control strategies. This has the capability to directly support local governments and environmental protection departments in carrying out intelligent supervision and precise pollution control in a business application. Furthermore, the data processing flow, assimilation and inversion method, source tracing modeling logic, and screening output system in the entire technical architecture can be fully reused in different cities and different industrial parks, demonstrating technical universality for cross-regional promotion. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating a method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion, provided in an embodiment of the present invention. The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0018] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0019] This application provides a method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion. The executing entity of this method includes, but is not limited to, at least one of the following electronic devices that can be configured to execute the method provided in this application: a server, a terminal, etc. In other words, the method can be executed by software or hardware installed on a terminal device or a server device. The server includes, but is not limited to, a single server, a server cluster, a cloud server, or a cloud server cluster. The server can be an independent server or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, content delivery networks (CDNs), and big data and artificial intelligence platforms.

[0020] Reference Figure 1 The diagram shown is a flowchart illustrating a method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion, according to an embodiment of the present invention. In this embodiment, the method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion includes: In this embodiment of the invention, when performing spatiotemporal resampling on the multi-source heterogeneous atmospheric environment data of the target process to obtain the multi-source fusion dataset of the target process, it is specifically used for: Acquire atmospheric multi-source monitoring data for the target process, including fixed-site monitoring data, mass spectrometer component monitoring data, UAV vertical profile data, meteorological observation data, satellite remote sensing inversion data, and enterprise online emission time series data; By statistical outlier detection, abnormal equipment values ​​and extreme meteorological interference values ​​are removed from the atmospheric multi-source monitoring data to obtain the cleaned multi-source data stream of the target process; Spatial interpolation is used to fill in the missing data points in the multi-source data stream being cleaned, thereby obtaining the gridded data field of the target process; Align the timestamps of the data layers in the gridded data field to the same time base, and project the spatial coordinates of the data layers to the same spatial reference system to obtain the multi-source fusion dataset of the target process.

[0021] Specifically, fixed-site monitoring data are obtained by continuous monitoring at a fixed sampling frequency by national-level ambient air quality automatic monitoring stations and atmospheric superstations distributed in the target process area. Mass spectrometer component monitoring data are obtained by continuous sampling at high temporal resolution by high-resolution online mass spectrometers deployed at key locations in the target process area. UAV vertical profile data are obtained by UAVs equipped with atmospheric component sensors performing vertical profile scanning in key parks and ecological functional zones in the target process area according to preset routes.

[0022] Specifically, statistical outlier detection is performed independently for each type of data in the multi-source atmospheric monitoring data. For each type of data, all measurement values ​​of each monitoring point within a continuous time window are first obtained. These measurement values ​​are arranged in chronological order to form an observation sequence. Then, the central tendency and dispersion of all measurement values ​​in the observation sequence are calculated. Based on this, the normal fluctuation range of the data type is determined. Measurement values ​​outside the normal fluctuation range in the observation sequence are marked as outliers. For measurement values ​​marked as outliers, the operating status identifier of the monitoring equipment corresponding to the measurement value and the meteorological observation record at the corresponding time are obtained simultaneously.

[0023] Specifically, spatial interpolation is performed separately for each type of data in the cleaning multi-source data stream. For each type of data, the effective measurement values ​​and corresponding spatial coordinates of all monitoring points in the target process area are first obtained. These monitoring points are discretely distributed in the target process area and the point density is uneven. The coverage of different monitoring points is different. Then, a regular grid system covering the entire area is constructed on the complete spatial range of the target process area. The regular grid system is composed of grid points formed by the intersection of crisscrossing grid lines, and adjacent grid points are distributed at equal intervals in space.

[0024] Specifically, the time base and timestamp information of each record in each data layer of the gridded data field are first obtained. Since data layers from different sources use their own independent time recording methods during the acquisition phase, the timestamps of different data layers may have deviations. Therefore, a unified target time base is determined from the timestamps of all data layers. Then, the actual timestamp of each record in each data layer is extracted, and the time offset between the actual timestamp and the target time base is calculated. Based on the time offset, the actual timestamp of each record is adjusted to the nearest time point under the target time base, so that after adjustment, records in all data layers at the same time point correspond to the same physical time.

[0025] Furthermore, meteorological observation data is obtained by continuous observation at a fixed sampling frequency by ground meteorological stations and boundary layer meteorological sounding equipment in the target process area. Satellite remote sensing inversion data is obtained by multispectral inversion of the target process area by atmospheric sounders carried by polar-orbiting and geostationary satellites. Enterprise online emission time series data is collected and uploaded in real time by the continuous emission monitoring system installed by key polluting enterprises in the target process area. After preliminary verification in their respective acquisition systems, the above six types of data are aggregated in real time to the target process data management platform in accordance with a unified data transmission protocol.

[0026] Furthermore, if the monitoring equipment operation status indicator corresponding to the measured value shows that the equipment has exceeded the calibration period, the sensor response is abnormal, or the communication is interrupted at the time of measurement, then the measured value is judged as an abnormal value of the equipment and is removed. If the meteorological observation record corresponding to the measured value shows that there is a sudden change in wind speed, a sudden change in temperature, or a drastic fluctuation in air pressure that exceeds the normal range at that time, then the measured value is judged as a meteorological extreme interference value and is removed. If the monitoring equipment operation status and meteorological observation record corresponding to the measured value are both normal, then the measured value is retained as a normal value. After outlier detection and outlier removal by category and point by point, the remaining effective measured values ​​of each monitoring point at each time are formed into a cleaned multi-source data stream according to their source data type and collection time.

[0027] Furthermore, for each grid point in the regular grid system, it is determined whether a monitoring point exists at that grid point location. If a monitoring point exists at that grid point location and a valid measurement value exists for that monitoring point in the cleaned multi-source data stream, then that valid measurement value is directly used as the value of that grid point. If no monitoring point exists at that grid point location or a valid measurement value is missing for that monitoring point in the cleaned multi-source data stream, then all known monitoring points within a preset radius around that grid point are searched, and the valid measurement values ​​at these known monitoring points are extracted. The spatial distance between each known monitoring point and the grid point is calculated, and different participation weights are assigned to the valid measurement values ​​of each known monitoring point based on the spatial distance. The closer the distance, the higher the participation weight. The higher the participation weight of the measurement point, the lower the participation weight of the monitoring point that is farther away. The effective measurement values ​​of all known monitoring points are weighted and summed according to their corresponding participation weights. The result of the weighted sum is used as the interpolation fill value for that grid point. After performing the above judgment and filling operation on each grid point in the regular grid system, each type of data source obtains the corresponding value on each grid point of the regular grid system. Thus, each type of data source is transformed from a set of measurement values ​​of discrete points into a set of continuous grid point values ​​covering the entire area space, and the grid point data layer corresponding to each type of data source is obtained. The grid point data layers corresponding to all data sources are superimposed according to the grid point position to form a grid point data field.

[0028] Furthermore, for records that still have missing timestamps after adjustment, time interpolation is performed based on the values ​​of adjacent times in the time series of the data layer containing that record. This ensures that the timestamps of each data layer in the gridded data field are all located under the same time reference. The spatial coordinate projection operation first obtains the spatial coordinate reference system used by each data layer in the gridded data field and the spatial coordinate values ​​of each grid point. Since data layers from different sources use their own independent spatial coordinate reference systems, the spatial positions of grid points in different data layers cannot be directly correlated. Therefore, a unified target spatial reference system is determined from the spatial coordinate reference systems of all data layers. Subsequently, the original spatial coordinates of each grid point in each data layer are extracted. Based on the projection transformation parameters between the source and target spatial reference frames where the original spatial coordinates are located, a coordinate mapping relationship is established from the source spatial reference frame to the target spatial reference frame. According to this coordinate mapping relationship, the original spatial coordinates of each grid point are converted into new spatial coordinates under the target spatial reference frame, so that records located at the same grid point coordinates in all data layers correspond to the same physical spatial location. After timestamp alignment and spatial coordinate projection, all data layers in the grid data field adopt a unified reference benchmark in both time and space dimensions. The values ​​of different data layers at the same spatial location at the same time can be directly jointly processed to obtain a multi-source fusion dataset.

[0029] In summary, by collecting six types of data from different sources and with different properties—fixed-site monitoring data, mass spectrometer component monitoring data, UAV vertical profile data, meteorological observation data, satellite remote sensing inversion data, and enterprise online emission time series data—this approach overcomes the limitations of traditional technologies that rely on single-site monitoring data. It enables subsequent source tracing analysis to simultaneously acquire near-ground concentration information, vertical component distribution information, large-scale regional pollution distribution information, and real-time emission dynamics information. This addresses the pain point of single monitoring dimensions in collaborative pollution source tracing from the data source itself, providing a complete data foundation for accurately capturing the spatiotemporal distribution characteristics of key precursors such as volatile organic compounds and nitrogen oxides. It avoids the problems of incomplete source tracing and missed identification of key pollution sources caused by incomplete data sources.

[0030] In summary, outlier detection was performed on six different types of raw monitoring data, each with distinct sources and characteristics. Normal fluctuation ranges were determined based on the inherent fluctuation characteristics of each data type. Measurements outside these normal fluctuation ranges were comprehensively evaluated in conjunction with their corresponding equipment operating status indicators and meteorological observation records. This effectively eliminated outlier values ​​caused by equipment calibration delays, abnormal sensor responses, or communication interruptions, as well as extreme meteorological interference values ​​caused by extreme weather conditions such as sudden wind speed changes, rapid temperature fluctuations, or drastic air pressure fluctuations. This ensured that the retained effective measurements more accurately reflected the actual emission characteristics of the emission sources, eliminating systematic errors introduced by outliers in subsequent assimilation and source tracing modeling from the source of data quality, and significantly improving source tracing accuracy.

[0031] In summary, considering the discrete distribution and uneven density of the six data sources within the target area, each data source is transformed from a set of discrete measurement points into a continuous set of grid point values ​​covering the entire area through spatial interpolation. This ensures that each data source obtains corresponding values ​​at every grid point in the regular grid system, resolving the issues of uneven spatial distribution of the original monitoring data and missing data in some areas. This guarantees that all grid points have complete data input in subsequent assimilation, inversion, and source tracing modeling, avoiding blind spots and result biases caused by missing data.

[0032] In summary, by adjusting the timestamps of all data layers in the gridded data field to a unified target time base, the six data layers, which originally had discrepancies due to different time recording methods of their respective acquisition systems, achieved accurate correspondence in the time dimension. By projecting the spatial coordinates of all data layers from their respective independent source space reference systems to a unified target space reference system, the six data layers, which originally could not be directly superimposed due to different coordinate projection methods, achieved precise registration in the spatial dimension. After the unification processing in both the temporal and spatial dimensions, the values ​​of different data layers at the same time and spatial location can be directly jointly processed. This provides high-quality input data with comparability and consistency in both the time and spatial dimensions for subsequent assimilation and inversion and meteorological constraint modeling using multi-source fusion datasets as input.

[0033] In this embodiment of the invention, when using the multi-source fusion dataset as an observation constraint and iteratively updating the emission inventory sample set of the target process through ensemble filtering to obtain the dynamic emission inventory of the target process, the specific method is as follows: Obtain an initial emission inventory sample set for the target process, the initial emission inventory sample set containing emission inventory samples generated based on the regional pollution source ledger; The ground component monitoring data, UAV vertical profile data and satellite remote sensing data of the multi-source fusion dataset are combined into the observation constraint vector of the target process; The residual between the simulated concentration values ​​of the samples in the initial emission inventory sample set on a unified spatial grid and the observation constraint vector is calculated to obtain the deviation sequence of the target process; The Kalman gain matrix of the initial emission inventory sample set is calculated based on the deviation sequence, and the emission amount of the sample is weighted and corrected using the Kalman gain matrix to obtain the assimilation sample set of the target process. The average value of the assimilated sample set is calculated along the sample dimension to obtain the dynamic emission inventory of the target process.

[0034] Specifically, when obtaining the initial emission inventory sample set of the target process, the latest regional pollution source ledger is first obtained from the ecological and environmental management department of the region where the target process is located. This ledger records the production process information, fuel consumption, raw and auxiliary material usage, pollution control facility operation parameters, and exhaust gas emission parameters of each pollution discharge link of all industrial enterprises in the region. At the same time, the ledger also includes the regional motor vehicle traffic statistics, the implementation status of road dust control measures, and the vegetation type distribution and area data of ecological natural sources.

[0035] Specifically, a gridded data layer corresponding to the ground component monitoring data is extracted from the multi-source fusion dataset. This data layer contains the concentration values ​​of various PM2.5 chemical components and various VOCs components at each grid point in the target process area. A gridded data layer corresponding to the UAV vertical profile data is extracted from the multi-source fusion dataset. This data layer contains the atmospheric component concentration values ​​at different vertical heights at each grid point in the target process area. A gridded data layer corresponding to the satellite remote sensing data is extracted from the multi-source fusion dataset. This data layer contains the aerosol optical thickness and trace gas column concentration values ​​at each grid point in the target process area.

[0036] Specifically, independent diffusion simulation calculations are performed for each emission inventory sample in the initial emission inventory sample set. For the emission inventory sample currently being processed, the emission data of each pollutant species at each grid point in the sample is input into the atmospheric chemical transport simulation system. The simulation system uses the emission data as the initial input and combines it with the topographic elevation data, land use type data, and boundary layer meteorological field data of the target process area to numerically simulate the advection transport process, turbulent diffusion process, chemical transformation process, and dry and wet deposition process of atmospheric pollutants on a unified spatial grid. After hourly integration calculations over all time steps, the simulation system outputs the simulated concentration values ​​of various pollutants at each grid point and at each time of the emission inventory sample on the unified spatial grid. The simulated concentration values ​​are organized according to the same spatial location, temporal order, and arrangement as the observation constraint vector to form the simulated concentration vector of the emission inventory sample. The simulated concentration vector corresponds completely to the observation constraint vector in terms of dimensional structure and arrangement order.

[0037] Specifically, the deviation sequences corresponding to all emission inventory samples are arranged in the order of the samples to form a set of deviation sequences. Based on this, the covariance between every two deviation sequences in the set of deviation sequences is calculated. The magnitude of the covariance reflects the degree of synchronization of the deviation of different samples at the same observation location with the sample. The larger the covariance, the more likely that the deviations of the two samples at the same location have the same trend of change. At the same time, the variance of each deviation sequence is calculated, which reflects the degree of dispersion of the deviation of the sample at different observation locations.

[0038] Specifically, for each pollutant species at each grid point in the assimilation sample set, the emission values ​​of all samples at that grid point for that species are extracted, and these values ​​are added together and divided by the total number of samples to obtain the average emission of that species at that grid point.

[0039] Furthermore, the parameter information of various emission sources recorded in the pollution source ledger is input into the emission factor calculation model. This model is configured with corresponding emission factors for each emission source category and each pollutant species. The emission source activity level data is multiplied by the corresponding emission factor to obtain the emission amount of the emission source per unit time. After calculating all emission sources category by category and species by species, the emission amount data of all emission sources in the region for all pollutant species are obtained. On this basis, multiple sets of random sampling are performed on the emission amount data by introducing uncertainty perturbation. Each set of sampling results constitutes an emission inventory sample. There are differences in the emission amount of the same species at the same grid point between the samples. These differences reflect the uncertainty in emission amount estimation caused by emission source ledger recording errors, uncertainty in emission factor selection, and statistical bias of activity level. All the emission inventory samples generated by sampling are summarized into an initial emission inventory sample set. Each sample in this sample set contains the emission amount data of each pollutant species at each grid point in the region where the target process is located.

[0040] Furthermore, the values ​​at the same grid point and at the same time in the above data layer are spliced ​​together according to a preset arrangement order, so that each grid point corresponds to a data vector at each time point, which is composed of ground component concentration value, vertical profile concentration value and satellite remote sensing inversion value. The data vectors of all grid points at all times are organized as a whole according to the spatial location and temporal order of the grid points to form an observation constraint vector. Each element in the observation constraint vector corresponds to a specific value of a specific grid point, a specific time and a specific observation type. The observation constraint vector as a whole represents the pollutant concentration observation information of the target process area at all grid points and all times provided by multi-source three-dimensional monitoring methods.

[0041] Furthermore, the difference between each observed value in the observation constraint vector and the simulated value of the same species at the same location and time in the simulated concentration vector is calculated to obtain the deviation between the simulated value and the observed value at that location. The deviations at all locations are arranged in the original order to form the deviation sequence corresponding to the emission inventory sample. The above diffusion simulation and deviation calculation operations are performed on each emission inventory sample in the initial emission inventory sample set to obtain the deviation sequence corresponding to each sample.

[0042] Furthermore, the covariances among all deviation sequences are organized into elements of a Kalman gain matrix according to the sample index and observation location index. The dimension of this Kalman gain matrix is ​​determined by the total number of samples and the total number of observation locations. The value of each element in the matrix represents the degree of common influence of the deviation of a certain sample at a certain observation location on the correction of emissions of all samples. Subsequently, for each emission inventory sample in the initial emission inventory sample set, the emission of each pollutant species of the sample at each grid point is weighted and corrected using the row vector corresponding to the sample in the Kalman gain matrix. The specific correction method is to subtract the deviation value corresponding to the sample in the deviation sequence of the corresponding species at each grid point from the current emission of the sample at each grid point by multiplying it by the product of the corresponding elements in the Kalman gain matrix. After calculating for each species at all grid points, the corrected emission data of the sample is obtained. The above weighted correction operation is performed on each sample in the initial emission inventory sample set to obtain the corrected emission data corresponding to each sample. All corrected samples are summarized into an assimilated sample set.

[0043] Furthermore, the above average calculation is performed on each grid point and each species of pollutant covered by the assimilation sample set to obtain the average emission data of each species at each grid point. The average emissions of all species at all grid points are organized according to the grid point arrangement order and species arrangement order of the unified spatial grid to form a dynamic emission inventory.

[0044] In summary, by inputting various emission source parameter information recorded in the regional pollution source ledger into the emission factor calculation model, and introducing uncertainty perturbations into the calculation results for multiple random sampling, multiple emission inventory samples are generated. This results in differences in the emission amounts of the same species at the same grid point among the samples in the initial emission inventory sample set. These differences cover the uncertainty in emission estimation caused by emission source ledger recording errors, uncertainty in emission factor selection, and statistical bias in activity levels. This provides an initial state set containing uncertainty information for subsequent iterative optimization of emissions under observation constraints through assimilation filtering. It fundamentally solves the problem of systematic deviation between traditional fixed emission inventories and actual regional emissions, making the final dynamic emission inventory more reflective of the true emission situation.

[0045] In summary, by integrating three types of three-dimensional monitoring data from different sources and with different observation dimensions—ground component monitoring data representing near-surface refined component concentration information, UAV vertical profile data representing atmospheric component concentration distribution information at different vertical altitudes, and satellite remote sensing data representing large-scale aerosol and trace gas column concentration information—into a unified observation constraint vector, this approach enables the observation constraint vector to simultaneously possess triple constraint capabilities: high-precision near-surface constraint, vertical structure constraint, and full regional coverage constraint. Compared to traditional technologies that rely solely on monitoring data from a single ground station as a constraint, this significantly enhances the ability of observation constraints to correct emission inventories.

[0046] In summary, by performing independent diffusion simulations on each emission inventory sample in the initial emission inventory sample set, the simulated concentration distribution of each sample at all grid points on a unified spatial grid is obtained. The simulated concentration values ​​are then compared point by point with the observed values ​​of the same species at the same location and time in the observation constraint vector to obtain the deviation sequence corresponding to each sample. This deviation sequence quantitatively characterizes the degree and pattern of deviation between the simulated concentration and the actual observed concentration of each emission inventory sample, providing quantitative input of deviation information for the subsequent calculation of the Kalman gain matrix. This allows the deviation differences between different samples to be accurately quantified and used to guide the direction and magnitude of emission correction.

[0047] In summary, by calculating the covariance and individual variances among all deviation sequences and organizing them into a Kalman gain matrix, each element in this gain matrix accurately reflects the degree of common influence of the deviation of a specific sample at a specific observation location on the correction of emissions of all samples. When using this gain matrix to perform weighted correction of emissions of each species at each grid point of each sample, the correction amount of each sample is jointly determined by its own deviation, the deviations of other samples, and the covariance between deviations. This achieves joint transmission and optimal fusion of multi-source observation information among all samples, ensuring that each sample in the assimilated sample set absorbs the information in the observation constraint vector. Compared with the traditional approach of not updating the emission inventory or only making simple adjustments, this significantly improves the accuracy of the emission inventory in representing the actual emission situation.

[0048] In summary, the average emission of a species at the same grid point is obtained by summing the emission values ​​of all samples in the assimilation sample set and dividing by the total number of samples. This average value reduces the residual random errors and uncertainties in each sample, while retaining the deterministic emission information corrected by multi-source observation constraints. At the same time, the dynamic emission inventory obtained by performing averaging calculations on all grid points and all species has both high spatiotemporal resolution and reduces the risk of extreme bias that may exist in a single sample through ensemble averaging. Compared with the static model of traditional fixed inventory, this dynamic emission inventory can continuously approach the actual emission situation as the observation data is continuously input and the iteration update cycle progresses.

[0049] In this embodiment of the invention, the step of normalizing the regulation intensity of the meteorological factors in the multi-source fusion dataset during a specific meteorological process to obtain the dynamic meteorological constraint matrix of the target process is specifically used for: The wind field vector, temperature stratification parameters, atmospheric boundary layer height, turbulence intensity time series observation sequence, and precipitation intensity time series observation sequence of the multi-source fusion dataset are used as the meteorological factor time series of the target process. The meteorological process type is identified by performing meteorological process type identification on the time series of the meteorological factors to obtain the specific meteorological process type of the target process. The specific meteorological process type includes stable weather process, drought and high temperature process and severe convective weather process. Extract time-period subsequences corresponding to the specific meteorological process type from the time-series sequence of the meteorological factors, and use the correlation coefficients between the meteorological factors in the time-period subsequences and the concentrations of PM2.5 and O3 pollutants in the corresponding time period as the regulation intensity value of the target process; Based on the total control intensity of the sub-meteorological factors, the control intensity value is transformed to obtain the allocation coefficient of the target process, and the allocation coefficient is mapped to the grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process.

[0050] Specifically, based on the data type labels carried by each type of data in the multi-source fusion dataset, the records labeled as wind field vectors are arranged in chronological order to form a wind field vector time series observation sequence; the records labeled as temperature stratification parameters are arranged in chronological order to form a temperature stratification parameter time series observation sequence; the records labeled as atmospheric boundary layer height are arranged in chronological order to form an atmospheric boundary layer height time series observation sequence; the records labeled as turbulence intensity are arranged in chronological order to form a turbulence intensity time series observation sequence; and the records labeled as precipitation intensity are arranged in chronological order to form a precipitation intensity time series observation sequence. Each record in the above five time series observation sequences contains the specific value of the meteorological factor at a specific grid point at a specific time, as well as the identification information of that time and that grid point.

[0051] Specifically, the meteorological factor values ​​at each time step in the meteorological factor time series are combined to form the meteorological state feature vector for that time step. The meteorological state feature vector of that time step is matched with a predefined static stable weather process feature template. If the ventilation coefficient of that time step is lower than a set range, the atmospheric stability parameter is higher than a set range, and the precipitation intensity is zero or close to zero, then that time step is determined to be a static stable weather process. The time step is then matched with a predefined drought and high temperature process feature template. If the precipitation intensity of that time step is lower than a set range, the temperature stratification parameter is higher than a set range, and the humidity parameter is lower than a set range, then that time step is determined to be a drought and high temperature process.

[0052] Specifically, the target type to be extracted is determined based on the specific meteorological process type. Then, the meteorological process type label of each time step in the meteorological factor time series is traversed. When the first time step labeled as the target type is encountered, this time step is taken as the starting position of the current time period. The traversal continues until the first time step not labeled as the target type is encountered. The time step before this time step is taken as the ending position of the current time period. All continuous time steps from the starting position to the ending position are taken as a complete time period subsequence. Then, the traversal continues from the next time step after the ending position. The next continuous time period labeled as the target type is found in the same way until all time steps are traversed, and all time period subsequences belonging to the target type are obtained. Then, each extracted time period subsequence is processed separately.

[0053] Specifically, the regulation intensity values ​​of all five sub-meteorological factors in the current specific meteorological process are obtained. These five regulation intensity values ​​are added together to obtain the total regulation intensity. Then, for each sub-meteorological factor, the regulation intensity value of the sub-meteorological factor is divided by the total regulation intensity to obtain the allocation coefficient of the sub-meteorological factor in the specific meteorological process. The value of the allocation coefficient is between zero and one, and the sum of the allocation coefficients of the five sub-meteorological factors is equal to one. The larger the allocation coefficient, the stronger the regulation effect of the meteorological factor on the change of pollutant concentration in the current specific meteorological process.

[0054] Furthermore, since all five types of data originate from multi-source fusion datasets and have undergone spatiotemporal registration processing, the five time-series observation sequences are consistent in terms of time reference and spatial reference frame. That is, the same grid point at the same time has a corresponding value in all five time-series observation sequences. Subsequently, these five time-series observation sequences are arranged and combined in a fixed order of wind field vector, temperature stratification parameter, atmospheric boundary layer height, turbulence intensity, and precipitation intensity to form an overall meteorological factor time-series sequence. Each time step of this meteorological factor time-series sequence contains the values ​​of all five meteorological factors at all grid points, providing complete meteorological information input for subsequent meteorological process type identification.

[0055] Furthermore, the time step is matched with a predefined feature template for severe convective weather processes. If the ventilation coefficient, atmospheric boundary layer height, and precipitation intensity of the time step are all above a set range, then the time step is determined to be a severe convective weather process. For each time step, the matching degree between the meteorological state feature vector of that time step and the three predefined feature templates is calculated. The meteorological process type corresponding to the feature template with the highest matching degree is taken as the preliminary identification result of that time step. Then, the temporal coherence of the preliminary identification results of all time steps is checked. If there is a sudden change in the preliminary identification result of a certain time step and the identification results of its immediate and adjacent time steps, the identification result of that time step is corrected according to the identification results of its immediate and adjacent time steps. Finally, the identification results of all time steps after the coherence check are summarized into the specific meteorological process type of the target process.

[0056] Furthermore, for each sub-meteorological factor in the currently processed time-segment subsequence, the values ​​of that sub-meteorological factor at each time step within the time-segment subsequence are used to form the time-segment numerical sequence of that factor. Simultaneously, the concentration values ​​of PM2.5 and O3 at each time step within the same time-segment subsequence are extracted from the multi-source fusion dataset to form the concentration time-segment numerical sequence. Then, the correlation coefficient between the time-segment numerical sequence of that sub-meteorological factor and the concentration time-segment numerical sequence is calculated. Specifically, the mean of the sequence is subtracted from each value in the time-segment numerical sequence to obtain the centered numerical sequence of that sequence, and the mean of the sequence is subtracted from each value in the concentration time-segment numerical sequence to obtain the centered coefficient of the concentration sequence. The numerator is the sum of the product of two centered values ​​at the same time step in the two centered numerical sequences. The denominator is the sum of the squares of each centered value in the centered numerical sequence of the meteorological factor. The denominator is the sum of the squares of each centered value in the centered numerical sequence of concentration. The denominator is the concentration value. The correlation coefficient of the meteorological factor in the time period is obtained by dividing the sum of the numerator by the product of the denominator factor value and the denominator concentration value. The sign and magnitude of the correlation coefficient reflect the same or opposite change relationship between the meteorological factor and the pollutant concentration and the degree of closeness. The correlation coefficient is used as the control intensity value of the meteorological factor in a specific meteorological process.

[0057] Furthermore, the measured value of each sub-meteorological factor at each grid point is multiplied by the allocation coefficient corresponding to that sub-meteorological factor to obtain the weighted contribution value of that sub-meteorological factor at that grid point. The weighted contribution values ​​of all five sub-meteorological factors at the same grid point are summed to obtain the comprehensive meteorological constraint value of that grid point. The comprehensive meteorological constraint value is calculated grid point by grid point according to the grid point arrangement order on the standard grid to obtain the comprehensive meteorological constraint value of each grid point at each time step. The comprehensive meteorological constraint values ​​of all grid points at all time steps are organized according to the spatial arrangement order and time step order of the standard grid to form a dynamic meteorological constraint matrix. Each element in this matrix corresponds to the comprehensive meteorological constraint value of a specific grid point at a specific time step. The matrix as a whole represents the distribution of meteorological constraint intensity formed by the combined action of multiple meteorological factors on all grid points and all time steps of the target process.

[0058] In summary, by selecting five categories of meteorological factors closely related to pollutant diffusion, transformation, and removal from multi-source fusion datasets and arranging them according to a unified time and spatial reference system, the horizontal transport conditions represented by wind field vectors, the atmospheric vertical structure characteristics represented by temperature stratification parameters, the vertical mixing spatial range represented by atmospheric boundary layer height, the vertical diffusion activity represented by turbulence intensity, and the wet removal capacity represented by precipitation intensity are fully aligned in both time and spatial dimensions. This provides a complete meteorological information input with a completely consistent time and space reference for subsequent meteorological process type identification and regulation intensity quantification, avoiding misjudgment of process types and deviations in regulation intensity calculation caused by inconsistent time references or spatial mismatches of various meteorological data. From the perspective of meteorological data input, this lays a data foundation for constructing a high-precision meteorological-pollution coupling constraint system.

[0059] In summary, by combining all meteorological factor values ​​at each time step into a meteorological state feature vector and matching it with feature templates for stable weather processes, drought and high temperature processes, and severe convective weather processes, the process type corresponding to the template with the highest matching degree is used as the preliminary identification result for that time step. After time series coherence verification and correction, the final identification result is obtained. This ensures that each time step of the target process is accurately classified into one of the following: stable weather process, drought and high temperature process, or severe convective weather process. This solves the problem that traditional source tracing models simply associate conventional meteorological parameters without distinguishing types for special meteorological scenarios. It provides an accurate process type classification basis for quantifying the intensity of meteorological factor regulation under different meteorological process types, avoiding the problems of distorted regulation intensity calculation and amplified source tracing bias caused by misjudgment of meteorological process types.

[0060] In summary, by extracting all continuous time steps belonging to the same specific meteorological process type from the time series of meteorological factors that have already been labeled with meteorological process types, a time-segment subsequence is constructed. This ensures that the meteorological conditions within each time-segment subsequence are homogeneous and consistent. Then, for each time-segment subsequence, the correlation coefficient between the numerical sequence of each type of meteorological factor within that time period and the numerical sequences of PM2.5 and O3 concentrations within the same time period is calculated. This correlation coefficient, obtained through hourly multiplication, accumulation, and normalization of the two centered sequences, accurately reflects the positive or negative driving strength of the meteorological factor on pollutant concentration changes under that specific meteorological process type. Compared with the traditional method of using uniform meteorological parameters for correlation throughout the entire time period, this significantly improves the relevance and accuracy of the quantitative results of meteorological factor regulation intensity for specific meteorological scenarios.

[0061] In summary, the total control intensity is obtained by summing the control intensity values ​​of the five meteorological factors in a specific meteorological process, and the allocation coefficient is obtained by dividing the control intensity value of each factor by the total control intensity. The sum of the five allocation coefficients is equal to one, and the value of each allocation coefficient accurately reflects the relative control contribution of the corresponding meteorological factor to the change of pollutant concentration among all meteorological factors in the specific meteorological process. Then, the measured value of each meteorological factor at each grid point is multiplied by the corresponding allocation coefficient and summed to obtain the comprehensive meteorological constraint value of that grid point. The comprehensive meteorological constraint value adapts to the spatial distribution differences of the measured values ​​of each factor among the grid points. After organizing the comprehensive meteorological constraint values ​​of all grid points and all time steps into a dynamic meteorological constraint matrix according to the spatial arrangement and time step order of the standard grid, this matrix realizes hourly and grid-point dynamic quantitative constraints on the impact of meteorological factors on the generation, accumulation, transmission and removal of pollutants under different specific meteorological scenarios such as static stability, drought and high temperature, and strong convection. This provides a dynamic constraint benchmark that can automatically adapt to changes in meteorological conditions for subsequent meteorological correction steps.

[0062] In this embodiment of the invention, when converting the control intensity value based on the total control intensity of the sub-meteorological factors to obtain the allocation coefficient of the target process, and mapping the allocation coefficient to grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process, the specific usage is as follows: Obtain the sequence of regulation intensity values ​​of the sub-meteorological factors in the specific meteorological process, and use the ratio of the regulation intensity value in the sequence of regulation intensity values ​​to the sum of the regulation intensity values ​​of the sub-meteorological factors as the allocation coefficient of the target process; The comprehensive meteorological constraint value of the target process is calculated based on the measured values ​​of the sub-meteorological factors and the allocation coefficients; The comprehensive meteorological constraint values ​​are combined according to the spatial arrangement and time step order of the standard grid to obtain the dynamic meteorological constraint matrix of the target process.

[0063] Specifically, when obtaining the sequence of regulation intensity values ​​of the sub-meteorological factors in the specific meteorological process, the regulation intensity values ​​of the five sub-meteorological factors in the specific meteorological process have been obtained from the previous steps. These five regulation intensity values ​​are arranged in a fixed order of wind field vector, temperature stratification parameter, atmospheric boundary layer height, turbulence intensity and precipitation intensity to form a sequence of regulation intensity values.

[0064] Specifically, measured values ​​of five meteorological factors at each time step are extracted from each grid point on the standard grid from the multi-source fusion dataset. For each grid point on the standard grid, the measured value of the wind field vector at that grid point is multiplied by the allocation coefficient corresponding to the wind field vector to obtain the weighted contribution value of the wind field vector. The measured value of the temperature stratification parameter at that grid point is multiplied by the allocation coefficient corresponding to the temperature stratification parameter to obtain the weighted contribution value of the temperature stratification parameter. The measured value of the atmospheric boundary layer height at that grid point is multiplied by the allocation coefficient corresponding to the atmospheric boundary layer height to obtain the weighted contribution value of the atmospheric boundary layer height. The measured value of the turbulence intensity at that grid point is multiplied by the allocation coefficient corresponding to the turbulence intensity to obtain the weighted contribution value of the turbulence intensity. The measured value of the precipitation intensity at that grid point is multiplied by the allocation coefficient corresponding to the precipitation intensity to obtain the weighted contribution value of the precipitation intensity.

[0065] Specifically, the spatial arrangement of the standard grid is determined, which specifies the fixed order of grid points from the first row and first column to the last row and last column in space. At the same time, the time step order is determined, which specifies the fixed arrangement of each time step from the start time to the end time. Then, the comprehensive meteorological constraint values ​​of all grid points and all time steps are organized in a fixed manner of time step first and grid point second.

[0066] Furthermore, the total control intensity is obtained by summing the five control intensity values ​​in the sequence. For each control intensity value in the sequence, the control intensity value is divided by the total control intensity, and the resulting ratio is the allocation coefficient corresponding to that sub-meteorological factor. This allocation coefficient reflects the contribution share of that sub-meteorological factor to the comprehensive meteorological constraint among all five meteorological factors. Each of the five sub-meteorological factors obtains an allocation coefficient, and the sum of the five allocation coefficients is equal to one. The above division operation is performed on each of the five control intensity values ​​in the sequence to obtain the allocation coefficient corresponding to each of the five sub-meteorological factors.

[0067] Furthermore, the weighted contribution values ​​of the above five meteorological factors at the same grid point are added together, and the sum is the comprehensive meteorological constraint value of that grid point at that time step. The above weighted summation calculation is performed on all grid points and all time steps on the standard grid to obtain the comprehensive meteorological constraint value of each grid point at each time step.

[0068] Furthermore, a spatial matrix is ​​assigned to each time step, and each position in the matrix corresponds to a grid point on the standard grid. The value at that position is the comprehensive meteorological constraint value of that grid point at that time step. The spatial matrices of all time steps are arranged sequentially according to the time step order to form a three-dimensional data structure. The first and second dimensions of this three-dimensional data structure correspond to the spatial rows and columns of the standard grid, respectively, and the third dimension corresponds to the time step order. Each element in this three-dimensional data structure corresponds to the comprehensive meteorological constraint value of a specific grid point at a specific time step. This three-dimensional data structure is used as a dynamic meteorological constraint matrix.

[0069] In summary, by arranging the regulation intensity values ​​of the five sub-meteorological factors in a specific meteorological process into a regulation intensity value sequence and calculating the sum of the five values ​​in the sequence, and then dividing each regulation intensity value in the sequence by the sum, each sub-meteorological factor obtains a distribution coefficient between zero and one, and the sum of the five distribution coefficients is exactly equal to one. The magnitude of this distribution coefficient directly reflects the relative regulatory contribution of the corresponding sub-meteorological factor to the change of pollutant concentration among all five meteorological factors in the specific meteorological process. The larger the distribution coefficient, the more significant the regulatory effect of the factor in the meteorological process. This conversion method unifies the five regulation intensity values, which originally had different dimensions and different numerical ranges, into dimensionless standard weight values ​​with clear relative meaning, completely eliminating the obstacle that the factors cannot be directly compared and combined due to their different original dimensions. At the same time, the distribution coefficient is calculated for a specific meteorological process type. When the meteorological process type changes from stable weather to severe convective weather, the distribution coefficients of each factor are automatically adjusted accordingly, realizing the adaptive configuration of meteorological factor weights for different meteorological scenarios.

[0070] In summary, for each grid point on the standard grid, the actual observed values ​​of the five meteorological factors at that grid point are multiplied by their respective allocation coefficients to obtain five weighted measured values. These five weighted measured values ​​are then summed at the same grid point to obtain the comprehensive meteorological constraint value for that grid point. This calculation method ensures that the comprehensive meteorological constraint value for each grid point is determined by the product of the actual observed values ​​of the five meteorological factors at that grid point and their respective allocation coefficients. Since there are spatial distribution differences in the measured values ​​of the five meteorological factors among different grid points, while the allocation coefficients remain constant throughout the entire region, the comprehensive meteorological constraint value adaptively and differentially distributes with the spatial distribution of the measured values ​​of each factor among the grid points. This achieves an organic combination of the relative contribution weights of meteorological factors at the process level and the measured spatial distribution of meteorological factors at the grid point level, generating a differentiated comprehensive constraint intensity for each grid point that reflects the actual meteorological conditions at that location.

[0071] In summary, by filling a two-dimensional spatial matrix with the comprehensive meteorological constraint values ​​of all grid points at each time step according to the fixed spatial arrangement order of the standard grid from the first row and first column to the last row and last column, the rows and columns of this matrix completely correspond to the spatial rows and columns of the standard grid. The value at each position in the matrix precisely represents the comprehensive meteorological constraint value of the corresponding grid point at that time step. Then, the two-dimensional spatial matrices of all time steps are arranged sequentially according to the fixed order of the time steps from the start time to the end time to form a three-dimensional data structure. Each element in this three-dimensional data structure corresponds to the comprehensive meteorological constraint value of a specific grid point at a specific time step. The meteorological constraint value, this dynamic meteorological constraint matrix, maintains a grid arrangement and spatial coverage that is completely consistent with the original pollutant concentration field and emission characteristic concentration field in the spatial dimension, and maintains a one-to-one correspondence with each time step in the temporal dimension. This ensures that in the subsequent meteorological correction steps, this matrix can be directly used as a constraint benchmark to perform grid-by-grid and time-step pairing operations with the pollutant concentration field. At the same time, the dynamic change characteristics of this matrix over time enable meteorological correction to capture the continuous regulatory effect of hourly evolution of meteorological conditions on the pollution process, solving the fundamental defect of the traditional static source tracing method where meteorological constraints are fixed and cannot adapt to dynamic pollution processes.

[0072] In this embodiment of the invention, the calculation formula for the comprehensive meteorological constraint value is specifically used for: in, The comprehensive meteorological constraint value is... For static and stable cumulative weighting coefficients, For atmospheric stability parameters, For temperature stratification parameters, For humidity parameters, Ventilation coefficient, The temperature and humidity saturation feedback coefficient is... For precipitation removal weighting coefficient, Given the current rainfall intensity, The attenuation coefficient for ventilation wet removal is... This is a reference value for the ventilation coefficient. The lagging precipitation coupling weighting coefficient, For lag Precipitation intensity The number of time steps is the lag time. This is the humidity suppression feedback coefficient.

[0073] Specifically, the static stability cumulative weighting coefficient, precipitation removal weighting coefficient, delayed precipitation coupling weighting coefficient, temperature and humidity saturation feedback coefficient, ventilation and moisture removal attenuation coefficient, and humidity suppression feedback coefficient are all derived from the sub-meteorological factor allocation coefficients calculated in the previous steps. These allocation coefficients are obtained by dividing the control intensity value of each sub-meteorological factor in a specific meteorological process by the sum of the control intensity values ​​of all five sub-meteorological factors. The allocation coefficient itself is the weighting coefficient. The atmospheric stability parameter, the temperature stratification parameter, the humidity parameter, the ventilation coefficient, the current precipitation intensity, and the delayed precipitation intensity are all derived from the measured values ​​of each sub-meteorological factor in the multi-source fusion dataset at each grid point and time step of the standard grid after dimensionless preprocessing. The ventilation coefficient reference value is derived from the statistical distribution characteristics of the ventilation coefficients of all grid points in the historical period of the target process area. The number of delayed time steps is derived from the time interval between the current calculation time and the precipitation occurrence time.

[0074] Furthermore, the comprehensive meteorological constraint value characterizes the comprehensive constraint strength on the diffusion, transformation, and removal of atmospheric pollutants at each grid point within the target process area at each time step, formed by the combined effects of multiple meteorological factors. The formula is calculated by dividing the meteorological factors into three components according to their different mechanisms of action on pollutants and then summing them up. The first part characterizes the combined constraint effect of atmospheric stability, temperature stratification, and humidity on the vertical diffusion inhibition and photochemical transformation of pollutants under the regulation of ventilation conditions. The second part characterizes the direct flushing and removal effect of current precipitation on pollutants under the attenuation regulation of wet removal efficiency under ventilation conditions. The third part characterizes the indirect transport effect of the residual humidity field after the end of the previous precipitation on the spatial distribution of pollutants under the drive of current ventilation conditions. The calculation results of the formula are directly used in the subsequent meteorological correction step to extract the meteorological contribution from the original pollutant concentration.

[0075] In summary, when the ventilation coefficient increases, the denominator of the first part increases, leading to a decrease in the overall value of the first part. The denominator of the second part, containing the ratio of the ventilation coefficient to its reference value, also decreases the overall value of the second part. The numerator of the third part, containing the product of the ventilation coefficient and the lagged precipitation intensity, increases the value of the third part. This indicates that the ventilation coefficient has different effects on the three types of meteorological constraints, suppressing direct scavenging and enhancing lagged transport. When the values ​​of the temperature stratification parameter and the humidity parameter increase simultaneously, the numerator of the first part increases, and the saturation feedback term in the denominator increases synchronously, causing the first part to approach its saturation upper limit. This indicates that the enhancing effect of temperature and humidity conditions on static and stable cumulative constraints has an upper limit rather than increasing indefinitely. When the current precipitation intensity... As the numerical value increases, the increase in the numerator of the second part leads to an increase in the overall value of the second part and approaches the upper limit of saturation, indicating that the precipitation removal effect is enhanced with the increase in precipitation intensity, but the rate of enhancement gradually slows down. When the value of the lagging precipitation intensity increases, the increase in the numerator of the third part leads to an increase in the overall value of the third part, indicating that the stronger the previous precipitation, the more significant the indirect transport effect of its residual humidity field under the current wind field. When the value of the current precipitation intensity increases and the value of the ventilation coefficient also increases, the increase in the numerator of the second part drives the second part to increase, while the ventilation attenuation term in the denominator simultaneously drives the second part to decrease. The two effects achieve an antagonistic dynamic trade-off through the fractional structure, and the final trend depends on the relative strength between the removal driving effect of precipitation enhancement and the removal efficiency attenuation effect of ventilation enhancement.

[0076] In this embodiment of the invention, when performing meteorological correction on the PM2.5 and O3 pollutant concentration fields of the multi-source fusion dataset and the dynamic emission inventory based on the dynamic meteorological constraint matrix to obtain the emission characteristic concentration field and net source inventory of the target process, it is specifically used for: Based on the multi-source fusion dataset, the concentration fields of PM2.5 and O3 pollutants on the baseline grid in the target process are extracted to obtain the original pollutant concentration field of the target process; The dynamic meteorological constraint matrix is ​​divided into transport constraint sub-matrices, diffusion constraint sub-matrices, accumulation constraint sub-matrices, and wet deposition constraint sub-matrices according to the mechanism of action of meteorological factors. Based on the transmission constraint submatrix, the diffusion constraint submatrix, the cumulative constraint submatrix, and the wet settlement constraint submatrix, the transmission contribution, diffusion contribution, cumulative contribution, and wet settlement reduction of the grid points on the reference grid are calculated respectively. The emission contribution concentration value of the target process is obtained by subtracting the weighted sum of the transport contribution, the diffusion contribution, the cumulative contribution, and the wet deposition reduction from the grid point concentration values ​​of the original pollutant concentration field. The emission contribution concentration value is then spatially reorganized to obtain the emission characteristic concentration field of the target process. Based on the comprehensive meteorological constraint value of the dynamic meteorological constraint matrix, the species emissions of the dynamic emission inventory are divided into a meteorological interference layer and an intrinsic layer. The emissions of the intrinsic layer are arranged according to the spatial organization of the baseline grid to obtain the net source inventory of the target process.

[0077] Specifically, the spatial coverage and grid division method of the baseline grid are determined. The baseline grid is completely consistent with the aforementioned standard grid in terms of spatial range, number of grid points, and grid spacing. Then, all records with PM2.5 concentration monitoring values ​​as the data type label are selected from the multi-source fusion dataset. According to the grid location identifier and time identifier carried by each record, the concentration values ​​in each record are filled into the corresponding time step of the corresponding grid point of the baseline grid, forming the concentration field of PM2.5 pollutants at all grid points and all time steps.

[0078] Specifically, the comprehensive meteorological constraint value at each time step of each grid point in the dynamic meteorological constraint matrix is ​​obtained. This comprehensive meteorological constraint value is a weighted composite value of five sub-meteorological factors: wind field vector, temperature stratification parameter, atmospheric boundary layer height, turbulence intensity, and precipitation intensity. Based on this, the sub-meteorological factors are classified according to their different mechanisms of action on the evolution of pollutants. The wind field vector is extracted separately to form a transport constraint sub-matrix. Each value in this sub-matrix represents the intensity of horizontal transport of pollutants driven by the horizontal wind field at that grid point at that time step. The atmospheric boundary layer height and turbulence intensity are combined and extracted to form a diffusion constraint sub-matrix.

[0079] Specifically, for each grid point and each time step on the reference grid, the transport constraint value for that grid point at that time step is extracted from the transport constraint submatrix. This value reflects the intensity of the horizontal wind field transporting pollutants from upstream to the current grid point, and is used as the transport contribution of that grid point. The diffusion constraint value for that grid point at that time step is extracted from the diffusion constraint submatrix. This value reflects the intensity of vertical turbulent mixing at the current grid point, which disperses pollutants from the ground layer to the upper atmosphere or mixes them from the upper atmosphere to the ground, and is used as the diffusion contribution of that grid point. The cumulative constraint value for that grid point at that time step is extracted from the cumulative constraint submatrix. This value reflects the intensity of stable meteorological conditions and photochemical reaction conditions promoting the local generation and retention of pollutants at the current grid point, and is used as the cumulative contribution of that grid point.

[0080] Specifically, for each grid point and each time step on the baseline grid, the original concentration value of that grid point at that time step is first obtained from the original pollutant concentration field. This original concentration value is considered as the result of the combined effects of meteorological and anthropogenic emission factors. Then, the transport contribution, diffusion contribution, cumulative contribution, and wet deposition reduction of that grid point at that time step are obtained. The total meteorological contribution is obtained by adding the transport contribution, diffusion contribution, cumulative contribution, and wet deposition reduction. This total meteorological contribution represents the concentration change at the current grid point at the current time step caused by the combined effects of four meteorological mechanisms: horizontal transport, vertical diffusion, local accumulation, and wet scavenging. The difference between the original concentration value and the total meteorological contribution is the emission contribution concentration value after removing meteorological influences.

[0081] Specifically, the emission data of each pollutant species at each grid point in the dynamic emission inventory are obtained. At the same time, the comprehensive meteorological constraint value at the grid point with the exact same location as each grid point is extracted from the dynamic meteorological constraint matrix. The emission data of each species at each grid point and the comprehensive meteorological constraint value of that grid point are registered according to the grid point location, so that the data record of each grid point contains both the comprehensive meteorological constraint value and the emission value of each species. For each pollutant species at each grid point, the emission value of that species at that grid point is split into two parts based on the comprehensive meteorological constraint value of that grid point as the dividing criterion.

[0082] Furthermore, all records labeled as O3 concentration monitoring values ​​were selected from the multi-source fusion dataset. The O3 pollutant concentration field was formed at all grid points and all time steps using the same filling method. The PM2.5 and O3 concentration fields were arranged one by one according to the grid point position and time step, and integrated into the original pollutant concentration field. Each element in the original pollutant concentration field corresponds to the set of PM2.5 and O3 pollutant concentration values ​​at a specific grid point and a specific time step. The field as a whole represents the original pollutant concentration distribution state of the target process area at all grid points and all time steps without meteorological correction.

[0083] Furthermore, each value in the submatrix represents the vertical diffusion intensity of pollutants driven by vertical turbulent mixing at that grid point at that time step. The temperature stratification parameter and humidity parameter are extracted to form an accumulation constraint submatrix. Each value in this submatrix represents the local generation and accumulation intensity of pollutants driven by both static and stable meteorological conditions and photochemical reaction conditions at that grid point at that time step. The precipitation intensity is extracted separately to form a wet deposition constraint submatrix. Each value in this submatrix represents the wet removal intensity of pollutants driven by precipitation scouring at that grid point at that time step. The above four submatrixes are consistent with the original pollutant concentration field in terms of the number of grid points, spatial range, and arrangement order. That is, the four submatrixes and the original pollutant concentration field have corresponding values ​​at each grid point location and at each time step.

[0084] Furthermore, the wet deposition constraint value at this time step of the grid point is extracted from the wet deposition constraint submatrix. This value reflects the intensity of the precipitation flushing and removing pollutants at the current grid point location, and is used as the wet deposition reduction amount of the grid point. The above extraction operation is performed on all grid points and all time steps on the reference grid one by one to obtain four types of values ​​for each grid point at each time step: transport contribution, diffusion contribution, cumulative contribution, and wet deposition reduction amount.

[0085] Furthermore, the emission contribution concentration value represents the concentration level contributed entirely by pollutants emitted by anthropogenic emission sources at the current grid point and current time step. The above difference calculation is performed on all grid points and all time steps on the reference grid to obtain the emission contribution concentration values ​​of PM2.5 and O3 at each grid point and each time step. The emission contribution concentration values ​​of all grid points and all time steps are reorganized according to the spatial arrangement order and time step order of the reference grid to form a concentration field that is completely consistent with the original pollutant concentration field in terms of spatial and temporal structure. This concentration field is used as the emission characteristic concentration field.

[0086] Furthermore, the portion proportional to the comprehensive meteorological constraint value is designated as the meteorological interference layer, and the remaining portion of the emission values ​​after deducting the meteorological interference layer is designated as the intrinsic layer. The emission values ​​marked as meteorological interference layer are removed from the emission records of that species at that grid point, while the emission values ​​marked as intrinsic layer are retained in the emission records of that grid point. After performing the above splitting, removal, and retention operations on all grid points and all pollutant species on the baseline grid, the intrinsic layer emission data retained by all grid points are placed into the corresponding spatial positions according to the spatial arrangement order of the baseline grid from the first row and first column to the last row and last column, forming a net source inventory.

[0087] In summary, by selecting all concentration monitoring records of PM2.5 and O3, two target pollutants, from the multi-source fusion dataset that has undergone spatiotemporal registration and gridding, and accurately filling the concentration values ​​into the corresponding positions of the baseline grid according to the grid location and time identifiers carried by each record, a complete concentration field was obtained that is completely consistent with the baseline grid in terms of spatial range, number of grid points, and grid spacing, and has concentration values ​​at all grid points and all time steps. This original pollutant concentration field directly serves as the benchmark for subsequent meteorological correction, providing a raw concentration data basis without any correction processing for removing meteorological interference, ensuring the traceability of concentration changes before and after correction and the quantitative assessability of correction effects. The baseline grid is completely consistent with the aforementioned standard grid in spatial structure, ensuring that the original pollutant concentration field and the dynamic meteorological constraint matrix and subsequent constraint sub-matrices can achieve accurate grid-by-grid pairing operations in the spatial dimension.

[0088] In summary, by extracting wind field vectors from the dynamic meteorological constraint matrix to form a transmission constraint sub-matrix, atmospheric boundary layer height and turbulence intensity to form a diffusion constraint sub-matrix, temperature stratification parameters and humidity parameters to form a cumulative constraint sub-matrix, and precipitation intensity to form a wet deposition constraint sub-matrix, the dynamic meteorological constraint matrix, originally a single comprehensive numerical value, is decomposed into four constraint sub-matrices with independent physical meanings according to the differentiated action mechanisms of different meteorological factors on the pollutant evolution process. Each sub-matrix represents the intensity distribution of the influence of a specific meteorological mechanism on the pollution process. The four sub-matrices are completely consistent with the original pollutant concentration field in terms of grid number, spatial range, and arrangement order. This provides independent constraint inputs that are spatially precisely corresponding to the original concentration field for subsequent calculation of different types of meteorological contributions, solving the fundamental problem in traditional techniques that mix multiple meteorological factors together and cannot distinguish their individual contributions.

[0089] In summary, for each grid point on the baseline grid, the constraint values ​​at the corresponding positions of the grid points are extracted from the four constraint sub-matrices as the transport contribution, diffusion contribution, accumulation contribution, and wet deposition reduction of that grid point. This allows each grid point to simultaneously obtain four types of values ​​representing the independent contribution intensity of four different meteorological mechanisms—horizontal transport, vertical diffusion, local accumulation, and wet scavenging—to the pollutant concentration at that grid point. The values ​​of the four types of contributions are directly derived from their respective constraint sub-matrices and have clear values ​​at the same grid point position and time step. This calculation method gives the originally abstract meteorological constraint matrix values ​​a clear meteorological process meaning. The relative magnitudes of the four types of contributions directly reflect the type of dominant meteorological mechanism at that grid point, providing a data foundation for the subsequent accurate separation of various meteorological contributions from the original concentration by mechanism-specific quantification, and avoiding correction errors caused by treating multiple meteorological contributions indiscriminately.

[0090] In summary, for each grid point on the baseline grid, the weighted sum of the concentration changes caused by the four meteorological mechanisms of transport, diffusion, accumulation, and wet deposition is subtracted from the original concentration value. The resulting difference is the concentration value contributed entirely by anthropogenic emission sources after removing all meteorological interference. This emission contribution concentration value eliminates the spurious concentration changes caused by meteorological fluctuations, ensuring that the concentration data used in subsequent source tracing modeling truly reflects the actual contribution level of emission sources at each grid point. The emission characteristic concentration field obtained by reorganizing the emission contribution concentration values ​​of all grid points according to the fixed spatial arrangement order of the baseline grid is completely consistent with the original pollutant concentration field in both spatial and temporal structure, but the information it contains has changed from the total concentration resulting from the combined effect of meteorology and emissions to the net concentration resulting from the effect of emissions alone. This emission characteristic concentration field serves as the concentration input in the subsequent construction of the spatiotemporal correlation matrix, fundamentally solving the technical defect of misjudging emission source contributions due to meteorological interference in traditional technologies.

[0091] In summary, by dividing the emissions of each species at each grid point in the dynamic emission inventory into a meteorological interference layer and an intrinsic layer based on the comprehensive meteorological constraint value, and removing the meteorological interference layer while retaining the intrinsic layer, the pseudo-change components caused by meteorological condition fluctuations in the emission data of each pollutant species at each grid point in the emission inventory are accurately identified and removed. Only the intrinsic layer data reflecting the inherent emission intensity of the emission source itself are retained. Since the comprehensive meteorological constraint value characterizes the comprehensive intensity of the influence of meteorological factors on each grid point under the current meteorological conditions and changes dynamically hourly, it is used as the division benchmark. This allows the intrinsic layer emissions to adaptively remove different proportions of meteorological interference spatially according to the differences in meteorological conditions at each grid point, and to dynamically adjust the removal magnitude temporally according to the changes in the comprehensive meteorological constraint value at each time step. This makes the net source inventory a set of independent emission datasets completely stripped of meteorological interference. This net source inventory and the emission characteristic concentration field constitute two sets of parallel correction products obtained by independently performing meteorological correction on the emission inventory and concentration field respectively, starting from the same set of dynamic meteorological constraint matrix and the same benchmark grid, ensuring that the meteorological correction operation has been completed in both the concentration data and emission data dimensions.

[0092] In this embodiment of the invention, the step of constructing a spatiotemporal correlation matrix based on the dynamic emission inventory and the emission characteristic concentration field, performing spatiotemporal convolution processing on the spatiotemporal correlation matrix to obtain the grid spatiotemporal convolution response value of the target process, and performing receptor-source analysis on the source grid data corresponding to the grid spatiotemporal convolution response value to obtain the collaborative pollution source tracing result of the target process is specifically used for: The gridded emissions of pollutant species in the dynamic emission inventory are expanded along the spatial dimension to obtain the emission space vector of the target process. The gridded concentration values ​​of PM2.5 and O3 in the emission characteristic concentration field are expanded along the spatial dimension to obtain the concentration space vector of the target process. The emission space vector and the concentration space vector are stacked in time series to obtain the emission time series space matrix and the concentration time series space matrix of the target process. Calculate the spatial cross-correlation matrix between the emission time series spatial matrix and the concentration time series spatial matrix, and extend the spatial cross-correlation matrix along the time step direction in the time dimension to obtain the spatiotemporal correlation matrix of the target process; The neighborhood association strength of the grid points in the spatiotemporal correlation matrix is ​​causally convolved in the time step direction to obtain the grid spatiotemporal convolution response value of the target process. Based on the distribution of the spatiotemporal convolution response values ​​of the grid points, the grid points whose response values ​​exceed the target and a set verification threshold are identified as pollution source grid points. Based on the location information of the pollution source grid points, the concentration time series of the pollutant species are extracted from the emission characteristic concentration field. The concentration time series is subjected to non-negative matrix decomposition to obtain the source contribution matrix and source spectrum matrix of the target process. Based on the source contribution matrix and the source spectrum matrix, the key precursor categories and pollution contribution of the target process are determined. The PM2.5 to O3 concentration ratio, temperature parameters, and radiation parameters of the pollution source grid points are obtained. Based on the temperature parameters and radiation parameters, the control type threshold of the target process is determined. Based on the control type threshold, the concentration ratio is judged to obtain the sensitivity type of the target process. The location and response value intensity of the pollution source grid points, the category of the key precursor, the pollution contribution, and the sensitivity type are combined to form the collaborative pollution source tracing result of the target process.

[0093] Specifically, the emission values ​​of all pollutant species at all grid points in the dynamic emission inventory are obtained. For each pollutant species covered by the dynamic emission inventory, the emission values ​​of the species at all grid points are extracted one by one according to the fixed grid point arrangement order of the baseline grid and arranged into a one-dimensional sequence. Each position of the one-dimensional sequence corresponds to a specific grid point on the baseline grid. The one-dimensional sequence is used as the emission space vector of the species. After performing the above expansion operation on all pollutant species, the emission space vectors corresponding to each species are obtained. These emission space vectors are arranged in parallel according to the species order to form a complete emission space vector set.

[0094] Specifically, for each pollutant species, the emission spatial vector contains only spatial dimension information and not temporal dimension information. The emission spatial vectors of the species at all time steps are obtained, and these emission spatial vectors are arranged in the order of time steps. The emission spatial vector of the first time step is used as the first row of the matrix, the emission spatial vector of the second time step is used as the second row of the matrix, and so on until all time steps are covered, forming a two-dimensional matrix. The row direction of the matrix corresponds to the time step order, and the column direction corresponds to the grid arrangement order. This matrix is ​​the emission time series spatial matrix of the species.

[0095] Specifically, for the emission time-series spatial matrix of each pollutant species, all columns of the matrix are taken, from the first column to the last column. Each column represents the emission time series of that species at a specific grid point across all time steps. Simultaneously, all columns of the concentration time-series spatial matrix are taken, each column representing the concentration time series of PM2.5 or O3 at a specific grid point across all time steps. Each column in the emission time-series spatial matrix is ​​paired with each column in the concentration time-series spatial matrix one by one. For each pair of paired column vectors, the two values ​​at the same time step in the two column vectors are multiplied and summed to obtain the correlation value of the pair. All the paired correlation values ​​are organized into a two-dimensional matrix according to the column indices of the emission spatial matrix and the concentration spatial matrix. The number of rows in this two-dimensional matrix is ​​equal to the number of columns in the emission time-series spatial matrix, i.e., the total number of grid points, and the number of columns is equal to the number of columns in the concentration time-series spatial matrix, i.e., the total number of grid points. This two-dimensional matrix is ​​the spatial cross-correlation matrix.

[0096] Specifically, for each grid point in the spatiotemporal correlation matrix, the sequence of spatial cross-correlation values ​​between the grid point and its neighboring grid points on the reference grid over all time offsets is extracted. For the grid point currently being processed, all spatially adjacent neighboring grid points on the reference grid are determined. The correlation values ​​between the current grid point and each neighboring grid point over all time offsets are extracted from the spatiotemporal correlation matrix. These correlation values ​​are arranged into a sequence in ascending order of time offset. Each position in this sequence corresponds to a specific time offset. This sequence represents the distribution of the influence intensity of the emission change of the current grid point on the concentration change of its neighboring grid points at different time lags.

[0097] Specifically, firstly, the spatiotemporal convolutional response values ​​of all grid points calculated in the previous steps are obtained. These response values ​​are then arranged into a complete response value spatial distribution matrix according to the position of their corresponding grid points on the reference grid. The central tendency and dispersion values ​​of all response values ​​in this matrix are calculated, and the central tendency and dispersion values ​​are added together to obtain the set test threshold.

[0098] Specifically, the location coordinates of all pollution source grid points are obtained. For each pollution source grid point, the concentration values ​​of all pollutant species at all time steps are extracted from the emission characteristic concentration field based on its location coordinates. These concentration values ​​are arranged into a two-dimensional concentration time series matrix according to the organization method of row direction corresponding to time steps and column direction corresponding to pollutant species. The number of rows of this matrix is ​​equal to the total number of time steps, and the number of columns is equal to the total number of pollutant species. The concentration time series matrices of all pollution source grid points are combined into a set of matrices to be decomposed.

[0099] Specifically, for each pollution source grid point, the PM2.5 concentration value and O3 concentration value at the same time step at that grid point are extracted from the emission characteristic concentration field based on its location coordinates. The PM2.5 concentration value is divided by the O3 concentration value to obtain the PM2.5 to O3 concentration ratio of that grid point. At the same time, the temperature observation value and radiation observation value at that grid point are extracted from the multi-source fusion dataset based on the location coordinates and the same time step as temperature parameters and radiation parameters. The temperature parameters and radiation parameters are input into the control type threshold conversion relationship for conversion to obtain the control type threshold.

[0100] Specifically, the location coordinates of all pollution source grid points obtained in the previous steps are organized into a location list according to the grid point arrangement order of the reference grid. The response value intensity values ​​of each pollution source grid point are organized into a response value intensity list according to the same grid point order as the location list. The key precursor categories of each pollution source grid point are organized into a key precursor list according to the same grid point order as the location list. The pollution contribution values ​​of each pollution source grid point are organized into a pollution contribution list according to the same grid point order as the location list.

[0101] Furthermore, the concentration values ​​of PM2.5 and O3 pollutants at all grid points and all time steps are extracted from the emission characteristic concentration field. For PM2.5 pollutant, the concentration values ​​of all grid points at each time step are extracted one by one and arranged into a one-dimensional sequence according to the fixed grid point arrangement order of the reference grid, which is exactly the same as the above. Each position of the one-dimensional sequence corresponds to a specific grid point on the reference grid. The one-dimensional sequence serves as the concentration space vector of PM2.5 at that time step. The O3 pollutant is extracted in the same way to obtain the concentration space vector of O3 at each time step. The concentration space vectors of PM2.5 and O3 at all time steps are organized according to the time step order to obtain the concentration space vector set.

[0102] Furthermore, the above-mentioned temporal stacking operation is performed on all pollutant species to obtain the emission temporal spatial matrix of each species. For the concentration spatial vector, PM2.5 and O3 each have a concentration spatial vector at each time step. These concentration spatial vectors are arranged in the order of the time steps, with the concentration spatial vector of the first time step as the first row of the matrix, the concentration spatial vector of the second time step as the second row of the matrix, and so on until all time steps are covered, forming a two-dimensional matrix. The row direction of this matrix corresponds to the time step order, and the column direction corresponds to the grid arrangement order. This matrix is ​​the concentration temporal spatial matrix, with PM2.5 and O3 each corresponding to a concentration temporal spatial matrix.

[0103] Furthermore, the spatial cross-correlation matrix is ​​extended along the time step direction in the time dimension. Multiple copies of the spatial cross-correlation matrix are made, each corresponding to a time offset. The spatial cross-correlation matrices under all time offsets are arranged in the order of the time offsets to form a three-dimensional data structure. This three-dimensional data structure is the spatiotemporal correlation matrix. The first and second dimensions of this matrix correspond to the spatial positional relationship between grid points, and the third dimension corresponds to the time offset.

[0104] Furthermore, a causal convolution operation is performed on the sequence in the time offset direction. The specific execution method of the causal convolution operation is as follows: for the current processing time offset position in the sequence, only the correlation values ​​corresponding to the time offsets before this position are used for weighted combination. The weights of the weighted combination are provided by a pre-set causal convolution kernel. The length of the convolution kernel is the same as the sequence length and has non-zero weights only at positions before the current processing position. The result of the weighted combination is used as the causal convolution output value at this time offset position. The above causal convolution operation is performed one by one for each time offset position in the sequence to obtain the causal convolution output sequence of the current grid point. The values ​​of each position in the causal convolution output sequence are added to obtain the spatiotemporal convolution response value of the current grid point. The above neighborhood correlation strength extraction and causal convolution operation are performed on each grid point on the reference grid to obtain the spatiotemporal convolution response value corresponding to each grid point.

[0105] Furthermore, each grid point on the baseline grid is traversed, and the spatiotemporal convolution response value of the grid point is compared with a set test threshold. If the response value of the grid point is greater than the set test threshold, the grid point is marked as a pollution source grid point. If the response value of the grid point is less than or equal to the set test threshold, the grid point is marked as a non-pollution source grid point and is not processed. After the traversal is completed, the location coordinates of all grid points marked as pollution sources and their corresponding response value intensities are extracted from the response value spatial distribution matrix. The response value intensity values ​​are recorded grid point by grid point using the grid point location coordinates as the index, forming a complete record set of pollution source grid points and their response value intensities.

[0106] Furthermore, for each concentration time series matrix in the set of matrices to be decomposed, a nonnegative matrix factorization operation is performed as input. Specifically, the concentration time series matrix is ​​approximately decomposed into the product of two low-dimensional nonnegative matrices under nonnegative constraints. In the first matrix, each row represents the contribution intensity of a factor at each time step, and each column represents the combination of contribution values ​​of each factor at a time step. In the second matrix, each row represents the loading distribution of a factor on each pollutant species, and each column represents the combination of loading values ​​of each factor on a pollutant species. After iterative optimization, the source contribution matrix and source spectrum matrix are obtained. VOCs-type factors and NO are identified from the source spectrum matrix. For each pollution source grid point, the contribution value of the grid point to the VOCs category factor is extracted from the source contribution matrix as the VOCs category load value, and the contribution value of the grid point to the NOx category factor is extracted as the NOx category load value. The VOCs category load value and the NOx category load value are compared, and the factor category with the higher value is taken as the key precursor category of the pollution source grid point. At the same time, the proportion of each factor load value in the sum of all factor load values ​​is taken as the pollution contribution degree of the emission source category corresponding to the pollution source grid point. The above decomposition, identification and calculation operations are performed on all pollution source grid points one by one to obtain the key precursor category and pollution contribution degree corresponding to each pollution source grid point.

[0107] Furthermore, the specific conversion method is as follows: the combined value of temperature parameter and radiation parameter is mapped to a numerical value through photochemical kinetic relationship. This numerical value is the control type threshold. The ratio of PM2.5 to O3 concentration is compared with the control type threshold. When the concentration ratio is less than the control type threshold, the sensitivity type of the pollution source grid point is determined to be VOC control type. When the concentration ratio is greater than the control type threshold, the sensitivity type of the pollution source grid point is determined to be NOx control type. The above ratio calculation, threshold conversion and discrimination operation are performed on all pollution source grid points one by one to obtain the sensitivity type corresponding to each pollution source grid point.

[0108] Furthermore, the sensitivity types of each pollution source grid point are organized into a sensitivity type list according to the same grid point order as the location list. Then, the location list, response value intensity list, key precursor list, pollution contribution list, and sensitivity type list are associated and combined according to the grid point location, so that the data record of each pollution source grid point simultaneously contains five types of information: the location coordinates of the grid point, the response value intensity of the grid point, the key precursor category of the grid point, the pollution contribution of the grid point, and the sensitivity type of the grid point. The complete data records of all pollution source grid points are summarized to form a collaborative pollution source tracing result.

[0109] In summary, by extracting and arranging emission and concentration data, originally stored in a two-dimensional grid matrix, into a one-dimensional sequence according to the fixed grid point arrangement of the baseline grid, each grid point's emission and concentration data has a unique and definite position index in the sequence. This position index forms a precise mapping relationship with the actual spatial position of the grid point. The emission spatial vector and concentration spatial vector are completely consistent in sequence structure, and the same position index corresponds to the same grid point. This establishes a data form that can be directly used for vector operations to calculate the spatial correlation between emissions and concentrations, solving the problem that cross-grid spatial correlation cannot be directly calculated under the original grid matrix form. At the same time, the time dimension is retained outside the vector sequence, allowing subsequent time-series stacking operations to be orderly expanded along the time direction while maintaining the spatial correspondence.

[0110] In summary, for the emission spatial vector of each pollutant species and the respective concentration spatial vectors of PM2.5 and O3, the vector of the first time step is used as the first row of the matrix, the vector of the second time step as the second row, and so on until all time steps are covered, forming a two-dimensional matrix in which the row direction corresponds to the time step sequence and the column direction corresponds to the grid arrangement sequence. This matrix structure simultaneously carries grid position information in the spatial dimension and evolution sequence information in the temporal dimension, so that emission data and concentration data that were originally independent at a single time step are organized into a data matrix with temporal continuity. The data changes between different rows in the same column of the matrix reflect the emission or concentration change trajectory of the same grid point over time, and the data differences between different columns in the same row reflect the spatial distribution differences of emission or concentration between different grid points at the same moment. This provides a data input structure that simultaneously contains information in both time and space dimensions for subsequent calculation of the spatiotemporal correlation between emissions and concentrations.

[0111] In summary, by pairing each column of the emission time-series spatial matrix of each pollutant species with each column of the concentration time-series spatial matrix, and multiplying the two values ​​at the same time step in each pair of column vectors and summing them, the correlation value of the pair is obtained. All the correlation values ​​of the pairs are organized into a two-dimensional spatial cross-correlation matrix with the number of rows and columns equal to the total number of grid points. This allows the spatial correlation between emission sources and affected concentration grid points to be quantitatively characterized. Subsequently, this spatial cross-correlation matrix is ​​extended along the time offset direction into a three-dimensional spatiotemporal correlation matrix, which fully preserves the characteristics of the change in correlation strength between different spatial locations with the degree of time lag. This spatiotemporal correlation matrix contains three dimensions: spatial location relationship information, time lag relationship information, and correlation strength information. It provides a data foundation containing complete spatiotemporal correlation information for subsequent extraction of pollution source grid points with temporal causal relationships through causal convolution.

[0112] In summary, for each grid point in the spatiotemporal correlation matrix, the correlation strength sequence between that grid point and all its spatially adjacent grid points at different time offsets is extracted. Then, a causal convolution operation is performed on this sequence in the direction of the time offset. This causal convolution operation only uses the correlation strength corresponding to the time offset before the current processing position for weighted combination, ensuring causal consistency in the time direction. That is, only the emission information of the past moment is used to predict the concentration response at the current moment, avoiding the non-causal distortion problem of using information from the future moment for retrospective prediction. The spatiotemporal convolution response value obtained by each grid point after the causal convolution operation integrates the dual information of spatial neighborhood relationship and temporal causal relationship. The magnitude of the response value accurately reflects the comprehensive influence intensity of the grid point as a pollution source on the concentration changes of surrounding grid points throughout the target process. Compared with the traditional approach of only using the spatial correlation of the same period or only using the temporal series correlation for analysis, this response value realizes the joint extraction and fusion expression of spatiotemporal correlation information.

[0113] In summary, the sum of the central tendency and dispersion of the spatiotemporal convolutional response values ​​of all grid points is used as the set test threshold. Grid points with response values ​​exceeding this threshold are marked as pollution source grid points. The determination of this test threshold is entirely based on the distribution characteristics of the current data itself. Under different regions, time periods, and meteorological conditions, this threshold can be automatically adjusted as the distribution of response values ​​changes. This effectively avoids false pollution source misjudgments or missed detections of real pollution sources caused by fixed thresholds. It also helps to screen out grid point objects with clear spatial orientation for subsequent chemical analysis steps, significantly narrowing the analysis scope of PMF and OBM from all grid points in the region to pollution source grid points that have been located by spatiotemporal convolution. This avoids ineffective calculations on non-pollution source grid points, improving overall computational efficiency and analytical focus.

[0114] In summary, by extracting the multi-species concentration time series at each grid point of the pollution source located by spatiotemporal convolution and performing non-negative matrix decomposition, multiple factors and their corresponding source contribution matrices and source spectrum matrices are separated from the mixed pollutant concentration signals, so that each factor corresponds to an emission source type. By identifying VOCs and NOx factors from the source spectrum matrix and comparing the loading values ​​of each pollution source grid point on the two types of factors, the factor category with the higher loading value is regarded as the key precursor category of the grid point. At the same time, the pollution contribution is calculated by using the proportion of each factor loading value in the total loading value. This allows the source tracing results, which originally only contained location and response intensity information, to obtain information on the material composition of emissions at the grid point, namely whether VOCs or NOx are dominant, and achieves independent quantification of the contribution share of industrial stationary sources, mobile sources of motor vehicles, ecological natural sources, and dust sources.

[0115] In summary, by acquiring the PM2.5 to O3 concentration ratio, temperature parameters, and radiation parameters of each pollution source grid point, and converting the temperature and radiation parameters into control type thresholds through photochemical kinetics, this method adaptively changes with the differences in temperature and radiation conditions at different grid points. By comparing the PM2.5 to O3 concentration ratio with this threshold, it determines whether the grid point is VOC-controlled or NOx-controlled sensitive. This provides each pollution source grid point in the source tracing results with control direction information on whether O3 generation at that grid point should prioritize VOCs or NOx control, providing a clear direction for precursor control in formulating differentiated and coordinated emission reduction strategies.

[0116] In summary, by linking and combining five independent data points—spatial location of pollution source grid points, response intensity, key precursor category, pollution contribution, and photochemical sensitivity type—according to grid point location, each pollution source grid point record simultaneously contains complete information across five dimensions: pollution source, pollution source contribution intensity, dominant emission type of the grid point, contribution levels of industrial, motor vehicle, dust, and natural sources, and control type of the grid point. This integrates information originally scattered across multiple independent data tables into a single complete record, providing comprehensive data support for developing coordinated emission reduction strategies based on the specific characteristics of each grid point in subsequent operational applications.

[0117] In this embodiment of the invention, when the collaborative pollution source tracing results are overlaid with the dynamic emission inventory in a grid-like manner to obtain the source tracing business results of the target process, it is specifically used for: The response values ​​of the source tracing grid points in the collaborative pollution source tracing results are matched with the pollutant species emission amounts at the same grid point in the dynamic emission inventory to obtain the pollution source matching table for the target process; According to the pollution source matching table, the pollutant proportion set of the source tracing grid points is extracted from the dynamic emission inventory, and the precursor category corresponding to the pollutant proportion set is used as the key precursor identifier of the target process; The location information of the source tracing grid points is spatially overlaid with the emission source category information in the dynamic emission inventory to identify the emission source category to which the source tracing grid points belong. The proportion of the response value corresponding to the emission source category in the total response value of the source tracing grid is taken as the pollution contribution of the target process; The location information of the source tracing grid points, the identification of the key precursors, the emission source category, and the pollution contribution are combined to form the source tracing business results of the target process.

[0118] Specifically, the location coordinates of all source tracing grid points and their corresponding response values ​​are extracted from the collaborative pollution source tracing results. At the same time, the emission data of all pollutant species at all grid points are extracted from the dynamic emission inventory. For each source tracing grid point, the emission records of pollutant species at the same grid point location are searched in the dynamic emission inventory using the location coordinates of the grid point. The response value of the source tracing grid point is matched one by one with the found emission records, with the response value record first and the emission record second to form a matching record. After performing the above search and matching operations on all source tracing grid points one by one, the matching records corresponding to each of the source tracing grid points are obtained. All matching records are summarized to form a pollution source matching table.

[0119] Specifically, for each source tracing grid point recorded in the pollution source matching table, the emission values ​​of all pollutant species at that grid point are obtained from the emission records corresponding to that grid point. The emission values ​​of various volatile organic compounds are summed to obtain the total volatile organic compound emission of that grid point, and the emission values ​​of various nitrogen oxide species are summed to obtain the total nitrogen oxide emission of that grid point. At the same time, the total emission of all species at that grid point is obtained.

[0120] Specifically, the location coordinates of all source tracing grid points are obtained from the collaborative pollution source tracing results, and the location coordinates of all grid points and their corresponding emission source category labels are extracted from the dynamic emission inventory. The emission source category label records the emission source category to which each grid point in the dynamic emission inventory belongs. The emission source categories include industrial stationary source category, motor vehicle mobile source category, dust source category and ecological natural source category.

[0121] Specifically, based on the correspondence between the full traceability grid points obtained in the previous step and their respective emission source categories, all traceability grid points are grouped according to their respective emission source categories. All traceability grid points belonging to the industrial stationary source category are assigned to the industrial stationary source group, all traceability grid points belonging to the motor vehicle mobile source category are assigned to the motor vehicle mobile source group, all traceability grid points belonging to the dust source category are assigned to the dust source group, and all traceability grid points belonging to the ecological natural source category are assigned to the ecological natural source group.

[0122] Specifically, the location coordinates of all source tracing grid points are organized into a source tracing grid point location list according to the grid point arrangement order of the baseline grid. The key precursor identifiers corresponding to each source tracing grid point are organized into a key precursor identifier list according to the same grid point order as the location list. The emission source categories corresponding to each source tracing grid point are organized into an emission source category list according to the same grid point order as the location list. The pollution contribution of each of the four emission source categories is organized into a pollution contribution list according to the fixed order of industrial stationary sources, motor vehicle mobile sources, dust sources, and ecological natural sources.

[0123] Furthermore, each row in the pollution source matching table corresponds to a source tracing grid point. Each row records the location coordinates of the grid point, the response value of the grid point, and the emission values ​​of various pollutant species at the grid point. This matching table establishes a direct correspondence between the source tracing response value and the emission of each species, providing a data foundation for the subsequent extraction of the pollutant composition ratio of each grid point.

[0124] Furthermore, the proportions of volatile organic compound (VOC) emissions and nitrogen oxide (NOx) emissions in the total emissions are calculated separately. These two proportions are used as the pollutant proportion set for that grid point. The above summation and proportion calculation operation is performed on each source tracing grid point in the pollution source matching table to obtain the pollutant proportion set corresponding to each source tracing grid point. Then, the pollutant proportion set of each source tracing grid point is identified to determine whether the proportion of VOCs in the proportion set is higher than that of NOx. The precursor category with the higher proportion is used as the key precursor identifier of that source tracing grid point. After performing the above identification operation on all source tracing grid points one by one, the key precursor identifier corresponding to each source tracing grid point is obtained.

[0125] Furthermore, each grid point in the dynamic emission inventory is pre-labeled with its corresponding emission source category. Then, for each source tracing grid point, the location coordinates of the source tracing grid point are compared one by one with the location coordinates of all grid points in the dynamic emission inventory. Grid points in the dynamic emission inventory with the same location coordinates as the source tracing grid points are found. The emission source category to which the grid point belongs is read from the emission source category label of the grid point. The emission source category is identified as the emission source category to which the source tracing grid point belongs. The above location comparison and emission source category reading operations are performed one by one for all source tracing grid points to obtain the emission source category corresponding to each source tracing grid point.

[0126] Furthermore, for each group, the response values ​​of all traceability grid points within the group are summed to obtain the total response value corresponding to that emission source category. At the same time, the response values ​​of all traceability grid points are summed to obtain the total response value of the traceability grid points. Then, for each emission source category, the total response value corresponding to that emission source category is divided by the total response value of the traceability grid points. The resulting ratio is the pollution contribution of that emission source category. The above grouping, summing, and division operations are performed on the four emission source categories respectively to obtain the pollution contribution of each emission source category.

[0127] Furthermore, the four lists of source tracing grid locations, key precursor identifiers, emission source categories, and pollution contribution are linked and combined to form a complete dataset containing all source tracing grids and their corresponding attribute information. Each record in this dataset contains the location coordinates of a source tracing grid, the key precursor identifier of that grid, the emission source category of that grid, and the pollution contribution of that emission source category. This complete dataset is used as the source tracing business outcome.

[0128] In summary, by utilizing the location coordinates of source tracing grid points to find emission records at the same grid point locations in the dynamic emission inventory and matching the response values ​​with the emission records to form matching records, the abstract response value of each source tracing grid point is anchored to the specific emission data of various pollutant species at that grid point. The response value represents the comprehensive contribution intensity of the grid point as a pollution source, while the emission data reveals the material composition of this contribution. After summarizing all matching records into a pollution source matching table, this matching table establishes a direct correspondence between source tracing response values ​​and the emission amounts of each species. This endows the source tracing results, which originally only contained location and response value intensity information, with an explanation of the material composition at the emission level. This solves the operational application problem of traditional source tracing results, which only provide spatial location but cannot answer what substances are mainly emitted by the grid point. It provides a precise matching foundation at the emission data level for subsequent extraction of pollutant composition ratios of each grid point and identification of key precursors.

[0129] In summary, for each source tracing grid, the total emissions of volatile organic compounds (VOCs) are obtained by summing the emissions of all VOC species and the total emissions of nitrogen oxides (NOx) species. Then, the proportion of each of these two types of precursors in the total emissions is calculated to obtain the pollutant proportion set. This proportion set accurately reflects the relative dominance of the two key precursors in the emission characteristics of this grid with two normalized proportion values. Subsequently, the precursor category with the higher proportion is used as the key precursor identifier for this grid, so that each source tracing grid has a clear precursor control direction.

[0130] In summary, by comparing the location coordinates of each source tracing grid point with the location coordinates of all grid points in the dynamic emission inventory one by one, and retrieving the category of the grid point from the emission source category label after finding the same location grid point, each source tracing grid point is accurately classified into a specific emission source category among industrial stationary sources, mobile vehicle sources, dust sources, or ecological natural sources. This classification operation transforms the source tracing grid points from abstract spatial coordinates and response values ​​into actual emission source forms with clear industry attributes and control entities. Industrial stationary sources correspond to specific industrial enterprises, mobile vehicle sources correspond to traffic control measures, dust sources correspond to construction dust control, and ecological natural sources correspond to vegetation management strategies. This provides category affiliation information that can be directly mapped to actual control objects for subsequent aggregation of contribution by emission source category and generation of classification control strategies.

[0131] In summary, by grouping all source tracing grid points according to their respective emission source categories, the response values ​​of all grid points within each group are summed to obtain the total response value corresponding to that category. This sum is then divided by the total response value of all source tracing grid points, resulting in a normalized pollution contribution value for each emission source category. This value accurately reflects the overall contribution share of that emission source category to co-polluting pollution across all emission source categories. The sum of the pollution contribution values ​​of the four emission source categories equals one, and the magnitude of each value directly indicates the relative contribution ranking of each source category. This quantitative result provides direct data support for precise emission reduction in regional coordinated prevention and control of air pollution.

[0132] In summary, the four independent data outputs—spatial location information from the source tracing grid location list, substance control guidance information from the key precursor identification list, control entity information from the emission source category list, and source category weight information from the pollution contribution list—are combined into a comprehensive dataset containing all source tracing grids and their complete attribute information. Each record in this dataset simultaneously includes the pollution source location, the precursor to be controlled at that location, the type of emission source at that location, and the contribution data of that type of source—providing complete decision-making information across four dimensions. This source tracing operational outcome is transformed from a simple source tracing analysis conclusion into an intelligent control tool that can be directly used by environmental management departments.

[0133] In this embodiment of the invention, when the proportion of the response value corresponding to the emission source category in the total response value of the source tracing grid is used as the pollution contribution of the target process, it is specifically used for: The response values ​​of the source grid points are accumulated over time steps to obtain the cumulative response intensity of the target process; The grid points on the reference grid in the target process are sorted according to the cumulative response intensity to obtain the cumulative response intensity sequence of the target process, and the grid points in the cumulative response intensity sequence whose response intensity reaches the cumulative contribution threshold are used as the hot spot screening grid point set of the target process. Spatially match the location of the hotspot screening grid points with the emission source categories marked in the dynamic emission inventory to obtain the actual emission source identifier of the target process; Extract the component emissions corresponding to the actual emission source identification from the dynamic emission inventory, and calculate the contribution value of the hot spot screening grid points to the common precursors of PM2.5 secondary formation and O3 photochemical pollution based on the component emissions; The proportion of the contribution value corresponding to the emission source category in the contribution value of the common precursors is used as the pollution contribution of the target process.

[0134] Specifically, the response values ​​of all source tracing grid points at each time step are obtained from the collaborative pollution source tracing results. For each source tracing grid point, the response value of the grid point at the first time step is used as the initial accumulation value. Then, the response value of the grid point at the second time step is added to the initial accumulation value to obtain the first accumulation result. Then, the response value at the third time step is added to the first accumulation result to obtain the second accumulation result. The response values ​​of the grid point at all time steps are accumulated one by one in the above manner until the response value of the grid point at the last time step is added to the accumulation result. The final sum is the cumulative response intensity of the grid point.

[0135] Specifically, all grid points on the reference grid are arranged in descending order of their respective cumulative response intensities. The grid point with the largest cumulative response intensity is placed at the first position of the sequence, the grid point with the second largest cumulative response intensity is placed at the second position of the sequence, and so on until the grid point with the smallest cumulative response intensity is placed at the last position of the sequence, thus forming a cumulative response intensity sequence.

[0136] Specifically, the location coordinates of each grid point are obtained from the hotspot screening grid point set, and the location coordinates of all grid points and their corresponding emission source category labels are obtained from the dynamic emission inventory. The emission source category label contains the emission source category information to which each grid point belongs. Then, for each grid point in the hotspot screening grid point set, the location coordinates of the grid point are compared one by one with the location coordinates of all grid points in the dynamic emission inventory. The grid point with the same location coordinates as the grid point is found in the dynamic emission inventory. The emission source category labeled by the grid point is read from the emission source category label of the grid point, and the category is used as the actual emission source identifier corresponding to the hotspot screening grid point.

[0137] Specifically, for each grid point in the hotspot screening grid set, the emission data of all pollutant species at the grid point location is retrieved from the dynamic emission inventory based on the actual emission source identifier corresponding to that grid point. From this emission data, the emission values ​​of various volatile organic compounds and various nitrogen oxide species that contribute to the secondary formation of PM2.5 and photochemical pollution of O3 are extracted. The emission values ​​of all volatile organic compounds are summed to obtain the total volatile organic compound emission of that grid point, and the emission values ​​of all nitrogen oxide species are summed to obtain the total nitrogen oxide emission of that grid point.

[0138] Specifically, based on all the hotspot screening grid points obtained in the previous step and their corresponding common precursor contribution values ​​and actual emission source identifiers, all hotspot screening grid points are grouped according to the emission source category represented by their actual emission source identifiers. All grid points belonging to the same emission source category are grouped into the same group. Then, for each group, the common precursor contribution values ​​of all grid points in the group are added together to obtain the contribution value corresponding to that emission source category. At the same time, the common precursor contribution values ​​of all hotspot screening grid points are added together to obtain the total precursor contribution value.

[0139] Furthermore, the cumulative response intensity reflects the total intensity accumulated by the spatial convolution response of the traceability grid point over the entire time span from the start to the end of the target process. The above time step accumulation operation is performed on all traceability grid points one by one to obtain the cumulative response intensity corresponding to each traceability grid point.

[0140] Furthermore, each position in the sequence corresponds to a grid point, and each grid point carries a specific value of its cumulative response intensity. Then, grid points are selected sequentially from the sequence in descending order of cumulative response intensity. Each time a grid point is selected, its cumulative response intensity is added to a temporary accumulator. When the accumulated value in the temporary accumulator reaches a preset cumulative contribution threshold, the selection stops. All selected grid points together constitute a hotspot screening grid point set. The grid points in this hotspot screening grid point set are the collection of all grid points whose cumulative response intensity contribution reaches the preset contribution threshold, representing key grid point areas with high cumulative response intensity and large pollution contribution in the target process.

[0141] Furthermore, the above-mentioned position comparison and category reading operations are performed on all grid points in the hotspot screening grid point set one by one to obtain the actual emission source identifier corresponding to each hotspot screening grid point. The actual emission source identifier represents the pollution source type of the grid point in actual emissions.

[0142] Furthermore, the weighted combination of total volatile organic compound (VOC) emissions and total nitrogen oxide (NOx) emissions is used as the common precursor contribution value for this grid point. This common precursor contribution value quantifies the potential contribution of the two key precursors, VOC and NOx, emitted at this grid point to the secondary formation of PM2.5 and O3 photochemical pollution. The above extraction and calculation operations are performed on all grid points in the hot spot screening grid point set one by one to obtain the common precursor contribution value corresponding to each hot spot screening grid point.

[0143] Furthermore, for each emission source category, the contribution value corresponding to that emission source category is divided by the total contribution value of common precursors. The resulting ratio is the pollution contribution of that emission source category. This pollution contribution reflects the share of that emission source category in the total contribution of common precursors of PM2.5 secondary formation and O3 photochemical pollution among all emission source categories. The above summation and division operation is performed on all emission source categories one by one to obtain the pollution contribution corresponding to each emission source category.

[0144] In summary, by summing all the response values ​​of each source grid point from the first time step to the last time step to obtain a comprehensive cumulative response intensity value, the response value information originally scattered across multiple time steps is integrated into a single value representing the total pollution contribution level of that grid point throughout the entire time span of the target process. This cumulative response intensity reflects the comprehensive characteristics of the grid point as a pollution source in two dimensions: duration and contribution intensity. A grid point with a high cumulative response intensity may be a short-term high-intensity emission or a long-term moderate-intensity emission. In either case, the value represents the overall degree of cumulative impact of that grid point on the surrounding pollution throughout the target process. Compared with the practice of using only a single time step response value or an average response value for hotspot screening, this cumulative response intensity more comprehensively characterizes the overall contribution level of each grid point throughout the target process, avoiding overestimation or underestimation caused by short-term fluctuations. It provides a comprehensive evaluation index that has been fully integrated in the time dimension for the subsequent extraction of hotspot screening grid points.

[0145] In summary, all grid points on the baseline grid are arranged in descending order of cumulative response intensity to form a cumulative response intensity sequence. Then, starting from the first point of the sequence, grid points are selected sequentially in descending order, and the cumulative response intensity of the selected grid points is continuously accumulated until the accumulated value reaches a preset cumulative contribution threshold. This method of accumulating from largest to smallest until the set threshold is reached ensures that the extracted hotspot screening grid point set covers the largest share of cumulative contribution with the fewest number of grid points. This ensures that the grid points selected into the hotspot screening grid point set are the group of grid points with the most prominent cumulative contribution and the highest proportion of overall pollution contribution to the target process. This screening method excludes many low-contribution grid points with cumulative response intensities below the threshold, allowing subsequent spatial matching and contribution value calculation operations to focus on key grid points that are truly of control significance. This significantly reduces the amount of data processed and solves the efficiency defect of traditional methods that treat all grid points equally and fail to highlight key control objects.

[0146] In summary, by comparing the location coordinates of each grid point in the hotspot screening grid with the location coordinates of all grid points in the dynamic emission inventory one by one, and reading the pre-labeled emission source category after finding grid points at the same location, each hotspot screening grid point is given a specific emission source category attribute instead of being an abstract grid coordinate. This actual emission source identification clearly answers the core operational question of "whether this key grid point corresponds to an industrial stationary source, a mobile source of motor vehicles, a dust source, or an ecological natural source." This makes the hotspot screening grid point no longer just a spatial point, but a manageable object that can be tracked to a specific emission source type. This provides a category classification basis for subsequent aggregation of contribution values ​​by emission source category and generation of source-specific control strategies. At the same time, this matching operation is entirely based on the existing labeling information in the dynamic emission inventory, without introducing additional classification uncertainty, thus ensuring the accuracy and reliability of the actual emission source identification.

[0147] In summary, for each grid point in the hotspot screening grid, the emissions of all volatile organic compounds (VOCs) and nitrogen oxides (NOx) at that location are extracted from the dynamic emission inventory based on the actual emission source identifier corresponding to that grid point. Then, the total VOC emissions for that grid point are summed, and the total NOx emissions are summed to obtain the total NOx emissions for that grid point. Finally, the common precursor contribution value for that grid point is obtained through the combination of the total VOC emissions and the total NOx emissions. This contribution value quantifies the potential contribution of VOCs and NOx emissions from that grid point to the secondary formation of PM2.5 and photochemical pollution of O3. This transforms the component emission data, which originally only contained emission quantity information, into a contribution indicator directly serving the goal of synergistic pollution control. This common precursor contribution value considers the synergistic effect of both VOCs and NOx precursors, rather than focusing solely on a single precursor, reflecting the core concept of tracing the source of synergistic pollution of PM2.5 and O3.

[0148] In summary, all grid points in the hotspot screening grid are grouped according to the emission source categories represented by the actual emission source identifiers. The common precursor contribution values ​​of all grid points in the same group are added together to obtain the contribution value corresponding to that emission source category. Then, the contribution value of each category is divided by the total common precursor contribution value of all hotspot screening grid points, so that each emission source category obtains a normalized pollution contribution. This pollution contribution accurately reflects the share of that emission source category in the total precursor contribution of PM2.5 secondary formation and O3 photochemical pollution among all emission source categories. The sum of the pollution contribution values ​​of the four emission source categories equals one, and each value directly gives the relative importance ranking of each type of source from the perspective of common precursors. This contribution value and the pollution contribution value calculated based on the response value in the previous step reflect the degree of contribution of emission sources from different perspectives. The two complement and verify each other, and together provide a quantitative decision-making basis for the formulation of precise emission reduction strategies.

[0149] In the several embodiments provided by this invention, it should be understood that the disclosed method can be implemented in other ways.

[0150] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0151] The embodiments of this application can acquire and process relevant data based on artificial intelligence technology. Artificial intelligence is the theory, method, and technology that uses digital computers or machines controlled by digital computers to simulate, extend, and expand human intelligence, perceive the environment, acquire knowledge, and use that knowledge to obtain optimal results.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A precise source tracing method for synergistic pollution of PM2.5 and O3 based on multi-source heterogeneous data fusion, characterized in that, The method includes: S1. Spatiotemporally resample the multi-source heterogeneous atmospheric environment data of the target process to obtain the multi-source fusion dataset of the target process; S2. Using the multi-source fusion dataset as an observation constraint, the emission inventory sample set of the target process is iteratively updated through integrated filtering to obtain the dynamic emission inventory of the target process; S3. Normalize the control intensity of the meteorological factors in the multi-source fusion dataset during a specific meteorological process to obtain the dynamic meteorological constraint matrix of the target process. S4. Based on the dynamic meteorological constraint matrix, perform meteorological correction on the PM2.5 and O3 pollutant concentration fields of the multi-source fusion dataset and the dynamic emission inventory to obtain the emission characteristic concentration field and net source inventory of the target process. S5. Construct a spatiotemporal correlation matrix based on the dynamic emission inventory and the emission characteristic concentration field, perform spatiotemporal convolution processing on the spatiotemporal correlation matrix to obtain the grid spatiotemporal convolution response value of the target process, and perform receptor source analysis on the source grid data corresponding to the grid spatiotemporal convolution response value to obtain the collaborative pollution source tracing result of the target process. S6. Overlay the collaborative pollution source tracing results with the dynamic emission inventory in a grid format to obtain the source tracing business results of the target process.

2. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, When performing spatiotemporal resampling on the multi-source heterogeneous atmospheric environment data of the target process to obtain the multi-source fused dataset of the target process, the following steps are included: Acquire atmospheric multi-source monitoring data for the target process, including fixed-site monitoring data, mass spectrometer component monitoring data, UAV vertical profile data, meteorological observation data, satellite remote sensing inversion data, and enterprise online emission time series data; By statistical outlier detection, abnormal equipment values ​​and extreme meteorological interference values ​​are removed from the atmospheric multi-source monitoring data to obtain the cleaned multi-source data stream of the target process; Spatial interpolation is used to fill in the missing data points in the multi-source data stream being cleaned, thereby obtaining the gridded data field of the target process; Align the timestamps of the data layers in the gridded data field to the same time base, and project the spatial coordinates of the data layers to the same spatial reference system to obtain the multi-source fusion dataset of the target process.

3. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, The step of using the multi-source fusion dataset as an observation constraint and iteratively updating the emission inventory sample set of the target process through ensemble filtering to obtain the dynamic emission inventory of the target process includes: Obtain an initial emission inventory sample set for the target process, the initial emission inventory sample set containing emission inventory samples generated based on the regional pollution source ledger; The ground component monitoring data, UAV vertical profile data and satellite remote sensing data of the multi-source fusion dataset are combined into the observation constraint vector of the target process; The residual between the simulated concentration values ​​of the samples in the initial emission inventory sample set on a unified spatial grid and the observation constraint vector is calculated to obtain the deviation sequence of the target process; The Kalman gain matrix of the initial emission inventory sample set is calculated based on the deviation sequence, and the emission amount of the sample is weighted and corrected using the Kalman gain matrix to obtain the assimilation sample set of the target process. The average value of the assimilated sample set is calculated along the sample dimension to obtain the dynamic emission inventory of the target process.

4. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, The step of normalizing the regulation intensity of meteorological factors in the multi-source fusion dataset during a specific meteorological process to obtain the dynamic meteorological constraint matrix of the target process includes: The wind field vector, temperature stratification parameters, atmospheric boundary layer height, turbulence intensity time series observation sequence, and precipitation intensity time series observation sequence of the multi-source fusion dataset are used as the meteorological factor time series of the target process. The meteorological process type is identified by performing meteorological process type identification on the time series of the meteorological factors to obtain the specific meteorological process type of the target process. The specific meteorological process type includes stable weather process, drought and high temperature process and severe convective weather process. Extract time-period subsequences corresponding to the specific meteorological process type from the time-series sequence of the meteorological factors, and use the correlation coefficients between the meteorological factors in the time-period subsequences and the concentrations of PM2.5 and O3 pollutants in the corresponding time period as the regulation intensity value of the target process; Based on the total control intensity of the sub-meteorological factors, the control intensity value is transformed to obtain the allocation coefficient of the target process, and the allocation coefficient is mapped to the grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process.

5. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 4, characterized in that, The process involves converting the control intensity value based on the total control intensity of the sub-meteorological factors to obtain the allocation coefficient of the target process, and mapping the allocation coefficient to grid points on the standard grid of the target process to obtain the dynamic meteorological constraint matrix of the target process, including: Obtain the sequence of regulation intensity values ​​of the sub-meteorological factors in the specific meteorological process, and use the ratio of the regulation intensity value in the sequence of regulation intensity values ​​to the sum of the regulation intensity values ​​of the sub-meteorological factors as the allocation coefficient of the target process; The comprehensive meteorological constraint value of the target process is calculated based on the measured values ​​of the sub-meteorological factors and the allocation coefficients; The comprehensive meteorological constraint values ​​are combined according to the spatial arrangement and time step order of the standard grid to obtain the dynamic meteorological constraint matrix of the target process.

6. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 5, characterized in that, The formula for calculating the comprehensive meteorological constraint value is as follows: in, The comprehensive meteorological constraint value is... For statically stable cumulative weighting coefficients, For atmospheric stability parameters, For temperature stratification parameters, For humidity parameters, Ventilation coefficient, The temperature and humidity saturation feedback coefficient is... For precipitation removal weighting coefficient, Given the current rainfall intensity, The attenuation coefficient for ventilation wet removal is... This is a reference value for the ventilation coefficient. The weighting coefficient for delayed precipitation coupling. For lag Precipitation intensity The number of time steps is the lag time. This is the humidity suppression feedback coefficient.

7. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, The process involves performing meteorological correction on the PM2.5 and O3 pollutant concentration fields and the dynamic emission inventory of the multi-source fusion dataset based on the dynamic meteorological constraint matrix to obtain the emission characteristic concentration field and net source inventory of the target process, including: Based on the multi-source fusion dataset, the concentration fields of PM2.5 and O3 pollutants on the baseline grid in the target process are extracted to obtain the original pollutant concentration field of the target process; The dynamic meteorological constraint matrix is ​​divided into transport constraint sub-matrices, diffusion constraint sub-matrices, accumulation constraint sub-matrices, and wet deposition constraint sub-matrices according to the mechanism of action of meteorological factors. Based on the transmission constraint submatrix, the diffusion constraint submatrix, the cumulative constraint submatrix, and the wet settlement constraint submatrix, the transmission contribution, diffusion contribution, cumulative contribution, and wet settlement reduction of the grid points on the reference grid are calculated respectively. The emission contribution concentration value of the target process is obtained by subtracting the weighted sum of the transport contribution, the diffusion contribution, the cumulative contribution, and the wet deposition reduction from the grid point concentration values ​​of the original pollutant concentration field. The emission contribution concentration value is then spatially reorganized to obtain the emission characteristic concentration field of the target process. Based on the comprehensive meteorological constraint value of the dynamic meteorological constraint matrix, the species emissions of the dynamic emission inventory are divided into a meteorological interference layer and an intrinsic layer. The emissions of the intrinsic layer are arranged according to the spatial organization of the baseline grid to obtain the net source inventory of the target process.

8. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, The process involves constructing a spatiotemporal correlation matrix based on the dynamic emission inventory and the emission characteristic concentration field, performing spatiotemporal convolution processing on the spatiotemporal correlation matrix to obtain the grid-point spatiotemporal convolution response values ​​of the target process, and performing receptor-source analysis on the source grid data corresponding to the grid-point spatiotemporal convolution response values ​​to obtain the collaborative pollution source tracing results of the target process, including: The gridded emissions of pollutant species in the dynamic emission inventory are expanded along the spatial dimension to obtain the emission space vector of the target process. The gridded concentration values ​​of PM2.5 and O3 in the emission characteristic concentration field are expanded along the spatial dimension to obtain the concentration space vector of the target process. The emission spatial vector and the concentration spatial vector are stacked in time series to obtain the emission time series spatial matrix and the concentration time series spatial matrix of the target process; the spatial cross-correlation matrix between the emission time series spatial matrix and the concentration time series spatial matrix is ​​calculated, and the spatial cross-correlation matrix is ​​extended along the time step direction in the time dimension to obtain the spatiotemporal correlation matrix of the target process; The neighborhood association strength of the grid points in the spatiotemporal correlation matrix is ​​causally convolved in the time step direction to obtain the grid spatiotemporal convolution response value of the target process. Based on the distribution of the spatiotemporal convolution response values ​​of the grid points, the grid points whose response values ​​exceed the target and a set verification threshold are identified as pollution source grid points. Based on the location information of the pollution source grid points, the concentration time series of the pollutant species are extracted from the emission characteristic concentration field. The concentration time series is subjected to non-negative matrix decomposition to obtain the source contribution matrix and source spectrum matrix of the target process. Based on the source contribution matrix and the source spectrum matrix, the key precursor categories and pollution contribution of the target process are determined. The PM2.5 to O3 concentration ratio, temperature parameters, and radiation parameters of the pollution source grid points are obtained. Based on the temperature parameters and radiation parameters, the control type threshold of the target process is determined. Based on the control type threshold, the concentration ratio is judged to obtain the sensitivity type of the target process. The location and response value intensity of the pollution source grid points, the category of the key precursor, the pollution contribution, and the sensitivity type are combined to form the collaborative pollution source tracing result of the target process.

9. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 1, characterized in that, The step of overlaying the collaborative pollution source tracing results with the dynamic emission inventory in a grid-like manner to obtain the source tracing business results of the target process includes: The response values ​​of the source tracing grid points in the collaborative pollution source tracing results are matched with the pollutant species emission amounts of the same grid point in the dynamic emission inventory to obtain the pollution source matching table for the target process; According to the pollution source matching table, the pollutant proportion set of the source tracing grid points is extracted from the dynamic emission inventory, and the precursor category corresponding to the pollutant proportion set is used as the key precursor identifier of the target process; The location information of the source tracing grid points is spatially overlaid with the emission source category information in the dynamic emission inventory to identify the emission source category to which the source tracing grid points belong. The proportion of the response value corresponding to the emission source category in the total response value of the source tracing grid is taken as the pollution contribution of the target process; The location information of the source tracing grid points, the identification of the key precursors, the emission source category, and the pollution contribution are combined to form the source tracing business results of the target process.

10. The method for precise source tracing of PM2.5 and O3 synergistic pollution based on multi-source heterogeneous data fusion as described in claim 9, characterized in that, The step of using the proportion of the response value corresponding to the emission source category in the total response value of the source tracing grid as the pollution contribution of the target process includes: The response values ​​of the source grid points are accumulated over time steps to obtain the cumulative response intensity of the target process; The grid points on the reference grid in the target process are sorted according to the cumulative response intensity to obtain the cumulative response intensity sequence of the target process, and the grid points in the cumulative response intensity sequence whose response intensity reaches the cumulative contribution threshold are used as the hot spot screening grid point set of the target process. Spatially match the location of the hotspot screening grid points with the emission source categories marked in the dynamic emission inventory to obtain the actual emission source identifier of the target process; Extract the component emissions corresponding to the actual emission source identification from the dynamic emission inventory, and calculate the contribution value of the hot spot screening grid points to the common precursors of PM2.5 secondary formation and O3 photochemical pollution based on the component emissions; The proportion of the contribution value corresponding to the emission source category in the contribution value of the common precursors is used as the pollution contribution of the target process.

Images