A method and system for quantitatively assessing the contribution of greenhouse gas sources

Abstract

This invention relates to the field of environmental monitoring and atmospheric composition source apportionment technology, specifically a method and system for quantitatively assessing the contribution of greenhouse gas sources. The method includes: extracting background concentrations and calculating anomalies based on hourly greenhouse gas concentration data from various stations using robust local weighted regression; combining anomalies with meteorological data to calculate a comprehensive event score and perform event stratification to identify strong event windows; calculating backward trajectories for strong event hours to generate a normalized spatial weight distribution; dividing observed anomalies into local and external contribution components according to the mass conservation principle based on this distribution; and further determining the source industries affecting observed concentrations based on the diagnostic ratios of characteristic species in the sampled data. This invention enables quantitative and auditable separation of the sources of concentration increments without relying on complex models and strong priors, effectively improving the accuracy, robustness, interpretability, and cross-site comparability of source tracing results.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring and atmospheric composition source apportionment technology, and more specifically to a method and system for quantitatively assessing the contribution of greenhouse gas sources. Background Technology

[0002] In the field of quantitative source tracing of greenhouse gases (GHG) at urban and regional scales, a core challenge has long existed: insufficient coupling between observational data and physical processes. Specifically, hourly data from a single station... Concentration records are superimposed with slow background changes such as seasonal cycles, diurnal cycles, and boundary layer evolution, as well as rapid disturbances such as local emission pulses, upwind transport, and mixing, resulting in significant non-stationary characteristics in the concentration series. Under these conditions, without relying on strong prior assumptions, it is difficult to repeatedly decompose concentration increments into the relative contribution shares of "local sources" and "external transport." Existing technologies mainly revolve around three technical pathways, but all have limitations.

[0003] The first type of approach uses empirical methods such as moving averages, fixed thresholds, or low-pass filtering to distinguish between background and anomalies. This method does not fully consider the robustness of outliers and the time-varying characteristics of background concentration, is prone to introducing seasonal biases, and leads to incomparability of results between different monitoring stations.

[0004] The second type of approach relies on isotope fingerprinting technology, such as using... Isotope endmember discrimination is used. Although this method has strong source identification specificity, the high cost of a single analysis, the cumbersome sample preparation process, and the existence of system uncertainties such as isotope fractionation limit its application in high-frequency, routine monitoring scenarios.

[0005] The third type of pathway is based on chemical transport models or inversion methods, such as source attribution models within the Lagrange or Eulerian framework. While these methods are consistent in terms of physical mechanisms, they are highly sensitive to boundary conditions and emission inventory data, require significant computational resources, have high professional barriers, and the comparability of results across sites and seasons is significantly affected by model parameter settings.

[0006] Furthermore, existing technologies generally lack explicit constraints on the organization of meteorological transport in the event identification stage. For example, they do not systematically consider key factors such as wind direction concentration and wind speed suitability, which can easily lead to misjudging pollutant accumulation under low wind speed conditions or rapid dilution caused by high wind speeds as real transport events. At the same time, the mapping process from "receiver hour" to "upwind grid" lacks transparency and an auditable weight allocation mechanism, resulting in insufficient uncertainty characterization and difficulty in forming a complete and verifiable chain of data evidence.

[0007] Given the uneven distribution of urban monitoring station density and human resources, the aforementioned technical shortcomings are further amplified, severely restricting the ability to provide direct and reproducible quantitative answers to the local and external contributions of concentration increases within a specific time period in environmental regulation and policy evaluation scenarios. Summary of the Invention

[0008] In view of this, the present invention provides a method and system for quantitatively assessing the contribution of greenhouse gas sources. It aims to address the problem in existing greenhouse gas source apportionment technologies that, due to the non-stationarity of observation sequences, the coupling of background and disturbances, and insufficient reliability of event identification, it is difficult to quantitatively allocate local and external contributions to concentration increments in a repeatable and auditable manner without relying on strong prior information or complex models. The present invention constructs an observation-driven and assumption-minimizing methodology to improve the robustness, interpretability, and cross-site and cross-seasonal comparability of greenhouse gas source apportionment results, providing direct and reproducible quantitative evidence for emission regulation and policy evaluation at the city and regional scales.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] A method for quantitatively assessing greenhouse gas source contributions includes the following steps:

[0011] Background extraction and outlier calculation steps: Based on hourly greenhouse gas concentration observation data of the station, the background concentration sequence that changes over time is extracted by a robust local weighted regression method, and the difference between the observed concentration and the background concentration is calculated as the outlier sequence.

[0012] Event scoring and stratification steps: Combine the abnormal quantity sequence with the synchronized wind speed and wind direction meteorological data to calculate the comprehensive event score, and divide each hour into background, weak event or strong event levels based on the percentile of the distribution of the comprehensive event score, and identify the continuous strong event window.

[0013] Spatial weight calculation steps: For each hour within the strong event window, calculate its backward trajectory, and based on the position and height information of the trajectory endpoints, generate a normalized spatial weight distribution for each hour that reflects the potential contribution probability of each upwind region;

[0014] Contribution allocation steps: Based on the spatial weight distribution, define the local influence area, and divide the observed anomalies of each strong event hour into local contribution and external contribution according to the principle of mass conservation.

[0015] In a specific feasible implementation, the robust locally weighted regression method in the background extraction and outlier calculation steps includes:

[0016] For each target time, a temporal neighborhood consisting of neighboring samples of a fixed time length is constructed;

[0017] Initial local weights are assigned to each sample in the neighborhood based on a cubic weighting function.

[0018] A weighted linear fit is performed within the neighborhood to obtain the initial background estimate and residual at the target time;

[0019] Based on the residual, robust weights are calculated using the Tukey dual-weight function, and the robust weights are combined with the initial local weights to form combined weights. The combined weights are then used for iterative refitting to finally obtain a robust background concentration estimate.

[0020] In one specific implementation, the calculation of the comprehensive event score in the event scoring and stratification step further includes:

[0021] Robust standardization based on the median and median absolute deviation was performed on the outliers of carbon dioxide and methane, respectively, to obtain dimensionless standardized scores.

[0022] Calculate the wind direction concentration based on the rolling time window as a bonus item for quantifying the organization of transportation;

[0023] Calculate the penalty for wind speed deviating from the ideal range as a deduction item to eliminate false anomalies;

[0024] The standardized scores of the two gases, the wind direction concentration, and the wind speed penalty are linearly combined to generate the comprehensive event score.

[0025] In one specific implementation, the spatial weight calculation step, generating a normalized spatial weight distribution, further includes:

[0026] Assign an original weight to each endpoint on the backward trajectory. The original weight is the product of an indicator function indicating that the endpoint height is lower than the hybrid layer height, a correction term based on the endpoint height, and an exponential decay term based on the backtrack time.

[0027] The geographic region is gridded, and the original weights of all trajectory endpoints falling into the same grid are summed to obtain the original weights of each grid.

[0028] The original weights of all grids are normalized so that for a single receptor, the sum of the weights of all grids is one, thus obtaining the normalized spatial weight distribution.

[0029] In a specific implementation scheme, the contribution sharing step divides the observed anomaly quantity into a local contribution portion and an external contribution portion by multiplying the observed anomaly quantity of each strong event hour by the sum of the local spatial weights and the sum of the external spatial weights corresponding to that hour.

[0030] In a specific feasible implementation, it also includes:

[0031] Uncertainty quantification step: Confidence interval estimation is performed on the quantitative results of the local contribution and the external contribution by repeated sampling and / or perturbation of key parameters using the bootstrap method.

[0032] In one specific implementation, following the event scoring and stratification steps, the method further includes:

[0033] VOC sample alignment and grouping steps: Align the volatile organic compound samples accumulated over 24 hours with the event stratification results within the sampling day; calculate the proportion of hours within the sampling day that were determined to be strong events; based on a preset proportion threshold, mark the volatile organic compound samples as either the background group or the strong event group.

[0034] In a specific feasible implementation, it also includes:

[0035] Diagnostic ratio analysis and comparison steps: Calculate the diagnostic ratios among the preset multiple volatile organic species in each group of volatile organic samples; statistically compare the distributions of the diagnostic ratios between the background group and the strong event group to identify changes in the compositional characteristics of volatile organic compounds under the influence of strong events.

[0036] A system for quantitatively assessing greenhouse gas source contributions includes:

[0037] The data acquisition and preprocessing module is used to access hourly greenhouse gas concentration observation data and synchronous meteorological data from the site;

[0038] The background extraction and outlier calculation module, connected to the data acquisition and preprocessing module, is configured to extract the background concentration sequence that changes over time based on the greenhouse gas concentration observation data using a robust local weighted regression method, and calculate the difference between the observed concentration and the background concentration as the outlier sequence.

[0039] The event recognition and stratification module, connected to the background extraction and anomaly calculation module, is configured to combine the anomaly sequence with wind speed and wind direction data to calculate a comprehensive event score, and to divide each hour into background, weak event, or strong event levels based on the percentile distribution of the comprehensive event score, thereby identifying a continuous strong event window.

[0040] The spatial tracing and weight calculation module, connected to the event recognition and hierarchical module, is configured to calculate the backward trajectory for each hour within the strong event window, and generate a normalized spatial weight distribution for each hour based on the position and height information of the trajectory endpoints, reflecting the potential contribution probability of each upwind region.

[0041] The contribution allocation and output module, connected to the spatial source tracing and weight calculation module, is configured to define the local influence area based on the spatial weight distribution, and divide the observed anomaly quantity of each strong event hour into a local contribution part and an external contribution part according to the mass conservation principle, and output the quantitative allocation result.

[0042] In a specific feasible implementation, it also includes:

[0043] The uncertainty quantification module, connected to the contribution allocation and output module, is configured to estimate the confidence interval of the quantitative results of the local contribution and the external contribution by means of bootstrap resampling and / or key parameter perturbation;

[0044] The VOC alignment and source parsing module, connected to the event recognition and layering module, is configured as follows:

[0045] Align the volatile organic compound samples collected cumulatively over 24 hours with the event stratification results within the sampling day;

[0046] Based on the proportion of hours within a sampling day that are identified as strong events, the volatile organic compound samples are labeled as either the background group or the strong event group.

[0047] Calculate and compare the diagnostic ratio distributions of volatile organic compound species in the background group and the strong event group to identify changes in source characteristics under the influence of strong events.

[0048] Compared with existing technologies, the method and system for quantitatively assessing greenhouse gas source contributions described in this invention are used for quantitative source tracing of greenhouse gases at the urban and regional scales. By constructing an observation-driven, assumption-minimizing, and mass-conserving methodology for allocating local and external contributions, combined with robust local weighted regression for background extraction, event identification based on meteorological organization constraints, backward trajectory spatial mapping, and bootstrap uncertainty quantification, this invention achieves quantitative and auditable separation of local and external contributions in observed concentration increments. This effectively improves the accuracy, robustness, interpretability, and comparability of source apportionment results across sites and seasons, and has the following beneficial effects:

[0049] 1. Improved accuracy and robustness: By employing robust local weighted regression to extract background over extended time scales, and utilizing standardization and quantile thresholds based on median absolute deviation to screen events, anomaly identification becomes insensitive to outliers and seasonal drift. Combined with explicit constraints on wind direction concentration and wind speed suitability, the identification of false positives or false negatives caused by low wind speed accumulation or high wind speed dilution is reduced, thereby improving the accuracy and reliability of contribution allocation results.

[0050] 2. Enhanced interpretability and auditability: By calculating the backward trajectory of each receptor hour and constructing normalized grid spatial weights, a transparent mapping relationship from the observation time to the upwind potential source region is provided. On this basis, the observed enhancement is divided into local and external shares in strict accordance with the principle of mass conservation, ensuring the closure of "the sum of local contribution and external contribution equals the total enhancement", making the entire analysis process and results traceable and verifiable.

[0051] 3. Superior efficiency and scalability: The method relies only on routine hourly greenhouse gas concentration monitoring data and basic meteorological field data from the stations. It does not require the introduction of complex emission prior inventories or reliance on computationally expensive chemical transport models. It also avoids high-cost, long-cycle isotope analysis, reducing the demand for computing resources and technical barriers. It is suitable for routine operational assessment and application promotion of multi-site, long-term series.

[0052] 4. Good comparability and consistency: By adopting a unified and robust statistical strategy and fixed quantile thresholds (such as the 70th and 90th percentiles) for event stratification throughout the entire process, the analysis results of different monitoring stations, different seasons, and even different years are ensured to have consistent benchmarks and comparability, providing a reliable basis for regional-scale pollution source comparison analysis and long-term trend assessment.

[0053] 5. Complete quantification of uncertainty: By adopting the trajectory endpoint bootstrap method combined with key parameter perturbation, a confidence interval estimate is given for the final contribution allocation result, realizing a quantitative characterization of the uncertainty of the method output. This provides a range-based basis including risk perception for regulatory decision-making and policy risk assessment based on the result, and improves the robustness of the result application. Attached Figure Description

[0054] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0055] Figure 1This is an overall flowchart of a method for quantitatively assessing the contribution of greenhouse gas sources according to the present invention.

[0056] Figure 2 This is a system architecture diagram of an embodiment of the system for quantitatively assessing the contribution of greenhouse gas sources according to the present invention.

[0057] Figure 3 This is a system architecture diagram of Embodiment 2 of the system for quantitatively assessing the contribution of greenhouse gas sources according to the present invention. Detailed Implementation

[0058] The technical solutions in the embodiments of the present invention will be clearly and completely described below. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0059] The present invention provides a method for quantitatively assessing the contribution of greenhouse gas sources. Addressing the technical bottlenecks in greenhouse gas source tracing, such as non-stationary background, event identification distortion, opaque spatial mapping, and lack of uncertainty, it constructs an observation-driven, assumption-minimizing, and mass-conserving method for allocating local and external shares.

[0060] The method first employs robust locally weighted regression, using temporal neighborhood tricube weighting and Tukey re-weighting to extract data over a 28-day timescale. The background concentration sequence was obtained, and its outlier was acquired.

[0061] Subsequently, the abnormal quantities of the two gases were normalized by MAD based on median absolute deviation, and then combined with the wind direction concentration based on a three-hour rolling window and the penalty for wind speed deviation from the ideal range of 3-8 m / s to generate a comprehensive event score.

[0062] Using the 70th and 90th percentiles of the comprehensive event score as quantile thresholds, we screened out strong event periods with strong transmission organization characteristics;

[0063] For the identified strong event hour, its 72-hour upwind trajectory is calculated. By constructing an endpoint weighting function that considers endpoint height and time decay and performing grid normalization, the normalized spatial weight distribution from the receptor moment to the upwind region is obtained.

[0064] Based on this definition, the area centered on the station with a radius of 40 kilometers is defined as the local influence circle, and the observed anomalies are divided into local contributions and external contributions under the mass conservation framework, strictly in accordance with the spatial weight share.

[0065] Finally, the results can be aggregated on a calendar day scale, and the trajectory endpoint bootstrapping method combined with key parameter perturbations can be used to give the confidence interval of the contribution, thereby realizing the quantification of uncertainty and robustness verification of the results.

[0066] The entire process relies on the site's hourly updates. The online monitoring data of the wind field, along with station location information and station air sampling data, form a complete closed-loop data evidence chain from 'time domain background stripping' to 'uncertainty quantification', ensuring the feasibility and auditability of the method described in this invention.

[0067] Accordingly, the present invention also provides a system for quantitatively assessing the contribution of greenhouse gas sources.

[0068] The system implements the aforementioned method and process through modular design. For example... Figure 2 As shown,

[0069] The system includes a data acquisition and preprocessing module, a background extraction and anomaly calculation module, an event recognition and hierarchical module, a spatial tracing and weight calculation module, and a contribution allocation and output module.

[0070] The data acquisition and preprocessing module is used to access site data and perform quality control.

[0071] The background extraction and outlier calculation module is connected to the data acquisition and preprocessing module, and extracts background sequences and outliers through robust local weighted regression.

[0072] The event recognition and layering module is connected to the background extraction and anomaly calculation module. Combining the anomaly sequence with wind speed and wind direction data, it calculates a comprehensive event score and divides each hour into background, weak event, or strong event levels based on the percentile distribution of the comprehensive event score, thus identifying a continuous strong event window.

[0073] The spatial tracing and weight calculation module is connected to the event recognition and hierarchical module. For each hour within a strong event window, it calculates the backward trajectory and generates a normalized spatial weight distribution for each hour based on the position and height information of the trajectory endpoints, reflecting the potential contribution probability of each upwind region.

[0074] The contribution allocation and output module is connected to the spatial source tracing and weight calculation module. Based on the spatial weight distribution, it defines the local influence area and divides the observed anomaly quantity of each strong event hour into the local contribution part and the external contribution part according to the mass conservation principle, and outputs the quantitative allocation result.

[0075] Furthermore, such as Figure 3 As shown, the system also includes two extended functional modules: an uncertainty quantification module and a VOC alignment and source resolution module.

[0076] The uncertainty quantification module is connected to the contribution allocation and output module. Through repeated sampling using the bootstrap method and / or perturbation of key parameters, it estimates the confidence interval of the quantitative results of the local contribution and the external contribution, thereby realizing the quantitative characterization of the uncertainty of the method output.

[0077] The VOC alignment and source apportionment module is connected to the event identification and stratification module to achieve refined source apportionment of volatile organic compounds. This module strictly aligns the volatile organic compound samples accumulated over 24 hours with the event stratification results within the sampling day; based on the proportion of hours within the sampling day that are identified as strong events, the volatile organic compound samples are marked as either background group or strong event group;

[0078] The diagnostic ratio distribution of volatile organic compound species in the background group and the strong event group was calculated and compared to identify changes in source characteristics under the influence of strong events, providing chemical fingerprint evidence for the identification of pollution source industries.

[0079] like Figure 1 As shown, the method for quantitatively assessing the contribution of greenhouse gas sources according to the present invention has the following logical chain: data acquisition (preparing raw input) → step 1 (processing greenhouse gas data) → step 2 (identifying events) → steps 3 and 4 (processing VOC data and correlating with events) → step 5 (in-depth analysis of greenhouse gas relationships) → steps 6, 7, and 8 (spatial source tracing and contribution allocation); wherein:

[0080] Data acquisition is the foundation and starting point of the entire method; based on the online monitoring instruments, the latitude and longitude coordinates of the station can be obtained, as well as... and Online monitoring data of wind speed and direction. Through air sampling and analysis, concentration data of volatile organic compounds can be obtained, including (benzene, toluene, ethylbenzene, m / p-xylene, isopentane, isobutane, n-butane, propane, ethane, acetylene, isoprene, etc.).

[0081] In the processing from step 1 to step 8, the output of each step becomes the input of the next step, including:

[0082] Step 1 produces background sequences and outlier sequences;

[0083] Step 2 uses anomaly sequences and meteorological data to generate a composite event score and a strong event window;

[0084] Step 3 aligns the strong event window with the VOC sampling day to produce grouped results (STR / BG).

[0085] Step 4: Perform diagnostic ratio analysis on the VOC data based on the grouping results;

[0086] Step 5 utilizes the outlier sequences and grouping results to perform robust regression and analyze source characteristics;

[0087] Steps 6, 7, and 8 are for calculating the trajectory and spatial weight of strong event windows, and finally completing the contribution allocation and province attribution.

[0088] Step 1: Baseline extraction and outlier calculation (background extraction and outlier calculation).

[0089] Step 1 extracts background concentration sequences and calculates outliers, targeting and The two types of gases are processed separately, as follows.

[0090] Step 1.1: Construct the temporal neighborhood;

[0091] Selecting an appropriate time neighborhood for each target time t This is used for subsequent local fitting.

[0092] For a certain time neighborhood Calculate its absolute time distance from the target time t:

[0093]

[0094] in:

[0095] : No. A point in time within an hour;

[0096] For the target hour, the target time point of the background baseline needs to be calculated;

[0097] : The time distance between sample i and the target time t, in hours.

[0098] Neighborhood set selection: Selecting neighborhood sets from the entire time series that are relevant to the target time. The shortest time distance The samples form a neighborhood set .

[0099] Step 1.1 involves the following fixed parameters:

[0100] Background Time Scale Days; 28 days is empirically optimal (a critical trade-off between baseline stability and fidelity to true seasonal signals).

[0101] Number of target samples in the neighborhood: The size of the neighborhood set (number of samples), calculated using the formula Nβ = 24 × , =28 days, then use the observation data of 672 hours before and after the target time t to fit the background value at time t;

[0102] The set consisting of the Nβ samples closest to t.

[0103] Maximum distance calculation: Determine the maximum time interval between the sample in the neighborhood and the target time.

[0104] To standardize the normalization scale, the maximum distance of the neighborhood is defined as follows:

[0105]

[0106] in: Indicates the distance from the target point in the neighborhood. The farthest sample and The time difference, i.e., the maximum neighborhood distance;

[0107] : The time distance between sample i and the target time t.

[0108] Step 1.2: Calculate the local weights;

[0109] We assign weights to each sample in the neighborhood so that samples that are closer to the target time in time have a greater impact on the fit.

[0110] To ensure that nearby time steps have greater weights and distant time steps have smaller weights during fitting, a cubic weighting function is introduced:

[0111]

[0112] in:

[0113] : The local weight of sample i at the target time t. That is, at time t. When calculating the baseline, the sample Local weights;

[0114] : The time distance between sample i and the target time t;

[0115] : This represents the time difference between the sample farthest from the target point t in the neighborhood and t, i.e., the maximum distance in the neighborhood.

[0116] The cubic weighting function ensures that within the neighborhood, samples closer to time t have a larger weight (maximum 1), while samples farther away from t have a smaller weight (dropping to 0 at the boundary), thus making the estimation of the background baseline mainly determined by the data near t.

[0117] Step 1.3: Locally weighted linear fitting;

[0118] Weighted linear regression is performed within the neighborhood to obtain the initial background concentration estimate and trend at the target time t.

[0119] In the neighborhood Within, the observed sequence is fitted using a weighted first-order straight line:

[0120]

[0121] in:

[0122] The local weight of sample i at the target time t, i.e., at time t. When calculating the baseline, the sample Local weights;

[0123] The intercept of the fitted line is directly used as the initial background concentration estimate at the target time t.

[0124] The slope of the fitted line represents the trend of background concentration change within the neighborhood time window;

[0125] In time The observed greenhouse gas concentration values.

[0126] Find the optimal value by minimizing the loss function J. and This minimizes the weighted sum of squared residuals, yielding the optimal result. then:

[0127]

[0128] The initial background and residuals can be obtained through fitting. Estimate the initial background concentration at the target time t; The residual for each sample i is the deviation between the observed value and the corresponding value on the fitted line. A large residual may indicate that the sample is an outlier (such as a short-term emission event).

[0129] By performing local weighted linear fitting in step 1.3, the first background line can be obtained, which depicts the slow changes on a 28-day scale, and the deviation of each sample can be calculated to facilitate the subsequent identification of outliers.

[0130] Step 1.4, Robust reweighting (iterative process);

[0131] Short-term spikes may lead to fitting bias. Therefore, iterative methods can be used to reduce the impact of outliers (such as short-term peak values) on the background baseline estimation, improve the robustness of the baseline, and suppress the influence of outliers.

[0132] Step 1.4.1, Robust Scaling Calculation:

[0133]

[0134]

[0135] in:

[0136] : The median of the neighborhood residuals;

[0137] Median absolute deviation measures the degree of dispersion of the residuals.

[0138] : Robust scaling estimate, used as an alternative to standard deviation. The constant 0.6745 is to ensure that S(t) is unbiased relative to the standard deviation under the assumption of normal distribution.

[0139] Step 1.4.1 involves the following fixed parameters:

[0140] The coefficient is 0.6745 (in a normal distribution, MAD / 0.6745 is equivalent to the standard deviation).

[0141] Step 1.4.2: Calculate the robust weights (Tukey double weight function):

[0142]

[0143] in:

[0144] Robust weights corresponding to residuals;

[0145] Cutoff constant, with a fixed value. This value determines how many samples with large residuals are considered outliers and their weights are reduced to zero. The absolute value of the residual is greater than... The samples will be completely excluded from the fit.

[0146] Step 1.4.3: Calculate the combined weights and refit;

[0147]

[0148] Use new combined weights Replace the original local weights and repeat the local weighted linear fitting process in step 1.3, minimizing the new weighted loss function to obtain the updated parameter estimates. ,in Indicates the number of iterations. Iteration After this, the final background baseline value at the target time t is the intercept after the third iteration:

[0149]

[0150] Through multiple iterations, the influence of outliers with persistently large residuals (such as short-lived emission peaks) is gradually weakened, thereby ensuring that the fitted background baseline better reflects slow, regional background changes without being affected by short-term local events.

[0151] Step 1.5: Calculate outliers;

[0152] The concentration deviation at each time point relative to the background baseline is obtained to identify enhancement events or cleaning periods.

[0153] Finally, the original observations are compared with the baseline:

[0154]

[0155] right and We obtained the following respectively:

[0156]

[0157]

[0158] in:

[0159] The concentration at this moment is higher than the background concentration, indicating an enhanced event;

[0160] The concentration at this moment is lower than the background concentration, which represents a clean period.

[0161] The original concentration is decomposed into "background + anomaly". The anomaly quantity is the basis for subsequent event scoring and daily-scale VOC alignment.

[0162] After completing step 1.5, the following two key time series will be obtained, laying the foundation for subsequent analysis (such as event identification and source analysis):

[0163] Background sequence: , , representing the slowly changing regional background concentration.

[0164] Outlier sequence: , , representing the deviation of the concentration from the background at each time point.

[0165] These two key time series will be used in the next step for robust standardization and event scoring.

[0166] Step 2: Event scoring and layering (including strong event window recognition);

[0167] Step 2 aims to label each hour's greenhouse gas monitoring data as "Background (BG)," "Weak Event," or "Strong Event," and identify persistent "strong event windows" to provide a foundation for subsequent VOC (volatile organic compound) alignment and analysis. The entire process comprises six core steps: robust normalization, wind direction stability calculation, wind speed penalty, event comprehensive scoring, quantile stratification, and strong event window identification. The following is a detailed explanation in execution order:

[0168] The input data for step 2 includes:

[0169] Inherited from the previous step: hourly increase in greenhouse gases;

[0170] (Unit: ppm)

[0171] (Unit: ppb)

[0172] Step 2 additionally requires: weather data for the same hour;

[0173] wind speed (Unit: m / s)

[0174] wind direction (Unit: degrees, 0-360°)

[0175] Step 2.1, Robust standardization (converting enhancements of different dimensions into comparable dimensionless fractions).

[0176] Different dimensions and The enhancement (Δ value) is converted into a comparable dimensionless fraction to avoid the influence of outliers.

[0177] right and Perform robust standardization separately to avoid being affected by outliers.

[0178]

[0179]

[0180] Will Substitute them separately and To obtain the corresponding standardized score and .

[0181] in:

[0182] Original enhancement sequence (ΔCO2(t) in ppm, ΔCH4(t) in ppb).

[0183] : The median of the data over the entire period, representing the typical level of concentration.

[0184] The median, which is the absolute deviation from the median, measures the dispersion of the data and is not sensitive to outliers.

[0185] Robust scaling estimate, with a constant of 0.6745 used to ensure that MAD conforms to the standard deviation under a normal distribution. Approximate standard deviation (a fixed constant of 0.6745 ensures consistency with the standard deviation under a normal distribution).

[0186] Standard deviation (only when) (Used as a fallback option).

[0187] Standardized fraction, dimensionless.

[0188] By eliminating dimensional differences, and The enhancements are comparable, and the standardization process is insensitive to extreme values ​​and more robust to outliers.

[0189] Step 2.2, Wind direction stability;

[0190] Quantify the stability of wind direction and identify organized air mass transport processes.

[0191] By discretizing the wind direction and calculating the concentration of the "most important wind direction sector" within a 3-hour rolling window, the organization of the transport process can be quantified.

[0192]

[0193]

[0194]

[0195] in:

[0196] Wind direction (degrees): 0° indicates northerly winds, 90° indicates easterly winds, etc.

[0197] : Wind direction sector numbering, dividing 360° into 8 equal 45° sectors.

[0198] :by The timeframe is a 3-hour window (fixed length) at the end of the timeframe.

[0199] : The window belongs to the first Hours per sector.

[0200] : Wind direction stability, the closer the value is to 1, the more consistent the wind direction.

[0201] A consistently stable wind direction may indicate organized air mass transport (such as industrial emissions). This indicator serves as a bonus for event identification, quantifying "wind direction stability" as a bonus from 0 to 1, highlighting potential events caused by continuous, single-direction wind transport.

[0202] Step 2.3, Wind speed penalty (deducting false anomalies caused by extremely low / high wind speeds);

[0203] Exclude false anomalies caused by excessively low wind speeds (pollutant accumulation) or excessively high wind speeds (rapid dilution).

[0204] Linear deductions are applied to wind speeds that deviate from the "ideal range (wind speed 3-8 m / s)"; no deductions are applied to wind speeds within the range.

[0205]

[0206] in:

[0207] Wind speed (m / s).

[0208] : Lower threshold; below this value, pollutants may accumulate due to calm winds.

[0209] : Upper limit threshold. Values ​​exceeding this may be diluted due to strong winds.

[0210] : Wind speed penalty (the penalty for wind speed deviating from the ideal range). It is a positive number when the wind speed deviates from the ideal range. The larger the value, the stronger the penalty.

[0211] This eliminates false concentration anomalies caused by excessively low wind speeds (leading to pollutant accumulation) or excessively high wind speeds (leading to rapid dilution), thereby improving the accuracy of event identification.

[0212] Step 2.4, Comprehensive event scoring (combining gas anomalies and meteorological factors);

[0213] By merging gas anomalies and meteorological factors, a unified event likelihood index is generated.

[0214] In short, the robust standardized scores of the two gases are combined with the two meteorological factors in a "plus-minus" pattern to form a scalar composite event score.

[0215]

[0216] in:

[0217] , : and Standardized scores.

[0218] Wind direction stability (bonus points).

[0219] Wind speed penalty (deduction).

[0220] : Overall event score (dimensionless), the higher the value, the more "like a real event" that hour.

[0221] By integrating information from multiple dimensions, bias from a single indicator is avoided. All weights are set to 1 to ensure objectivity, and the four summation terms are dimensionless. The weights are fixed at 1 to prevent subjective parameter adjustments.

[0222] We obtained a unified metric for hourly “event likelihood” that takes into account both greenhouse gas concentration anomalies and meteorological transport conditions.

[0223] Step 2.5, Quantile Stratification (mark each hour as BG / Weak / Strong);

[0224] Throughout the research period Using fixed quantiles as thresholds in the distribution (the 70th and 90th percentiles of the composite event scores calculated in step 2.4 are used as thresholds respectively), the continuous composite event scores are discretized into three categories, and ES(t) is defined. The label is Strong, which will The tag is Weak, which will The label is BG;

[0225]

[0226]

[0227] in:

[0228] : The 70th and 90th percentiles of the composite event score (fixed stratification points).

[0229] : The category label for each hour, with values ​​of BG, Weak, or Strong.

[0230] Transforming continuous composite event scores into discrete categories facilitates subsequent statistical analysis and matching with downstream VOC sampling data.

[0231] Step 2.6: Identify strong event continuous windows (merge consecutive strong hours into event segments);

[0232] Merge consecutive Strong hours to form event segments and filter out transient anomalies. Combine adjacent hours marked as Strong into consecutive time periods and retain only windows with a duration of at least 3 hours to filter out isolated outliers.

[0233] Define the Strong indicator function:

[0234]

[0235] Identify all that meet the requirements Continuous value of 1 and duration The window.

[0236] in:

[0237] :Strong indicator function (1 indicates Strong, 0 otherwise).

[0238] Fixed minimum duration: Hours, minimum window duration, filtering isolated peak values.

[0239] Output: Several strong event windows And its duration in hours.

[0240] Seeking all satisfaction For intervals where each value is 1 and the duration is ≥3 hours, output the start and end times [Tstart, Tend] and duration of each window. Ensure that those identified as "strong events" are continuous processes, not transient spikes, to facilitate meaningful alignment with 24-hour VOC sampling.

[0241] After step 2.6 is completed, the following output will be obtained:

[0242] Output per hour:

[0243] Overall Event Score

[0244] Standardized gas anomalies ,

[0245] meteorological factors ,

[0246] Layered tags

[0247] Strong event window collection:

[0248]

[0249] The output of step 2.6 will be used in step 3 (VOC-Event Alignment) to calculate indicators such as "Percentage of Strong Hours within the Sampling Day" and "Average Enhancement within the Sampling Day", thereby marking 24-hour VOC samples as BG or STR (strong events) and extracting the summary quantities for subsequent regression, ratio and trajectory analysis.

[0250] Step 3: Align and group VOC (sampling date) with event sequence;

[0251] The core objective of step 3 is to rigorously align a daily collected VOC sample with an hourly event sequence, calculate key indicators for this sample, and determine its group (STR or BG), providing a foundation for subsequent comparative analysis and attribution studies. The entire process sequentially includes five core steps: defining the hourly set of the sampling day, calculating strong event statistics, calculating daily averages, group determination, and identifying representative moments. The following is a detailed explanation in execution order:

[0252] Step 3 objective: To strictly align a bottle of VOC sample collected continuously from 00:00 to 24:00 each day with the hourly event sequence obtained in the previous step, calculate key indicators for this bottle of VOC sample and determine its group for subsequent analysis.

[0253] Input data (hourly sequence inherited from the previous step):

[0254] (Unit: ppm)

[0255] (Unit: ppb)

[0256] (Comprehensive event score, dimensionless)

[0257]

[0258] Step 3.1: Hourly aggregation of the sampling day;

[0259] Objective: To determine the corresponding calendar day for each VOC sample and to establish the hourly data range for alignment and analysis. Define the first... The corresponding calendar day for the bottled VOC sample is The set consists of all the hours within that natural day. :

[0260]

[0261] Symbol explanation:

[0262] : No. The calendar date corresponding to the bottle of VOC sample.

[0263] : Hourly time point (time resolution: 1 hour).

[0264] : The set of time points for all hours within this natural day.

[0265] Establish a one-to-one correspondence between 24-hour VOC samples and hourly event data (such as ΔCH4(tj), ΔCO2(tj), ES(tj), layer(tj)) to ensure that subsequent calculations are based on the same time range.

[0266] Step 3.2: Calculation of strong event hours and overlap ratio;

[0267] In step 3.1, within the sampling day, we can know which category (strong, weak, bg) each hour corresponds to. We can then calculate the number of strong hours within that sampling day and quantify the time coverage of the VOC sample affected by "strong events" within the sampling day.

[0268] Hours of strong events (total number of hours marked "Strong" within the day):

[0269]

[0270] An indicator function that takes the value 1 when the condition is met, and 0 otherwise.

[0271] The number of hours on that day that were judged as "Strong".

[0272] Valid hours (total event score for the day) (Defined, non-missing total number of hours)

[0273]

[0274] Strong event overlap ratio (the proportion of strong event hours to valid hours):

[0275]

[0276] Symbol explanation:

[0277] : No. The number of hours in the "Strong" category within the day.

[0278] : The number of hours of valid data for the day (excluding hours with missing or invalid data).

[0279] : The proportion of hours with strong events to the total effective hours, ranging from [0,1]. This proportion directly reflects the "time coverage" of VOC samples affected by strong events.

[0280] Data coverage lower limit (fixed rule): If the number of valid hours on a certain day If the coverage is less than 1 hour, the sample bottle is recorded as "insufficiently covered". The system will still retain its statistic, but will not participate in subsequent significance comparisons and regression analyses to ensure statistical reliability.

[0281] Step 3.3: The daily average of the daily increment and the comprehensive event score;

[0282] Calculate the average of the greenhouse gas anomaly intensity and the composite event score within the sampling day to generate an index representing the overall anomaly level for that day. Within this process, the arithmetic mean of all valid hourly data is calculated to obtain an index representing the average anomaly intensity for that day. Based on the data obtained in the previous step, the sampling day's... and the results obtained in step 1.5 , The results obtained from 2.4 and 2. .

[0283] Daily average increment:

[0284]

[0285] Daily average increment:

[0286]

[0287] Daily average of comprehensive event scores:

[0288]

[0289] , : Represents the day and The average level of abnormal concentration reflects the overall intensity of the abnormality on the sampling day.

[0290] This represents the average level of the comprehensive event score for the day, comprehensively reflecting the overall "event likelihood" of gas anomalies and meteorological conditions.

[0291] Summation is performed only for valid hours (as defined by ES(tj)). A "daily mean" is generated for this 24-hour sample, consistent with the overall anomaly intensity of its sampling day. These daily means will serve as the basis for subsequent group comparisons (STR vs. BG) and regression analyses, such as for calculating the VOC to greenhouse gas ratio or building regression models.

[0292] Step 3.4: Determine the group of VOCs in this bottle (BG / STR);

[0293] Objective: To divide VOC samples into two groups, "affected by strong events" and "background," based on the coverage ratio of strong events.

[0294] Use a fixed threshold High event overlap ratio To determine the final grouping of the sample bottle:

[0295]

[0296] Parameters and Decision Explanation:

[0297] (Fixed): This threshold indicates that if at least approximately 4 hours (4 / 24≈0.17) within the sampling day are determined to be strong events, then the VOC sample is considered "significantly affected by strong events" and is marked as STR.

[0298] The final label (STR or BG) for this sample.

[0299] All 24-hour VOC samples were clearly divided into two groups: "samples affected by strong events (STR)" and "background samples not significantly affected by strong events (BG)" for subsequent grouped statistical and comparative analysis. The threshold selection was based on an empirical balance: if the threshold was too low (e.g., <0.10), weak events might be misclassified as STR, introducing noise; if the threshold was too high (e.g., >0.20), some true events might be missed. 0.15 is a commonly used compromise value to ensure high purity in the STR group.

[0300] Step 3.5, the representative time of the day (used for fixed purposes in diagrams and trajectories);

[0301] Identify sampling day and The hour of peak occurrence is used for subsequent graphical and trajectory analysis. A representative peak hour is defined for each VOC sample within the sampling day (this definition is standardized and used uniformly):

[0302] Peak value represents the moment:

[0303]

[0304] Peak value represents the moment:

[0305]

[0306] , : These represent the times when the CH4 and CO2 enhancements reach their maximum values ​​within the sampling day, respectively.

[0307] argmax: Returns the argument (i.e., the time point) that maximizes the function value.

[0308] When drawing the attached figures, peak hours are used to mark them; in the trajectory analysis stage, the backward trajectory of the STR sample will be calculated hourly for all Strong hours within its sampling day; and the two representative moments mentioned above will be uniformly retained for fixed graphical comparison.

[0309] Step 3.6: Calculate and output the following parameters for direct use later;

[0310] For each bottle of VOC sample (the first one) (Day) Calculate and produce the following data:

[0311] Strong event statistics: (Hours of major events) (Valid hours) (Proportion of strong event overlap).

[0312] Daily average quantity: ( (incremental daily average) ( (incremental daily average) (Daily average of comprehensive event scores).

[0313] Group: .

[0314] Representative moment: ( (peak moment) ( Peak time).

[0315] These outputs will be used directly for:

[0316] Step 1.4 (Diagnostic Ratio Comparison): Compare the VOC samples from the STR and BG groups and analyze their ratio differences. Trajectory Analysis: For the STR group samples, calculate the hourly forward trajectory of all Strong hours within the sampling day, and perform attribution analysis by combining the daily average amount from Step 3.6 and the grouping results.

[0317] Step 3 successfully aligned and grouped VOC samples with high-frequency event data, laying a solid foundation for subsequent in-depth analysis.

[0318] Step 4: Calculation of VOC diagnostic ratio and comparison with BG / STR grouping;

[0319] The core objective of step 4 is to identify specific changes in VOC composition on strong event days by calculating diagnostic ratios among VOC species and conducting systematic statistical comparisons between the BG (background day) and STR (strong event day) groups, thereby providing crucial "chemical fingerprint" evidence for pollution source analysis. The entire process includes four core steps in sequence: diagnostic ratio calculation, group-based summary statistics, significance testing, and physical interpretation.

[0320] Known preconditions (from step 3):

[0321] For each bottle of 24-hour VOC sample (corresponding to the first... (Day) has received:

[0322] (Based on the threshold of strong event proportion) determination).

[0323] (Unit: ppb).

[0324] (Unit: ppm).

[0325] (Effective coverage hours).

[0326] Fixed rule: Only when Only after a certain number of hours will the sample be included in the statistical comparison and testing in step 4 (samples with insufficient coverage will be recorded but not included in the comparison).

[0327] Step 4.1: 8 diagnostic ratios (dimensionless) for each sample bottle.

[0328] Objective: To transform complex VOC spectral information into dimensionless ratios with clear chemical indicative significance, highlighting the characteristics of different pollution sources or atmospheric chemical processes.

[0329] Candidate species and basic data:

[0330] Each 24-hour VOC sample provides the volume fraction of the following species (units: ppt or ppb, units cancel each other out during calculation):

[0331] Aromatic hydrocarbons: benzene Toluene ethylbenzene The sum of m- and p-xylene Alkanes: Isopentane (iC5ᵢ) n-Pentane (nC5ᵢ) Isobutane n-Butane propane ethane ;

[0332] Other: Acetylene Isoprene ;

[0333] (Benzene);

[0334] (Toluene);

[0335] (Ethylbenzene);

[0336] (The sum of m- and p-xylene, (m+p)-Xylene);

[0337] (isopentane);

[0338] (n-pentane);

[0339] (isobutane);

[0340] (n-butane);

[0341] (Propane);

[0342] (Ethane);

[0343] (Acetylene);

[0344] (Isoprene);

[0345] Eight Diagnostic Ratios and Their Chemical Indicative Significance: The table below defines eight key diagnostic ratios and explains their chemical indicative significance. These ratios were selected based on their recognized source tracing role in environmental science. Definition of the Eight Diagnostic Ratios:

[0346] The main chemical indicator is the contribution intensity of traffic emissions (especially gasoline vehicle exhaust and evaporative emissions). An increased ratio usually indicates a stronger impact from motor vehicle sources.

[0347] The main chemical indicator is the degree of atmospheric photochemical aging of aromatic hydrocarbons. A decrease in the ratio indicates a longer residence time of the air mass or stronger photochemical activity (because MPX reacts more readily with OH radicals than EB).

[0348] The main chemical indicator is the contribution of oil and gas evaporation sources (especially gasoline evaporation). An increase in the ratio indicates a significant impact from this source type.

[0349] Its main chemical indication is that it indicates a combustion process (such as biomass combustion) or a specific industrial source.

[0350] Key chemical indications: the contribution intensity and "freshness" of natural gas / liquefied petroleum gas (LPG) sources. An increased ratio indicates a significant recent impact from emissions of that source.

[0351] The main chemical indications further refine the compositional characteristics of natural gas / LPG sources.

[0352] Its main chemical indicative significance is its indication of incomplete combustion processes (such as motor vehicle exhaust and coal combustion). An increased ratio indicates a relatively stronger contribution from the combustion source.

[0353] The main chemical indicator is the contribution of biological sources (vegetation emissions) relative to anthropogenic sources (benzene mainly comes from combustion and solvent use). The ratio typically increases significantly during the daytime in summer.

[0354] Calculation rules and data validity:

[0355] Denominator validity: The ratio is calculated only if the concentration of the denominator species is >0 and not missing; otherwise, the ratio for that sample is recorded as missing.

[0356] No data transformation: The ratio is calculated directly using the original concentration without logarithmic transformation or truncation. Robust statistics (median, quartiles) are then used to naturally reduce the impact of potential outliers.

[0357] Data coverage requirement: only when the sample's effective hours The sample was only included in the statistical comparison in step 4.1 after a certain number of hours.

[0358] Step 4.1 Projects the complex VOC spectral information onto eight specific "chemical axes," for example... Indicating the degree of aging, Indicators such as the freshness of the natural gas source facilitate systematic comparisons between the BG and STR groups.

[0359] Step 4.2: Perform summary statistics by grouping by BG / STR;

[0360] Objective: To quantify the distribution characteristics of the BG and STR groups at each diagnostic ratio and to calculate the differences between the groups.

[0361] For a specified site and a specified diagnostic ratio Perform the following operations:

[0362] Define this ratio in the sample sets of the BG and STR groups:

[0363]

[0364] Calculate the median and quartiles for each group:

[0365] Group BG: Sample size lower quartiles , median Upper quartile .

[0366] STR group: Sample size lower quartiles , median Upper quartile .

[0367] For each group, the following statistics are calculated. These statistics are not sensitive to extreme values ​​and better describe the central tendency and dispersion of the data:

[0368] Median: The middle value after sorting the data by size. If the sample size is odd, the median is used; if the sample size is even, the average of the two middle values ​​is used.

[0369] Quartiles: The lower quartile (P25) is the 25th percentile, meaning that 25% of the values ​​in the data are below this point; the upper quartile (P75) is the 75th percentile. They can be calculated by sorting the data and taking the values ​​at the corresponding positions, or by using linear interpolation.

[0370] The sample size (N{BG} and N{STR}) is the number of observations in each group.

[0371] Quantifying differences between groups:

[0372] The absolute difference (ΔMedian) is the absolute difference of the median: .

[0373] The relative percentage difference (%Δ) is the same as the relative percentage difference of the median. .

[0374] Minimum sample size requirement for each group:

[0375] To ensure the reliability of statistical comparisons, only when and Only when the sample size of any group is insufficient will subsequent significance tests and difference interpretations be conducted. If the sample size of any group is insufficient, the ratio for that site is considered to lack sufficient statistical power for comparison.

[0376] Step 4.3, Significance test: Mann-Whitney U test (two-tailed);

[0377] Purpose of the test: To determine whether the observed difference in median values ​​between the BG and STR groups is statistically significant (i.e., unlikely to be due to chance), and to determine the effect size. , It provides statistical significance evidence to avoid misjudgment due to human factors.

[0378] Test method:

[0379] Merging ranks: (BG Group) and (STR group) All data from the two groups are merged, sorted uniformly from smallest to largest, and assigned ranks (for equal values, the average of their ranks is taken, i.e., the joint rank).

[0380] Calculate the rank sum: Calculate the sum of the ranks of all samples in the STR group, denoted as . .

[0381] Calculate the U statistic: Let the sample size of group BG be... The sample size of the STR group is .

[0382] Calculate the U statistic for the STR group: .

[0383] Calculate the U statistic for group BG: .

[0384] Take the smaller of the two as the final test statistic. : .

[0385] Null hypothesis The distribution of the BG group is the same as that of the STR group.

[0386] in, : Sample size of group BG (i.e., the number of valid samples in this group).

[0387] : Sample size of the STR group.

[0388] : The sum of the ranks of all samples in the STR group in the mixed data.

[0389] Calculate the p-value: In the case of a large sample (usually the sample size of each group is >5), the U statistic approximately follows a normal distribution.

[0390] Calculate the expected value of U: .

[0391] Calculate the variance of U: .

[0392] in , It is the first The number of identical values ​​in a group of parallel pairs.

[0393] Calculate the standardized Z-score: .

[0394] Based on the standard normal distribution table, find the two-sided p-value corresponding to the Z-value.

[0395] The number of samples with the same value in the k-th "tiegroup". For example, if there are 3 identical values ​​in the mixed data, they form a tiegroup. =3.

[0396] T: Parallel correction term, calculated using the following formula: , used to adjust variance calculations (when there are identical values).

[0397] Interpretation threshold (fixed):

[0398] The difference between the two groups is statistically significant (reliable), and the observed difference is unlikely to be accidental.

[0399] The difference is marginally significant, suggesting that a difference may exist, but the strength of the evidence is slightly weak.

[0400] Statistically, the differences between the two groups cannot be distinguished. The observed differences are likely due to random fluctuations.

[0401] Step 4.4, Readings and Physical Direction (provide a directional interpretation for each ratio);

[0402] Transforming statistical differences into directional interpretations of chemical mechanisms for subsequent applications... Regression analysis and interpretation of trajectory-source region linkage:

[0403] The decrease indicates that the aromatic components are more aged (m / p xylene is more easily photochemically consumed than ethylbenzene).

[0404] , The increase indicates that the traffic / solvent emission end-units are "fresher" or that their proportion in the mixed source has increased.

[0405] , , The increase indicates that the end-members of the natural gas / LPG transmission and distribution chain are enhanced and have higher "freshness".

[0406] An increase indicates a relative enhancement in the composition of combustion components.

[0407] An increase (common during summer days) indicates an increased contribution from bio-derived isoprene or a change in its mixing state with gasoline vehicle endmembers.

[0408] Step 4.5: Calculate and output statistical information for direct use later;

[0409] For each "site-ratio" combination Output:

[0410] BG Group Statistics: , , , .

[0411] STR group statistics: , , , .

[0412] Between-group differences: (Absolute difference of median) (Median relative percentage difference).

[0413] Value (only when each of the two groups) (To be given at the time).

[0414] Step 5 Robust regression (divided into BG / STR groups);

[0415] The core objective of step 5 is to perform robust regression analysis on the BG and STR groups using Theil-Sen. Enhancement amount ( )and Enhancement amount ( The proportional relationship is used to determine the impact of strong events. Compared to Whether there has been a systemic change in the emission ratio provides crucial evidence for identifying the dominant emission sources (such as combustion sources, natural gas leaks, etc.). The entire process strictly follows the logical sequence of data preparation → core parameter estimation → model evaluation → robustness verification.

[0416] Input data:

[0417] Hourly increase: (t) (unit: ppb) (t) (unit: ppm) — from step 3.

[0418] Group labels: layer(t)∈{BG,Weak,Strong} (only the BG and Strong groups are used in step 5).

[0419] Methodology (using Theil-Sen robust regression):

[0420] Step 5.1, Data filtering within the group (basic preparation);

[0421] Objective: To prepare a high-quality, effectively paired BG and STR dataset for subsequent regression analysis. Step 5.1 is the foundation for all subsequent calculations.

[0422] For each station and each group G∈{BG,Strong}, select a sample set where both gas data are valid at the same time:

[0423]

[0424] Parameter explanation:

[0425] : The set of all valid data points within group G, used for subsequent regression analysis.

[0426] Filtering logic: Define a valid sample set for each site(s) and each group (G∈{BG,STR}).

[0427] Explanation of parameters and operations:

[0428] Validity criteria: Data at a given time (t) will be included in the set only if both of the following conditions are met simultaneously:

[0429] The group label at this moment is consistent with the group being analyzed.

[0430] The sum at this moment is defined (i.e., not missing, not invalid value).

[0431] Independent screening: The screening of the BG group and the STR group is carried out completely independently to ensure that the data of the two groups are pure and do not mix with each other.

[0432] Output: Two clean datasets are obtained, each containing sample points that belong only to the group and are valid for both gas data.

[0433] Step 5.1 forms the basis of "core parameter estimation." Subsequent slope and intercept calculations, as well as all evaluation and validation steps, are strictly limited to the TG set generated in this step. High-quality data screening is a prerequisite for obtaining reliable regression results.

[0434] Step 5.2, Theil-Sen slope and intercept estimation (core step);

[0435] Objective: To robustly estimate Increase amount The proportional relationship of the increase in magnitude (slope β) and the baseline level (intercept α).

[0436] Let the independent variable be... dependent variable .

[0437] The computational logic and advantages are as follows:

[0438] Pairwise calculation: Calculate the slope of all possible pairs of points in the set (that satisfy AND).

[0439] Median: Sort all calculated slope values ​​in ascending order and take the median as the final slope estimate.

[0440] Robustness stems from the fact that the median is insensitive to outliers. Even if there are individual points with large deviations in the data, the final result is determined by the trend of most "normal" point pairs, thus ensuring the robustness of the estimate.

[0441] Slope calculation formula:

[0442]

[0443] Parameter explanation:

[0444] Theil-Sen slope estimation indicates Enhancement amount and The proportional relationship of enhancement (unit: ppb / ppm). It is robust to outliers by calculating the median of the slope for all points.

[0445] : Time point index, representing different hours.

[0446] : The difference in the amount of enhancement.

[0447] : The difference in the amount of enhancement.

[0448] Intercept estimation:

[0449]

[0450] In obtaining the slope Then, calculate based on all sample points. Then, take the median of these values ​​as the intercept. This is also a robust estimation method.

[0451] Parameter explanation:

[0452] : Intercept, representing when When the enhancement amount is 0 Baseline estimation of enhancement.

[0453] Step 5.3, Goodness of Fit and Residual Evaluation (Model Diagnosis);

[0454] Objective: To evaluate the fit of the Theil-Sen regression line to the actual data and to diagnose the effectiveness of the model. Residual calculation:

[0455] Residual calculation:

[0456]

[0457] Parameter explanation:

[0458] Residual: Residual represents the deviation between the observed value and the fitted value. The residual is the deviation between the observed value y(t) and the predicted value of the regression model, reflecting the part that the model failed to explain.

[0459] Median Absolute Residuals (MAR):

[0460]

[0461] Parameter explanation:

[0462] MAR (Median Absolute Residual): MAR measures the average level of error in the fit. The smaller the MAR value, the smaller the overall deviation between the regression line and the data points, and the higher the fitting accuracy.

[0463] Robust fit goodness of fit:

[0464]

[0465] Parameter explanation:

[0466] This metric is similar to the traditional R², but it is calculated based on the L1 norm (sum of absolute values), making it more robust to outliers. The closer its value is to 1, the stronger the model's ability to explain data variation and the better the fit.

[0467] : The median of the y-values.

[0468] Step 5.3 is a model diagnostic based on the obtained core parameters (α^, β^). It evaluates and supplements the Theil-Sen slope and intercept estimation results, assessing the goodness of fit of the straight line.

[0469] Step 5.4, Robustness verification (result confirmation);

[0470] Backward regression (mutual verification): exchanging independent and dependent variables (assuming...) , Repeat the above steps to obtain another set of slopes. ,intercept , and .

[0471] In an ideal situation, the positive slope (Unit: ppb / ppm) and reverse slope (Unit: ppm / ppb) should be reciprocals of each other.

[0472] If the two are reciprocals of each other, then it means and If the cooperative relationship is stable and symmetrical, then the positive regression results are reliable. Step 5.4 and "core parameter estimation" are parallel calculation processes, but they serve the same goal—verifying robustness. It is an important supplement and verification step for Theil-Sen slope and intercept estimation.

[0473] Confidence interval estimation:

[0474] Randomly draw samples with replacement from the original dataset TG to generate a Bootstrap sample with the same sample size as TG. Repeat this process B=1000 times. For each Bootstrap sample, calculate its Theil-Sen slope estimate. 1000 The values ​​are arranged in ascending order, and their 2.5% and 97.5% quantiles constitute the slope. The 95% confidence interval (CI).

[0475] Relationship with step 5.2: This step quantifies the uncertainty of the Theil-Sen slope and intercept estimation. It relies on the Theil-Sen slope and intercept estimation method but provides crucial information about the stability of the results, thus deepening and expanding upon step 5.2.

[0476] Step 5.4 Core Output:

[0477] For each site and group (BG / STR), the following key results are output for the final inter-group comparison:

[0478] Slope and intercept: , .

[0479] Goodness-of-fit metric: MAR (lower is better). (The closer to 1, the better).

[0480] Robustness verification indicator: Back regression parameters 95% confidence interval for Bootstrap.

[0481] The scientific significance of the results:

[0482] The value approximately represents " and The ratio of enhancement levels. By comparing the BG group and the STR group. The value can be used to determine under strong events. Compared to Whether there has been a systematic change in the proportion of emissions, thereby inferring whether the dominant emission source of strong events has changed (e.g., from combustion sources to natural gas leak sources).

[0483] If STR group Significantly higher than the BG group, which may indicate a higher proportion of days with strong events. Emission sources are activated (e.g., increased contribution from natural gas leaks, because natural gas is primarily composed of...). ).

[0484] If STR group If there is no significant difference or a decrease compared to the BG group, it may suggest that the strong event is dominated by the combustion process (because combustion simultaneously produces...). 2 and And the proportion is relatively fixed.

[0485] Step 6: Calculation of backward trajectory and spatial weights for strong events;

[0486] Input data:

[0487] High-intensity events hourly collection: (Steps 2.5 and 2.6)

[0488] The latitude and longitude coordinates of the station (data acquisition part, which is known).

[0489] Hourly increase: (t) (ppb), (t) (ppm) (Step 1.5).

[0490] Step 6.1: Trajectory Calculation;

[0491] Based on the obtained strong event time points and known station latitude and longitude, trajectory data can be obtained from public trajectory calculation websites simply by inputting specific parameters and time. The following are the specific settings of this invention: for each strong event hour Calculate a 72-hour backward trajectory. The trajectory includes... Each endpoint (time step) (hours). Regarding the first endpoints ( ),Record:

[0492] Position coordinates .

[0493] high (meters, height above ground AGL).

[0494] Rewind duration (Hour).

[0495] If boundary layer height data is unavailable, a fixed blending layer height is used. Meters are used as an approximation.

[0496] Step 6.2: Calculate the original endpoint weights;

[0497] Considering high-altitude shielding, time decay, and boundary layer correction, weights are assigned to each endpoint of the trajectory obtained in the previous step (one trajectory can be obtained at each calculation moment, and each trajectory has 72 endpoints, representing the air mass path tracing back 72 hours from that moment), reflecting its potential contribution to receptor point contamination.

[0498]

[0499] Parameter explanation:

[0500] Hour: Time step.

[0501] : Indicator function, takes the value 1 when the endpoint height is lower than the mixing layer height, otherwise takes the value 0 (masking the high-altitude endpoint).

[0502] Meters: Reference height.

[0503] Meters: Minimum height threshold to avoid division by zero errors.

[0504] : High adjustment index.

[0505] Hour: Time decay constant, which controls the rate at which the weight decays over time.

[0506] The time decay term has a higher weight for endpoints that are closer to the receptor in time.

[0507] This is a correction term based on endpoint height. This is an exponentially decaying term based on the backtracking time.

[0508] Step 6.3: Gridding and Normalization;

[0509] The trajectory movement path and trajectory endpoint weights obtained in the previous two steps are visualized using coordinates. The discrete trajectory endpoint weights are then aggregated onto a regular geographic grid (grid resolution). (Approximately 55km × 55km), facilitating spatial analysis and visualization.

[0510] Original grid weights:

[0511]

[0512] Parameter explanation:

[0513] The original weights of grid g are obtained by summing the weights of all endpoints falling into grid g.

[0514] Grid normalized weights:

[0515]

[0516] Parameter explanation:

[0517] The normalized grid weights, for a fixed t, have: .

[0518] Step 6.4: Calculation of local vs. external share;

[0519] Objective: To quantify the relative importance of strong events being influenced by local sources and external transmission.

[0520] Local impact zone definition: The area within a radius of 40km centered on the station.

[0521] Define local grid collection .

[0522] Define the set of external meshes For all non-local grids.

[0523] Local share:

[0524]

[0525] Parameter explanation:

[0526] Local share: This indicates the weight of the air mass's influence by its local source.

[0527] External share:

[0528]

[0529] Parameter explanation:

[0530] External share: This indicates the weight of an air mass influenced by external transport. Approximately 1: Strong events are mainly affected by local emissions. Approaching 1: Strong events are mainly affected by regional transport.

[0531] Fixed parameter values:

[0532] Trajectory: , Hour.

[0533] Weight calculation: rice, rice, , Hour.

[0534] Grid: .

[0535] Local radius: kilometer.

[0536] Output result:

[0537] Every hour of a major event Grid weights .

[0538] Local share and foreign share .

[0539] By quantifying the physical process of "which geographical regions influence an air mass at a given hour" into a spatial weighted distribution, a foundation is laid for subsequent source region contribution decomposition. Specifically, ΔCH4(t) and ΔCO2(t) are weighted and distributed across each grid cell.

[0540] Step 7: Allocate observational enhancements according to local / external share;

[0541] Input data:

[0542] The weights output in step 6: , (For strong events) ). This represents the probability ratio of an air mass to local and external influences at each hour t of a strong event.

[0543] Observational augmentation: (t) (ppb), (t) (ppm). The actual concentration anomaly relative to the background baseline for each strong event hour, representing the pollution intensity.

[0544] Procedure (distributed hourly and gas-by-gas, following the principle of mass conservation):

[0545] Step 3 may allow us to obtain information about the local and external regions during a strong event, specifically regarding the time period. (t), Contribution ratio of (t) , Therefore, (t), (t) is allocated to the local and external regions to obtain their absolute contribution to this increment, calculated as follows:

[0546]

[0547] Parameter explanation:

[0548] CH4 enhancement amount allocated to the local source.

[0549] : The amount of CH4 enhancement allocated to external transport.

[0550] : The amount of CO2 enhancement allocated to the local source.

[0551] : The amount of CO2 enhancement allocated to external transport.

[0552] Output result:

[0553] Enhancement after decomposition: , , , .

[0554] By quantitatively dividing the total enhancement into "local source contribution" and "external transport contribution" while ensuring the conservation of the observed enhancement, the source contribution is initially separated.

[0555] Step 8: Attributing the external component by province (overlaying emission inventory intensity);

[0556] Step 8 is the refinement stage of source apportionment. It further decomposes the external transport contributions obtained in Step 7 to specific upstream provinces, identifying key source regions that significantly contribute to strong events at the observation points. This achieves precise source tracing "from region to province." Step 6 reveals which provinces influence the trajectory and the magnitude of their influence. Step 8 obtains the absolute contribution values ​​of all external provinces to the observed increment. Based on this, emission inventory data for each province is introduced, which can be obtained from the national emission database. According to the eight characteristic ratios in Step 4.1, the specific industries affected can be determined based on the corresponding ratio ranges. Combining this with the provincial industry emission inventory data, the contribution ratio of these provincial industries to the observed increment can be obtained.

[0557] Input data:

[0558] Step 6.1: External grid weights obtained in the trajectory calculation section (External grid only) ).

[0559] Publicly available data shows the relationship between provincial boundaries and grid affiliation.

[0560] By breaking down provincial emission data into grids, emission inventory intensity can be obtained. , (each grid) of and Emissions).

[0561] practice:

[0562] Calculate the signal strength of a province:

[0563] for :

[0564]

[0565] Parameter explanation:

[0566] : Province p relative to time t Enhanced signal strength, based on trajectory weights and emissions inventory.

[0567] : A grid that belongs to province p and is in the set of external grids.

[0568] for :

[0569]

[0570] Parameter explanation: similar .

[0571] Calculate the contribution ratio of provinces:

[0572] for :

[0573]

[0574] Parameter explanation:

[0575] : The relative contribution of province p to the current hourly increase in external CH4.

[0576] for :

[0577]

[0578] Parameter explanation: similar .

[0579] Increased allocation to provinces:

[0580] for :

[0581]

[0582] Parameter explanation:

[0583] : Foreign elements assigned to province p Increased quantity.

[0584] for :

[0585]

[0586] Parameter explanation: similar .

[0587] Output result:

[0588] Province breakdown of external enhancement: , .

[0589] The above steps allow us to allocate the observational increments caused by the combined effects of other provinces to specific provinces and industries, i.e., the contribution of a particular industry in a specific province to the observational increments.

[0590] This invention achieves verifiability and auditability of results by constructing a complete data processing chain. Specifically, the method sequentially generates and outputs a series of clear intermediate data results: First, robust local weighted regression is used to extract background concentration and anomaly sequences, providing a benchmark for concentration decomposition; then, based on the robust standardized scores of gas anomalies and wind direction stability and wind speed penalty factors, a comprehensive event score is synthesized, and event classification labels and persistent strong event window identifiers are output, providing a quantitative basis for event identification; furthermore, the backward trajectory is calculated for strong events, generating a normalized spatial weight matrix from the receptor time to the upwind grid, and the quantitative ratio of local and external contributions is output accordingly; finally, under the mass conservation framework, the spatial weights are used to divide the observation enhancement into absolute values ​​of local and external contributions, and the uncertainty confidence intervals of key parameters can be given through bootstrapping. The above-mentioned interconnected and traceable intermediate data together constitute a transparent and auditable data evidence chain, ensuring that the entire process from the original observation to the final contribution allocation can be verified and verified.

[0591] The various embodiments described in this specification are presented in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method of quantitatively assessing the contribution of sources of greenhouse gases, characterized in that, The method comprises the following steps: a background extraction and anomaly amount calculation step: based on hourly greenhouse gas concentration observation data of a site, a time-varying background concentration sequence is extracted by a robust local weighted regression method, and the difference between the observation concentration and the background concentration is calculated as an anomaly amount sequence; an event scoring and stratification step: combined with the anomaly amount sequence and synchronous wind speed and wind direction meteorological data, a comprehensive event score is calculated, and each hour is divided into background, weak event or strong event levels based on the distribution percentile of the comprehensive event score, so as to identify a continuous strong event window; a spatial weight calculation step: for each hour in the strong event window, a backward trajectory is calculated, and a normalized spatial weight distribution reflecting the potential contribution probability of each region in the upwind direction is generated for each hour based on the position and height information of the trajectory endpoint; a contribution allocation step: according to the spatial weight distribution, a local influence area is defined, and the observation anomaly amount of each strong event hour is divided into a local contribution part and an external contribution part according to the principle of mass conservation.

2. The method of claim 1, wherein, The robust local weighted regression method in the background extraction and anomaly amount calculation step comprises: a time neighborhood composed of adjacent samples with a fixed time length is constructed for each target time; initial local weights are assigned to each sample in the neighborhood based on a cubic weight function; a weighted linear fitting is performed in the neighborhood to obtain an initial background estimate and a residual for the target time; robust weights are calculated using a Tukey double weight function based on the residual, and the robust weights and the initial local weights are combined to form combined weights, and the combined weights are used for iterative re-fitting to finally obtain a robust background concentration estimate.

3. The method of claim 1, wherein, In the event scoring and stratification step, the calculation of the comprehensive event score further comprises: robust standardization based on the median and the median absolute deviation is performed on the anomaly amounts of carbon dioxide and methane respectively to obtain dimensionless standardized scores; a wind direction concentration degree based on a rolling time window is calculated as a bonus item for quantifying transport organization; a penalty amount for wind speed deviation from an ideal interval is calculated as a deduction item for excluding false anomalies; the standardized scores of the two gases, the wind direction concentration degree and the penalty amount are linearly combined to generate the comprehensive event score.

4. The method of claim 1, wherein, In the spatial weight calculation step, generating the normalized spatial weight distribution further comprises: an original weight is assigned to each endpoint on the backward trajectory, which is the product of an indicator function of the endpoint height being lower than the mixing layer height, a correction term based on the endpoint height and an exponential decay term based on the backtracking time; a geographical area is gridded, and the original weights of all trajectory endpoints falling into the same grid are summed to obtain the original weight of each grid; the original weights of all grids are normalized so that the sum of the weights of all grids for a single receptor hour is one, thereby obtaining the normalized spatial weight distribution.

5. The method of claim 1, wherein, In the contribution apportionment step, the observed anomaly is divided into local contribution and external contribution by multiplying each hour's observed anomaly by the sum of local spatial weight and external spatial weight corresponding to the hour.

6. The method of claim 1, wherein, Further comprising: Uncertainty quantification step: confidence interval estimation of the quantitative results of the local contribution and external contribution through bootstrap resampling and / or key parameter perturbation.

7. The method of claim 1, wherein, After the event scoring and stratification step, further comprising: VOC sample alignment and grouping step: aligning the 24-hour cumulative sampled volatile organic compound sample with the event stratification results within its sampling day; calculating the proportion of hours judged as strong events within the sampling day; based on the preset proportion threshold, labeling the volatile organic compound sample as background group or strong event group.

8. The method of claim 7, wherein, Further comprising: Diagnostic ratio analysis and comparison step: calculating the diagnostic ratio between a plurality of preset volatile organic compound species in each group of volatile organic compound samples; Statistically comparing the diagnostic ratio distribution of the background group and the strong event group to identify the changes in the composition characteristics of volatile organic compounds under the influence of strong events.

9. A system for quantitatively assessing the contribution of sources of greenhouse gases for performing the method of quantitatively assessing the contribution of sources of greenhouse gases according to any one of claims 1 to 8, characterized in that Comprising: A data acquisition and preprocessing module for accessing hourly greenhouse gas concentration observation data and synchronous meteorological data of the site; A background extraction and anomaly calculation module connected to the data acquisition and preprocessing module, configured to extract the time-varying background concentration sequence based on the greenhouse gas concentration observation data through robust local weighted regression method, and calculate the difference between the observed concentration and the background concentration as the anomaly sequence; An event identification and stratification module connected to the background extraction and anomaly calculation module, configured to calculate the comprehensive event score in combination with the anomaly sequence and wind speed and direction data, and divide each hour into background, weak event, or strong event levels based on the distribution percentile of the comprehensive event score, identifying the strong event window; A spatial tracing and weight calculation module connected to the event identification and stratification module, configured to calculate the backward trajectory for each hour within the strong event window, and generate a normalized spatial weight distribution reflecting the potential contribution probability of each region in the upwind direction based on the location and height information of the trajectory endpoint; A contribution apportionment and output module connected to the spatial tracing and weight calculation module, configured to define the local impact area based on the spatial weight distribution, divide each strong event hour's observed anomaly into local contribution and external contribution according to the principle of mass conservation, and output the quantitative apportionment results.

10. The system for quantitatively assessing contributions of sources of greenhouse gases according to claim 9, wherein, Further comprising: An uncertainty quantification module connected to the contribution apportionment and output module, configured to estimate the confidence interval of the quantitative results of the local contribution and external contribution through bootstrap resampling and / or key parameter perturbation; A VOC alignment and source analysis module connected to the event identification and stratification module, configured to: align the 24-hour cumulative sampled volatile organic compound sample with the event stratification results within its sampling day; Based on the proportion of hours determined as strong events within the sampling day, the volatile organic samples are labeled as background group or strong event group; The diagnostic ratio distributions of volatile organic species in the background group and the strong event group are calculated and compared to identify the source signature changes under the influence of strong events.

Images