Methods and Systems for Locating and Analyzing Atmospheric Pollution Sources Based on Monitoring Stations

Abstract

This invention discloses a method and system for locating and analyzing atmospheric pollution sources based on monitoring stations, belonging to the field of environmental monitoring and pollution source apportionment. The method involves acquiring pollutant concentration sequences and spatial coordinates of monitoring stations; extracting the main disturbance characteristic parameters of each station to construct a disturbance time-series vector and a propagation delay matrix; generating a candidate vector set of pollution propagation paths based on the relationship between delay and space, and performing path consistency assessment to select the path with the highest coherence; combining meteorological data to infer the initial location of the pollution source; establishing a local anti-diffusion model based on the estimation results; reconstructing concentration data from neighboring stations; and verifying the effectiveness of source location through error feedback. This invention can achieve high-precision identification of pollution source locations in complex pollution scenarios, possessing technical advantages such as strong timeliness, high adaptability, and verifiable location results.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring and pollution source apportionment, specifically to a method and system for locating and analyzing atmospheric pollution sources based on monitoring stations. Background Technology

[0002] The rapid diffusion and uncertainty of air pollutants have placed higher demands on urban air quality management, especially in areas such as the edges of urban built-up areas, areas with complex terrain, or areas with significantly different local micro-meteorological conditions. Air pollution sources are often difficult to track and locate in a timely manner, which has become a key challenge in environmental law enforcement and emergency response.

[0003] Current methods for locating air pollution sources largely rely on diffusion models to infer pollution pathways or on analyzing source points by retrieving pollutant concentration gradients. However, in areas with high building density and irregular wind field disturbances (such as urban canyons, under viaducts, and areas with abrupt topographic changes), the accuracy of traditional models decreases significantly, and they require a large amount of prior meteorological data, making them unsuitable for the real-time location of source points during sudden pollution events.

[0004] Especially when monitoring stations are sparse and irregularly distributed, the propagation trajectory of pollutants is easily affected by local wind fields and building disturbances, resulting in enhanced nonlinearity of pollution paths and a sharp increase in errors from traditional uniform grid interpolation methods. In addition, some pollution sources release pollutants intermittently, exhibiting strong temporal variability and spatial randomness, rendering methods based on the steady-state diffusion assumption ineffective. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for locating and analyzing atmospheric pollution sources based on monitoring stations, in order to address the shortcomings of the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for locating and analyzing atmospheric pollution sources based on monitoring stations, comprising:

[0007] S100. Obtain the set of atmospheric pollutant concentration change sequences of each monitoring station in the target area within a preset time period, C={C1,C2,...,Ci,...,Cn}, and the set of spatial coordinates of the corresponding stations, P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi.

[0008] S200. Perform multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C for each station, and extract its main disturbance characteristic parameter set Fi, including peak change amplitude Vi, change time Ti and disturbance duration Di.

[0009] S300. Based on Ti in Fi, establish the perturbation time-series vector T={T1,T2,...,Tn} and construct the perturbation propagation delay matrix;

[0010] S400. Based on the spatial coordinate differences between sites and the disturbance propagation delay matrix, construct a candidate vector set of pollution propagation paths L={L1,L2,...,Lk}, where each Lk includes a possible pollution propagation path and its estimated propagation rate.

[0011] S500. Perform path consistency evaluation on each candidate path Lk, calculate the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and select the path with the maximum coordination L{max}.

[0012] S600. Based on the coupling results of the temporal characteristics of each node in L{max} and meteorological parameters, the possible initial location S{est} of the pollution source is inferred.

[0013] S700. Based on S{est}, construct a local anti-diffusion model and reconstruct and simulate the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, output the location of the pollution source; otherwise, return to step S500, update the candidate vector set L of the pollution propagation path, and re-evaluate.

[0014] Preferably, S200 specifically includes:

[0015] S201. Perform wavelet multi-scale decomposition on the pollutant concentration time series Ci, extract high-frequency components and construct perturbation response curves;

[0016] S202. Identify the maximum point of the first derivative of the concentration gradient in the disturbance response curve and determine the abrupt change time Ti;

[0017] S203. Using Ti as the center, search for continuous time intervals where the concentration change exceeds a set threshold to obtain the duration of the disturbance, Di.

[0018] S204. Calculate the maximum amplitude of concentration change within the time interval, which is defined as the peak abrupt change amplitude Vi, and then construct the set of main perturbation characteristic parameters Fi={Ti,Di,Vi}.

[0019] Preferably, S300 specifically includes:

[0020] S301. Based on the abrupt change time Ti in the set of main disturbance characteristic parameters of each monitoring station, construct a disturbance time series vector T={T1,T2,...,Tn}, which represents the time when the pollution disturbance occurs at each station;

[0021] S302. For any two stations Pi and Pj, calculate the absolute value of their abrupt change time difference, and construct the delay unit in the disturbance propagation delay matrix. ;

[0022] S303. Combining the spatial coordinate set P={P1,P2,...,Pn} of the stations, calculate the propagation path vector and theoretical propagation rate for the coordinate difference corresponding to any ΔT{ij};

[0023] S304. Under the condition of satisfying the maximum propagation rate threshold, select effective delay pairs to construct the perturbation propagation delay matrix ΔT.

[0024] Preferably, S400 specifically includes:

[0025] S401. Extract all effective delay units that satisfy the propagation rate constraint from the disturbance propagation delay matrix ΔT, and construct the corresponding site pair set.

[0026] S402. Based on the spatial coordinate difference of the stations and the time sequence of the abrupt change, sort the effective station pairs in time sequence, construct several path sequences that satisfy the disturbance propagation logic, and obtain the pollution propagation path candidate vector set L={L1,L2,...,Lk};

[0027] S403. For each candidate path Lk, calculate the local propagation rate segment by segment based on the spatial distance and time delay between adjacent stations in the path, and estimate the average propagation rate of the entire path.

[0028] Preferably, S500 specifically includes:

[0029] S501. For all stations in each candidate path vector Lk, extract their main perturbation feature parameter set Fi={Ti,Di,Vi} and construct the path perturbation feature matrix.

[0030] S502. Calculate the disturbance time consistency index based on the difference in disturbance duration and abrupt change time interval between adjacent stations in the path.

[0031] S503. Normalize the peak change amplitude of each station in the path, evaluate the amplitude synchronicity, and form a comprehensive coordination score γk by combining the time consistency index.

[0032] S504. Select the path vector L{max} with the highest synergy score from all candidate paths and use it as the priority path input for pollution source localization and reverse inference.

[0033] Preferably, S600 specifically includes:

[0034] S601. Obtain the mutation time Ti and spatial coordinate Pi corresponding to each monitoring station in the path with the maximum synergy L{max}, and construct the path spatiotemporal sequence.

[0035] S602. Synchronously acquire ground meteorological data within the corresponding time period, including wind speed, wind direction and temperature, and perform meteorological interpolation processing on the path sequence according to the time step;

[0036] S603. Starting from the end point of the path, the pollutant diffusion trajectory is constructed in reverse by combining meteorological data, and then gradually traced back to the possible area where the pollution mutation first occurred.

[0037] S604. Based on the dense area of ​​trajectory intersections or the area with the smallest back-calculation error in the propagation rate, determine the possible initial location S{est} of the pollution source.

[0038] Preferably, S700 specifically includes:

[0039] S701. Using the possible initial location of the pollution source S{est} as the initial point, and combining meteorological data and geographical boundary conditions within the corresponding time period, a local anti-diffusion model is constructed.

[0040] S702. The Gaussian anti-diffusion kernel function is used to simulate the spatiotemporal propagation process of pollutants in a local area to generate a concentration estimation field within the target time period.

[0041] S703. Extract the simulated concentration values ​​of the corresponding nearby monitoring station locations from the concentration estimation field and compare them one by one with the actual observed concentration data.

[0042] S704. Calculate the mean square error based on the concentration estimation error of each site. If the simulation error is lower than the standard threshold, then confirm S{est} as the effective pollution source estimation location and output the final pollution source location result; otherwise, return to step S500, update the pollution propagation path candidate vector set L and re-evaluate.

[0043] This invention also provides an atmospheric pollution source location analysis system based on monitoring stations, comprising:

[0044] The monitoring data acquisition module acquires the set of atmospheric pollutant concentration change sequences C={C1,C2,...,Ci,...,Cn} for each monitoring station in the target area within a preset time period, and the set of spatial coordinates of the corresponding stations P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi;

[0045] The disturbance feature extraction module performs multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C at each station, and extracts its main disturbance feature parameter set Fi, including peak abrupt amplitude Vi, abrupt time Ti, and disturbance duration Di.

[0046] The propagation delay analysis module establishes the perturbation time series vector T={T1,T2,...,Tn} based on Ti in Fi, and constructs the perturbation propagation delay matrix;

[0047] The pollution propagation path construction module constructs a pollution propagation path candidate vector set L={L1,L2,...,Lk} based on the spatial coordinate difference between sites and the disturbance propagation delay matrix. Each Lk includes a possible pollution propagation path and its estimated propagation rate.

[0048] The path consistency assessment module performs path consistency assessment on each candidate path Lk, calculates the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and selects the path with the maximum coordination L{max}.

[0049] The pollution source inversion and localization module inverts the possible initial location S{est} of the pollution source based on the coupling results of the temporal characteristics of each node in L{max} and meteorological parameters.

[0050] The anti-diffusion simulation verification module, based on S{est}, constructs a local anti-diffusion model and reconstructs the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, the location of the pollution source is output; otherwise, it returns to the path consistency assessment module to update the candidate vector set L of the pollution propagation path and reassess.

[0051] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0052] 1. The atmospheric pollution source location analysis method based on monitoring stations provided by this invention overcomes the problem of inaccurate source identification in complex scenarios such as sparse station distribution, strong wind field disturbances, and intermittent pollution source release by existing methods through a series of steps including multi-scale disturbance feature extraction, propagation delay matrix construction, pollution path candidate vector identification, and synergy screening. By introducing the quantitative analysis of disturbance time series and main disturbance parameters, the accuracy of pollution propagation path discrimination and the response efficiency of pollution source location are effectively improved, exhibiting strong adaptability and real-time performance.

[0053] 2. This invention constructs an anti-diffusion model by coupling the optimal propagation path with dynamic meteorological parameters and introduces a concentration simulation reconstruction and error feedback mechanism to form a closed-loop optimization process for pollution source location. This mechanism not only achieves the physical verifiability of the pollution source location but also has adaptive correction capabilities. When the simulation error exceeds the standard, it can dynamically return to the path evaluation stage for repeated optimization and deduction, significantly improving the credibility and engineering applicability of the pollution source tracing results. Attached Figure Description

[0054] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

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

[0056] Figure 2 This is a flowchart of the system modules of the present invention. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, 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.

[0058] Example 1, please refer to Figure 1 As shown in this embodiment, the atmospheric pollution source location analysis method based on monitoring stations includes:

[0059] S100. Obtain the set of atmospheric pollutant concentration change sequences C={C1,C2,...,Ci,...,Cn} for each monitoring station in the target area within a preset time period, and the set of spatial coordinates of the corresponding stations P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi.

[0060] In this embodiment, the pollutant concentration is the PM2.5 concentration value. The preset time period is a continuous 24 hours, and the concentration value is recorded every 10 minutes, resulting in a concentration time series Ci={c{i,1},c{i,2},...,c{i,144}} containing 144 data points for each station, where c{i,t} represents the PM2.5 concentration value of station Pi at time point t, in units of PM2.5. .

[0061] The spatial coordinates Pi=(xi,yi) are the two-dimensional planar projection coordinates of station Pi in the geographic coordinate system. The Mars coordinate system (GCJ-02) is adopted and converted into planar coordinates in meters through a unified map projection, which facilitates subsequent path vector calculation and propagation delay estimation.

[0062] To ensure data integrity and representativeness, the selected monitoring stations cover high-density traffic areas, residential areas, industrial parks, and green buffer zones within the target area, exhibiting typical spatial heterogeneity. If a station has more than three consecutive missing data points, the corresponding sequence Ci for that station is removed and not included in subsequent calculations.

[0063] S200. Perform multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C for each station, and extract the set of main disturbance characteristic parameters Fi, including peak change amplitude Vi, change time Ti, and disturbance duration Di.

[0064] S201: For the pollutant concentration time series Ci={c{i,1},c{i,2},...,c{i,t},...,c{i,N}} of monitoring station Pi, where N represents the total number of time steps within the observation period, N=144 in this embodiment, and the sampling interval is 10 minutes, with units of... First, the sequence is decomposed into multiple scales using discrete wavelet transform.

[0065] The wavelet basis function selected is the Daubechies fourth-order wavelet, which is used for three-level wavelet decomposition to obtain the corresponding low-frequency approximate components and multi-level high-frequency detail components. All high-frequency detail components are reconstructed to obtain a high-frequency perturbation component sequence, denoted as the perturbation response curve Ri={r{i,1},r{i,2},...,r{i,N}}. This curve reflects the fluctuation characteristics of short-term abrupt changes in pollutant concentration.

[0066] The first-order difference sequence in the time dimension of the disturbance response curve Ri is calculated, which represents the rate of change of concentration fluctuation at each adjacent point, to identify the maxima of the concentration change rate. Specifically, for any time t, the following steps are performed: The maximum positive value of the first difference is searched throughout the entire sequence, and this time point is denoted as the mutation time Ti, which represents the starting moment when the pollutant concentration increases significantly.

[0067] S203: Centered on the abrupt change time Ti, search forward and backward in the time series Ci to find time intervals where the continuous concentration change amplitude is greater than a set threshold. The concentration change amplitude is defined as the absolute value of the concentration difference between adjacent time points, and the threshold is set to... It is set based on local atmospheric fluctuation background concentration experience.

[0068] When the concentration difference at each step within a continuous time period exceeds the threshold and includes Ti, the total duration of that time period is denoted as the perturbation duration Di, in minutes. The final perturbation interval can be expressed as... , where Di = (d1 + d2) is the sampling interval.

[0069] S204: During the aforementioned disturbance duration interval In this study, the difference between the maximum and minimum concentration values ​​within a given time period is extracted from the pollutant concentration time series Ci, and this difference is used as the peak abrupt change amplitude Vi of the pollution disturbance event.

[0070] Therefore, a set of main perturbation characteristic parameters Fi={Ti,Di,Vi} for monitoring station Pi is constructed to describe the spatiotemporal characteristics of abrupt events caused by pollutants at the station.

[0071] S300. Based on Ti in Fi, establish the perturbation time series vector T={T1,T2,...,Tn} and construct the perturbation propagation delay matrix.

[0072] S301: First, organize the times of abrupt changes in pollutant concentrations at all monitoring stations into an ordered vector according to their spatial numbers, constructing a disturbance time series vector T={T1,T2,...,Tn}, where Ti represents the time point at which station Pi detected the pollutant disturbance, in minutes, with the time starting from a unified reference time (e.g., 0:00 of the current day). This vector is used to represent the chronological order of the pollution disturbances at each station, reflecting the initial propagation trend of the pollution event.

[0073] S302: For any two monitoring stations Pi and Pj, calculate the absolute value of the time difference between their corresponding abrupt change times Ti and Tj, which is defined as the disturbance propagation delay time between the station pair. The formula for calculating the delay unit is: Delay Unit The unit is minutes. This delay represents the time difference it takes for a pollution disturbance to propagate from one site to another, and constitutes an element in the disturbance propagation delay matrix. For all site combinations, the disturbance propagation delay matrix for n is calculated sequentially.

[0074] S303: Based on the disturbance propagation delay unit ΔT{ij} obtained in step S302, and combined with the spatial coordinate set P={P1,P2,...,Pn} of the corresponding stations, where each coordinate Pi=(xi,yi), calculate the Euclidean distance D{ij} between the station pairs. Then, divide this distance by the propagation delay time to obtain the theoretical propagation rate S{ij} of the corresponding station pair. The calculation method is: Theoretical Propagation Rate This rate reflects the speed at which pollution disturbances propagate spatially between the site pairs and is used for subsequent validity assessments.

[0075] S304: To exclude non-genuine disturbance propagation paths and abnormal disturbance responses, a theoretical maximum pollution propagation rate threshold S{max} is set. In this embodiment, the value of S{max} is 500 m / min, which is obtained by converting the typical wind speed (about 8 m / s) and includes appropriate redundancy to adapt to local wind field variations. When S{ij} is less than or equal to S{max}, the propagation path is determined to be a valid path, and the corresponding delay unit ΔT{ij} is retained; otherwise, it is regarded as an abnormal response or a non-direct propagation path and is excluded. Finally, all delay units that meet the rate constraint conditions are retained to form a corrected disturbance propagation delay matrix ΔT, which serves as the core input data for pollution path identification and consistency evaluation.

[0076] S400. Based on the spatial coordinate differences between stations and the disturbance propagation delay matrix, a candidate vector set L = {L1, L2,..., Lk} of pollution propagation paths is constructed. Each Lk includes a possible pollution propagation path and its propagation rate estimate.

[0077] S401: First, all elements ΔT{ij} in the disturbance propagation delay matrix ΔT are traversed to identify all pairs of stations that satisfy the condition that the propagation rate does not exceed the set threshold S{max}. The propagation rate is calculated as: the Euclidean spatial distance between stations Pi and Pj divided by the corresponding time delay ΔT{ij}. If the result is less than or equal to the maximum propagation rate S{max}, the path is considered physically feasible.

[0078] Let the maximum propagation rate S{max} be 500 m / min, which is set according to the typical wind speed experience value of the city. For all pairs of stations (Pi, Pj) that meet the conditions, an effective pair of stations set is constructed to provide a data basis for subsequent path sorting.

[0079] S402: For each pair of stations in the effective pair of stations set, based on the magnitude relationship between the pollution disturbance mutation times Ti and Tj, the pairs of stations are sorted in chronological order according to Ti < Tj to ensure that the pollutant propagation direction conforms to the time causality. On the basis of the sorting, according to the topological connection principle, the effective pairs of stations sharing the same front and rear stations are connected into a path sequence. For example, if (Pi, Pj) and (Pj, Pk) both belong to the effective pair of stations set E and Ti < Tj < Tk, a three-segment propagation path Lk = [Pi → Pj → Pk] can be constructed. In this way, all possible multi-segment pollution propagation paths are constructed to obtain a candidate vector set L = {L1, L2,..., Lk} of pollution propagation paths. Each vector Lk represents a pollution disturbance propagation sequence from the starting station to the target station.

[0080] S403: For each path Lk in the path candidate vector set L, calculate the local propagation rate between each adjacent station segment in the path. Specifically, for any consecutive station pair (Pm, Pn), calculate the spatial distance D{mn} and time delay ΔT{mn} between the two points, then divide the distance by the time difference to obtain the local propagation rate. .

[0081] The arithmetic mean of the propagation rates of all segments along the entire path Lk is defined as the average propagation rate Sk of the path, which characterizes the overall pollution diffusion trend and speed level of the path. This ultimately forms a path-rate correlation set {(Lk,Sk)}, providing input for subsequent pollution path screening and pollution source location inversion.

[0082] S500. Perform path consistency evaluation on each candidate path Lk, calculate the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and select the path with the maximum coordination L{max}.

[0083] In this embodiment, based on the pollution propagation path candidate vector set L={L1,L2,...,Lk} constructed in step S400, a path consistency assessment is performed on each candidate path, and the temporal and amplitude coordination relationship of pollution disturbances at each monitoring station within the path is quantitatively analyzed, thereby selecting the most representative propagation path.

[0084] S501: For all monitoring stations Pi contained in the path candidate vector Lk, extract their main disturbance feature parameter set Fi={Ti,Di,Vi}, organize the Fi of all stations in the path into a three-column matrix to form the path disturbance feature matrix Mk, which is used for subsequent calculation of the coordination index.

[0085] S502: For each pair of adjacent stations (Pi, Pj) in path Lk, calculate the following two indicators: Disturbance duration difference: defined as... Mutation time interval: defined as .

[0086] Suppose there are m stations in the path, then the number of adjacent station pairs is m-1. Summing all differences and averaging them yields the average disturbance time deviation value deltat(k) for the path. This value is then normalized and used as the disturbance time consistency index Ut(k), which is calculated as follows: Where thetat is the maximum allowable threshold for time difference, and its value is set according to the actual pollution response delay characteristics. In this embodiment, it is set to 30 minutes. The closer Ut(k) is to 1, the more consistent the disturbance time of the stations in the path.

[0087] S503: Normalize the peak mutation amplitude Vi of each site in path Lk. The method is to subtract the minimum amplitude value within the path from each Vi and then divide by the difference between the maximum and minimum amplitude values within the path to obtain the normalized value. . Use the standard deviation of V'i as the measurement index for the degree of amplitude change. The amplitude synchronization index Uv(k) is defined as: Amplitude synchronization ; where σv represents the standard deviation of the normalized amplitude . The smaller the value, the more consistent the amplitude. The closer Uv(k) is to 1, the stronger the amplitude synchronization.

[0088] The final synergy score γk is composed of the time consistency index and the amplitude synchronization index through weighted combination. The calculation method is: Synergy score ; where the weighting coefficients α and β satisfy α + β = 1. In this embodiment, α = 0.6 and β = 0.4 are set to emphasize the priority of time causality.

[0089] S504: Calculate the corresponding synergy score γk for all candidate paths Lk respectively. Define the path vector with the highest score as the optimal path L{max}. This path is used as the main reference path for the reverse deduction of pollution source location and is used for subsequent pollution source estimation.

[0090] S600. Based on the coupling result of the time series characteristics of each node in L{max} and meteorological parameters, reverse-deduce the possible initial position S{est} of the pollution source.

[0091] S601: Extract all monitoring stations Pi in the path vector L{max}, obtain their corresponding pollution mutation times Ti and spatial coordinates Pi = (xi, yi), and arrange them in ascending order of the mutation time Ti to construct a path spatio-temporal sequence: Q = {(P1, T1), (P2, T2),..., (Pm, Tm)}, where m is the number of stations in the path, satisfying T1 < T2 <... < Tm. This path sequence describes the propagation process of pollution perturbation in space and time and serves as the basic input for reverse deduction.

[0092] S602: Obtain the meteorological observation data within the time period corresponding to the path sequence Q, including the ground wind speed V{wind} (unit: m / s), wind direction θ{wind} (unit: degree), and air temperature T{air} (unit: °C). The time resolution of the data is 10 minutes. Align the meteorological data with Ti in the path sequence at each time step, and use linear interpolation to perform time alignment and spatial local smoothing processing on the meteorological parameters so that each path node (Pi, Ti) has a corresponding meteorological parameter triple: Mi = {V{wind}(Ti), θ{wind}(Ti), T{air}(Ti)}; this interpolation process is used to ensure the continuity and timeliness of the meteorological input for pollution trajectory deduction.

[0093] S603: Starting from the end point of the path, Pm, and combining its abrupt change time Tm and corresponding meteorological parameter Mm, the meteorological-guided reverse trajectory method is used to gradually trace the diffusion path backward.

[0094] Within each time step Δt, the possible upward propagation distance and direction of the pollutants are calculated based on wind speed and direction. The spatial displacement vec{Δx} at each step of the reverse path is calculated as follows: Where vec{d}(θ{wind}) is the unit wind direction vector. This reverse displacement is accumulated to the path point of the previous moment to obtain the progressive spatial trajectory of the possible sources of pollutants. Continuing to backtrack to the earliest path node P1, a complete pollution back-diffusion trajectory 'back' is obtained, representing the reverse path of the possible directions of pollutant propagation.

[0095] S604: Perform spatial overlap analysis on multiple back-diffusion trajectories to identify the intersection areas of multiple paths during the backtracking process, and record the area with the highest density of trajectory intersections as the hot zone where pollution sources may occur.

[0096] Meanwhile, the error between the propagation time from each candidate source point to the actual observation point and the predicted time is calculated in the reverse path. The least squares method is used to evaluate the propagation rate error corresponding to each source point, and the region with the smallest propagation error is selected as a supplementary criterion.

[0097] Taking into account both trajectory intersection density and velocity back-calculation error, the initial location S{est} for pollution source estimation was finally determined as the input location for subsequent local diffusion simulation and pollution source confirmation.

[0098] S700. Based on S{est}, construct a local anti-diffusion model and reconstruct and simulate the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, output the location of the pollution source; otherwise, return to step S500, update the candidate vector set L of the pollution propagation path, and re-evaluate.

[0099] S701: Using the two-dimensional coordinates of S{est} as the model source point, extract the corresponding abrupt change time T{est}, and simultaneously acquire meteorological data within this time period, including wind speed V{wind} (unit: m / s), wind direction θ{wind} (unit: angle), and atmospheric stability level. The spatial scope is set as a circular area centered on S{est} with a radius of 2000 meters, and the boundary is processed using a flux-free condition.

[0100] When constructing the local back-diffusion model, a two-dimensional Gaussian back-diffusion kernel function is used to describe the upstream propagation distribution of pollutants in a non-uniform wind field. The lateral and longitudinal diffusion coefficients in the kernel function are determined by looking up a table based on the Pasquill stability level. The wind field is dynamically updated using hourly interpolation to simulate the back-diffusion trajectory and impact range of the pollution.

[0101] S702: Within the aforementioned model area, spatiotemporal discretization is performed with a spatial resolution of 10 meters × 10 meters and a time step of 10 minutes. For each spatial grid point G{ij}, the pollutant concentration value C{ij}(t) is calculated within the target time period [T{est}, T{end}], forming a three-dimensional concentration estimation field C(x,y,t), which represents the dynamic distribution estimation result of pollutants within a local area. In the simulation calculation, the pollution source intensity is set to unit emission (1 microgram / second) for relative concentration distribution simulation. The final simulated concentration value will be normalized and proportionally fitted to the measured value.

[0102] S703: Extract the simulated concentration time series C{i}(t) at the corresponding coordinates of all neighboring monitoring stations from the concentration estimation field, and compare it one-to-one with the actual observed concentration series Ci(t) of each station within the same time period. For each station, calculate the mean square error Ei between its simulated and observed values, using the following formula: Where N represents the number of time steps, and t represents each time point.

[0103] S704: Average the errors Ei of all stations to obtain the overall simulation error E{avg}. Set the standard threshold ε for error judgment, which is 25 in this embodiment (unit: micrograms per cubic meter squared), i.e., the upper limit of the allowable average mean square error.

[0104] If E{avg} is lower than ε, then the current source point estimated location S{est} is considered to be able to effectively reconstruct the concentration field of neighboring sites, confirming it as an effective pollution source location, and the output is the final location result.

[0105] If E{avg} is higher than ε, it indicates that there is a deviation in the current path or source estimation. Return to step S500, re-perform path consistency assessment from the pollution propagation path candidate vector set L, screen a new path with the highest degree of synergy L{max}, and repeat the pollution source back-inference and simulation verification process until the error meets the threshold requirement.

[0106] Example 2, please refer to Figure 2 As shown in this embodiment, the atmospheric pollution source location and analysis system based on monitoring stations includes:

[0107] The monitoring data acquisition module acquires the set of atmospheric pollutant concentration change sequences C={C1,C2,...,Ci,...,Cn} for each monitoring station in the target area within a preset time period, and the set of spatial coordinates of the corresponding stations P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi;

[0108] The disturbance feature extraction module performs multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C at each station, and extracts its main disturbance feature parameter set Fi, including peak abrupt amplitude Vi, abrupt time Ti, and disturbance duration Di.

[0109] The propagation delay analysis module establishes the perturbation time series vector T={T1,T2,...,Tn} based on Ti in Fi, and constructs the perturbation propagation delay matrix;

[0110] The pollution propagation path construction module constructs a pollution propagation path candidate vector set L={L1,L2,...,Lk} based on the spatial coordinate difference between sites and the disturbance propagation delay matrix. Each Lk includes a possible pollution propagation path and its estimated propagation rate.

[0111] The path consistency assessment module performs path consistency assessment on each candidate path Lk, calculates the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and selects the path with the maximum coordination L{max}.

[0112] The pollution source inversion and localization module inverts the possible initial location S{est} of the pollution source based on the coupling results of the temporal characteristics of each node in L{max} and meteorological parameters.

[0113] The anti-diffusion simulation verification module, based on S{est}, constructs a local anti-diffusion model and reconstructs the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, the location of the pollution source is output; otherwise, it returns to the path consistency assessment module to update the candidate vector set L of the pollution propagation path and reassess.

[0114] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for locating and analyzing atmospheric pollution sources based on monitoring stations, characterized in that: include: S100. Obtain the set of atmospheric pollutant concentration change sequences of each monitoring station in the target area within a preset time period, C={C1,C2,...,Ci,...,Cn}, and the set of spatial coordinates of the corresponding stations, P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi. S200. Perform multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C for each station, and extract its main disturbance characteristic parameter set Fi, including peak change amplitude Vi, change time Ti and disturbance duration Di. S300. Based on Ti in Fi, establish the perturbation time-series vector T={T1,T2,...,Tn} and construct the perturbation propagation delay matrix; S400. Based on the spatial coordinate differences between sites and the disturbance propagation delay matrix, construct a candidate vector set of pollution propagation paths L={L1,L2,...,Lk}, where each Lk includes a possible pollution propagation path and its estimated propagation rate. S500. Perform path consistency evaluation on each candidate path Lk, calculate the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and select the path with the maximum coordination L{max}. The S500 specifically includes: S501. For all stations in each candidate path vector Lk, extract their main perturbation feature parameter set Fi={Ti,Di,Vi} and construct the path perturbation feature matrix. S502. Calculate the disturbance time consistency index based on the difference in disturbance duration and abrupt change time interval between adjacent stations in the path. S503. Normalize the peak change amplitude of each station in the path, evaluate the amplitude synchronicity, and form a comprehensive coordination score γk by combining the time consistency index. S504. Select the path vector L{max} with the highest synergy score from all candidate paths and use it as the priority path input for pollution source localization back-inference. S600. Based on the coupling results of the temporal characteristics of each node in L{max} and meteorological parameters, the possible initial location S{est} of the pollution source is inferred. S700. Based on S{est}, construct a local anti-diffusion model and reconstruct and simulate the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, output the location of the pollution source; otherwise, return to step S500, update the candidate vector set L of the pollution propagation path, and re-evaluate.

2. The method for locating and analyzing air pollution sources based on monitoring stations according to claim 1, characterized in that: Specifically, S200 includes: S201. Perform wavelet multi-scale decomposition on the pollutant concentration time series Ci, extract high-frequency components and construct perturbation response curves; S202. Identify the maximum point of the first derivative of the concentration gradient in the disturbance response curve and determine the abrupt change time Ti; S203. Using Ti as the center, search for continuous time intervals where the concentration change exceeds a set threshold to obtain the duration of the disturbance, Di. S204. Calculate the maximum amplitude of concentration change within the time interval, which is defined as the peak abrupt change amplitude Vi, and then construct the set of main perturbation characteristic parameters Fi={Ti,Di,Vi}.

3. The method for locating and analyzing air pollution sources based on monitoring stations according to claim 1, characterized in that: Specifically, S300 includes: S301. Based on the abrupt change time Ti in the set of main disturbance characteristic parameters of each monitoring station, construct a disturbance time series vector T={T1,T2,...,Tn}, which represents the time when the pollution disturbance occurs at each station; S302. For any two stations Pi and Pj, calculate the absolute value of their abrupt change time difference, and construct the delay unit in the disturbance propagation delay matrix. ; S303. Combining the spatial coordinate set P={P1,P2,...,Pn} of the stations, calculate the propagation path vector and theoretical propagation rate for the coordinate difference corresponding to any ΔT{ij}; S304. Under the condition of satisfying the maximum propagation rate threshold, select delay units to construct the perturbation propagation delay matrix ΔT.

4. The method for locating and analyzing atmospheric pollution sources based on monitoring stations according to claim 3, characterized in that: Specifically, S400 includes: S401. Extract all effective delay units that satisfy the propagation rate constraint from the disturbance propagation delay matrix ΔT, and construct the corresponding site pair set. S402. Based on the spatial coordinate difference of the stations and the time sequence of the abrupt change, sort the effective station pairs in time sequence, construct several path sequences that satisfy the disturbance propagation logic, and obtain the pollution propagation path candidate vector set L={L1,L2,...,Lk}; S403. For each candidate path Lk, calculate the local propagation rate segment by segment based on the spatial distance and time delay between adjacent stations in the path, and estimate the average propagation rate of the entire path.

5. The method for locating and analyzing air pollution sources based on monitoring stations according to claim 1, characterized in that: The S600 specifically includes: S601. Obtain the mutation time Ti and spatial coordinate Pi corresponding to each monitoring station in the path with the maximum synergy L{max}, and construct the path spatiotemporal sequence. S602. Synchronously acquire ground meteorological data within the corresponding time period, including wind speed, wind direction and temperature, and perform meteorological interpolation processing on the path sequence according to the time step; S603. Starting from the end point of the path, the pollutant diffusion trajectory is constructed in reverse by combining meteorological data, and then gradually traced back to the possible area where the pollution mutation first occurred. S604. Based on the dense area of ​​trajectory intersections or the area with the smallest back-calculation error in the propagation rate, determine the possible initial location S{est} of the pollution source.

6. The method for locating and analyzing air pollution sources based on monitoring stations according to claim 1, characterized in that: The S700 specifically includes: S701. Using the possible initial location of the pollution source S{est} as the initial point, and combining meteorological data and geographical boundary conditions within the corresponding time period, a local anti-diffusion model is constructed. S702. The Gaussian anti-diffusion kernel function is used to simulate the spatiotemporal propagation process of pollutants in a local area to generate a concentration estimation field within the target time period. S703. Extract the simulated concentration values ​​of the corresponding nearby monitoring station locations from the concentration estimation field and compare them one by one with the actual observed concentration data. S704. Calculate the mean square error based on the concentration estimation error of each site. If the simulation error is lower than the standard threshold, then confirm S{est} as the effective pollution source estimation location and output the final pollution source location result; otherwise, return to step S500, update the pollution propagation path candidate vector set L and re-evaluate.

7. An atmospheric pollution source location analysis system based on monitoring stations, used to implement the atmospheric pollution source location analysis method based on monitoring stations as described in any one of claims 1-6, characterized in that: include: The monitoring data acquisition module acquires the set of atmospheric pollutant concentration change sequences C={C1,C2,...,Ci,...,Cn} for each monitoring station in the target area within a preset time period, and the set of spatial coordinates of the corresponding stations P={P1,P2,...,Pi,...,Pn}, where n is the number of monitoring stations and Ci represents the pollutant concentration time series of station Pi; The disturbance feature extraction module performs multi-scale time-series fluctuation analysis on the set of atmospheric pollutant concentration change sequences C at each station, and extracts its main disturbance feature parameter set Fi, including peak abrupt amplitude Vi, abrupt abrupt time Ti, and disturbance duration Di. The propagation delay analysis module establishes the perturbation time series vector T={T1,T2,...,Tn} based on Ti in Fi, and constructs the perturbation propagation delay matrix; The pollution propagation path construction module constructs a pollution propagation path candidate vector set L={L1,L2,...,Lk} based on the spatial coordinate difference between sites and the disturbance propagation delay matrix. Each Lk includes a possible pollution propagation path and its propagation rate estimate. The path consistency assessment module performs path consistency assessment on each candidate path Lk, calculates the coordination index of the set of main disturbance characteristic parameters Fi of each node in the path, and selects the path with the maximum coordination L{max}. The pollution source inversion and localization module inverts the possible initial location S{est} of the pollution source based on the coupling results of the temporal characteristics of each node in L{max} and meteorological parameters. The anti-diffusion simulation verification module, based on S{est}, constructs a local anti-diffusion model and reconstructs the concentration data of neighboring sites. If the simulation error is lower than the standard threshold, the location of the pollution source is output; otherwise, it returns to the path consistency assessment module to update the candidate vector set L of the pollution propagation path and reassess.

Images