Method and system for pollution source tracing of chemical industry park based on machine learning

Abstract

This invention relates to a machine learning-based method and system for tracing pollution sources in chemical industrial parks, belonging to the field of pollution source tracing. It utilizes building outline information and real-time wind field data from GIS to perform physical masking filtering and Gaussian plume weight reconstruction on a static distance matrix, generating a dynamic adjacency matrix sequence that characterizes airflow around-the-flow effects and dynamic causal relationships. Based on this, a graph convolutional neural network embedding long short-term memory units is constructed, synchronously inputting standardized sensor spatiotemporal tensors and the dynamic adjacency matrix to extract spatiotemporal causal features in a non-Euclidean space. Finally, by inverting the source probability distribution heatmap and combining it with continuity constraints for coordinate correction, accurate location of pollution sources is achieved under complex interference environments.

Description

Technical Field

[0001] This invention relates to the field of pollution source tracing, specifically to a method and system for pollution source tracing in chemical industrial parks based on machine learning. Background Technology

[0002] With the increasing agglomeration of the chemical industry, chemical industrial parks are home to numerous enterprises with complex processes, leading to a growing risk of emissions of various hazardous chemicals and pollutants. To effectively address potential sudden environmental pollution incidents, the widespread deployment of various gas sensors and the establishment of an efficient pollution source tracing system within the parks have become an industry consensus. Traditional source tracing methods often rely on Gaussian plume models or computational fluid dynamics (CFD) simulations, but these methods typically require precise source parameters and are computationally time-consuming, making it difficult to meet the stringent real-time requirements of emergency response. Therefore, utilizing massive amounts of historical monitoring data to construct data-driven source tracing models based on machine learning, due to their fast response speed and strong nonlinear fitting capabilities, has gradually become a research hotspot in the field of smart environmental protection in chemical industrial parks.

[0003] However, while existing machine learning-based pollution tracing schemes are theoretically feasible, they still face numerous challenges and significant technical shortcomings in the complex application scenarios of chemical industrial parks. First, existing deep learning models, such as standard convolutional neural networks (CNNs), assume that the input data has a regular Euclidean structure and extract features based on fixed geometric distances. However, in chemical industrial parks, the distribution of fixed sensor sites is sparse and unstructured, and this static modeling approach cannot accurately represent the true topological relationships of the sensor network. Second, existing technologies generally ignore the dynamic characteristics and physical constraints of pollutant transport under micrometeorological conditions. Chemical industrial parks are filled with factories, and the airflow field is greatly affected by building obstruction and disturbance, forming an extremely irregular and unsteady micrometeorological environment. The causal relationships between sensors exhibit strong dynamic connectivity with changes in real-time wind direction and speed. Existing static graph networks or statistical models cannot capture this wind-changing topology, making it difficult for the models to identify transient causal relationships between upstream and downstream sensors when facing variable wind fields. Furthermore, purely data-driven models lack physical perception of geographical information and ignore the flow effect of fluids around buildings, which often leads to the "wall-penetration fallacy" in the inversion process. That is, the predicted pollution transmission path violates the common sense of physical isolation, thus seriously affecting the accuracy and credibility of the source tracing results.

[0004] Therefore, there is an urgent need for a new spatiotemporal source tracing method that can integrate geographic information physical constraints, adapt to unsteady micro-meteorological changes, and process sparse unstructured data. Summary of the Invention

[0005] This application provides a machine learning-based method and system for tracing pollution sources in chemical industrial parks, which at least solves the above-mentioned technical problems existing in the prior art.

[0006] According to a first aspect of this application, a machine learning-based method for tracing pollution sources in chemical industrial parks is provided, comprising: S1: acquiring raw sensor data streams and geographic information system (GIS) data from the chemical industrial park; S2: performing spatiotemporal alignment and sliding window tensor quantization on the raw sensor data streams based on the GIS data to obtain a standardized input feature tensor and a static distance matrix; S3: using real-time wind field information in the standardized input feature tensor and building outline information in the GIS data to perform physical mask filtering and Gaussian plume weight estimation on the static distance matrix to obtain a dynamic adjacency matrix sequence; S4: synchronously inputting the standardized input feature tensor and the dynamic adjacency matrix sequence into a graph convolutional neural network model embedded with long short-term memory units for spatiotemporal causal feature extraction to obtain a spatiotemporal latent feature vector; S5: performing pollution source probability heatmap inversion on the spatiotemporal latent feature vector to obtain a source probability distribution heatmap; S6: based on the node connection logic in the dynamic adjacency matrix sequence, performing coordinate correction based on continuity constraints on high-probability areas in the source probability distribution heatmap to obtain the final pollution source coordinates.

[0007] According to a second aspect of this application, a machine learning-based pollution source tracing system for chemical industrial parks is provided, comprising: a chemical industrial park data acquisition module for acquiring raw sensor data streams and geographic information system (GIS) data from the chemical industrial park; a data alignment and tensor quantization module for performing spatiotemporal alignment and sliding window tensor quantization on the raw sensor data streams based on GIS data to obtain a standardized input feature tensor and a static distance matrix; and a dynamic adjacency matrix generation module for using real-time wind field information in the standardized input feature tensor and building outline information in the GIS data to perform physical mask filtering and Gaussian plume generation on the static distance matrix. The system includes a weight estimation module to obtain a dynamic adjacency matrix sequence; a spatiotemporal latent feature extraction module to synchronously input the standardized input feature tensor and the dynamic adjacency matrix sequence into a graph convolutional neural network model embedded with long short-term memory units for spatiotemporal causal feature extraction to obtain spatiotemporal latent feature vectors; a source probability distribution heatmap analysis module to invert the pollution source probability heatmap from the spatiotemporal latent feature vectors to obtain a source probability distribution heatmap; and a pollution source coordinate determination module to perform coordinate correction based on continuity constraints on high-probability regions in the source probability distribution heatmap based on the node connection logic in the dynamic adjacency matrix sequence to obtain the final pollution source coordinates.

[0008] Compared with existing technologies, this application constructs an expert question-and-answer assisted interactive system centered on a large-scale model in the water resources field. Addressing the problem of traditional models failing due to the sparse, unstructured distribution of sensors and unsteady micro-meteorological interference in chemical industrial parks, this solution abandons traditional static Euclidean distance modeling. Instead, it utilizes building outline information and real-time wind field data from GIS to perform physical mask filtering and Gaussian plume weight reconstruction on the static distance matrix, generating a dynamic adjacency matrix sequence that characterizes airflow around-the-wall effects and dynamic causal relationships. This effectively avoids the wall-penetrating fallacy and captures the transient topological structure with the wind. Based on this, a graph convolutional neural network with embedded long short-term memory units is constructed. The standardized sensor spatiotemporal tensor and the dynamic adjacency matrix are input synchronously to extract spatiotemporal causal features in a non-Euclidean space. Finally, by inverting the source probability distribution heatmap and combining it with continuity constraints for coordinate correction, accurate location of pollution sources is achieved under complex interference environments. Attached Figure Description

[0009] The above and other objects, features, and advantages of exemplary embodiments of this application will become readily understood by reading the following detailed description with reference to the accompanying drawings. Several embodiments of this application are illustrated in the drawings by way of example and not limitation, wherein: in the drawings, the same or corresponding reference numerals denote the same or corresponding parts.

[0010] Figure 1 This is a flowchart of a machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application; Figure 2 This is a data flow diagram illustrating the machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application; Figure 3 This is a flowchart of step S3 in the machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application; Figure 4 This is a block diagram of a machine learning-based pollution tracing system for chemical industrial parks according to an embodiment of this application. Detailed Implementation

[0011] To further illustrate the technical means and effects adopted by this application in order to achieve the intended purpose, the following detailed description of the specific implementation methods, structures, features and effects of this application is provided in conjunction with the accompanying drawings and preferred embodiments.

[0012] To address the shortcomings in the aforementioned technical fields, this application provides a machine learning-based method for tracing pollution sources in chemical industrial parks, such as... Figure 1 and Figure 2 As shown, Figure 1 This is a flowchart of a machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application. Figure 2This is a schematic diagram of data flow for a machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application.

[0013] Specifically, step S1 involves acquiring raw sensor data streams and geographic information system (GIS) data from the chemical industrial park. Specifically, the raw sensor data streams include pollutant concentrations, wind speeds, and wind directions at each monitoring station; the GIS data includes metadata for the monitoring stations and geographic information about the park. It should be understood that in the extremely complex physical environment of a chemical industrial park, the atmospheric diffusion of pollutants is not driven solely by a single concentration gradient, but is constrained by both micro-meteorological conditions and the topology of surface buildings. The dense factory facilities within the park create highly irregular airflow channels, rendering traditional Euclidean distance metrics ineffective, and sparsely distributed sensors struggle to independently characterize continuous pollution fields without physical spatial context. To enable the subsequently constructed machine learning model to perceive physical boundaries and dynamic airflow causality, this application first constructs a multi-source heterogeneous data foundation that integrates real-time environmental conditions and static spatial geometric features, thereby providing accurate node attributes and edge connection criteria for the spatiotemporal graph neural network. Based on this, step S1 is implemented.

[0014] In one specific implementation, step S1 of this application is performed as follows: First, the acquisition of the raw sensor data stream is achieved by accessing the IoT sensing layer interface deployed at various key locations in the chemical industrial park. This sensing layer covers fixed gas monitoring stations and micro weather stations, and the data stream is aggregated to the data center in real time via an industrial wireless network (such as LoRaWAN or 4G / 5G private network). Logically, the raw sensor data stream consists of a series of data packets with timestamps. Each monitoring station's data packet contains three core physical quantities: the real-time concentration value of the target pollutant, the on-site wind speed, and the on-site wind direction. For example, in a scenario of tracing the source of a certain volatile organic compound (VOCs), the preset sampling frequency may be set to transmit one frame of data every 30 seconds. Taking the sensor numbered Node_01 as an example, the data frame uploaded at a certain moment contains: a VOCs concentration value of 0.45 mg / m³. 3 The wind speed was 2.1 m / s, and the wind direction was 270 degrees, i.e., westerly (airflow from west to east). These continuously uploaded data frames constituted an unfiltered and aligned time series stream, which could reflect the air quality and micrometeorological conditions at various points in the park over time.

[0015] Meanwhile, the acquisition of GIS data focuses on constructing the static physical framework of the park, which forms the geometric basis for subsequent calculations of physical masking and wall penetration. This data is mainly obtained by parsing high-precision vector map files of the park (such as Shapefile or GeoJSON format), and is specifically divided into two parts: metadata of monitoring stations and geographic information of the park. The metadata of monitoring stations mainly refers to the mapping table between sensor IDs and geographic coordinates. By querying the equipment ledger, the precise planar coordinates of each sensor node on the map (e.g., X and Y coordinates under the universal transverse Mercator projection) are extracted. The geographic information of the park contains vector outline data of all major buildings in the park. This data is stored in the form of polygons, each polygon being defined by a series of closed vertex coordinate sequences, accurately describing the footprint and geometric boundaries of factories, warehouses, and large installations. For example, a factory building is described as a rectangular polygon in GIS data, and its vertex coordinates establish that the area is a physical entity that airflow cannot penetrate. During the acquisition process, it is necessary to ensure that the sensor coordinate system is consistent with the building outline coordinate system. If there are differences between latitude and longitude coordinates and projected coordinates, a standard geographic coordinate transformation algorithm must be applied at this stage to unify them, thereby outputting a set of geographic information data with a unified benchmark.

[0016] Specifically, step S2 involves performing spatiotemporal alignment and sliding window tensor quantization on the raw sensor data stream based on Geographic Information System (GIS) data to obtain a standardized input feature tensor and a static distance matrix. Correspondingly, the original form of multi-source heterogeneous data within chemical industrial parks often exhibits significant spatiotemporal discreteness and dimensional differences. This is mainly reflected in the non-uniform topological distribution of sensor nodes in physical space, and the asynchronous transmission and sampling frequency jitter of monitoring data in the time dimension. Since the subsequent graph convolutional neural network and long short-term memory network models require input data to have a strict tensor structure and aligned time steps, and the models are sensitive to numerical scale, directly using the raw unstructured data can lead to difficulties in gradient descent convergence and failure in spatial feature extraction. To transform discrete physical monitoring points into graph structure nodes understandable by deep learning models, and to transform asynchronous time series into synchronous spatiotemporal feature maps, it is necessary to construct a unified spatiotemporal benchmark and feature engineering pipeline, thereby providing a standardized, high-quality data foundation for subsequently capturing dynamic causal relationships under unsteady micro-meteorological conditions. Therefore, step S2 is implemented.

[0017] In one specific implementation, step S2 includes: S21, parsing the planar coordinates of each monitoring station from the geographic information system data; S22, constructing a static distance matrix based on the planar coordinates of each monitoring station; S23, performing unified time axis resampling and missing value imputation on the original sensor data stream to obtain an aligned sensor data stream; S24, performing orthogonal vector decomposition on the wind direction angle in the aligned sensor data stream and performing maximum and minimum normalization processing on the pollutant concentration, and performing sliding slicing and stacking encapsulation on the processed multidimensional data based on a preset time window length to obtain a standardized input feature tensor.

[0018] Step S2 of this application is executed as follows: First, sub-step S21 is executed. In this process, the input data comes from the Geographic Information System (GIS) data obtained in step S1, which is stored as a vector file in ESRI Shapefile or GeoJSON format. Since the original GIS data uses a latitude and longitude coordinate system (such as WGS84), and its unit is degrees, it cannot be directly used for accurate Euclidean distance calculations. Therefore, a map projection transformation is first required. The layer containing the monitoring station metadata is read, and the spherical coordinates are converted to metric coordinates in a Cartesian coordinate system using the Gauss-Kruger projection or the general transverse Mercator projection algorithm. For example, if there are N fixed sensor nodes within a chemical industrial park, ... This is indicated for each sensor node. Extract its unique device ID from the GIS attribute table and parse the planar coordinates of its geometric center point. For example, for node Node_01 mentioned in step S1, the planar coordinates obtained after projection transformation are (45020.5, 32100.8), in meters. This process traverses all registered monitoring stations within the park, generating a coordinate list containing N nodes. This list not only defines the absolute position of the sensor in the physical world but also provides the geometric basis for the subsequent construction of the adjacency relationship of the graph neural network.

[0019] Next, sub-step S22 is executed. Based on the planar coordinates of each node output in sub-step S21... Construct an N×N square matrix, denoted as the static distance matrix. Each element in the matrix This represents the straight-line Euclidean distance between the i-th sensor node and the j-th sensor node. It is calculated using the Pythagorean theorem. This reflects the pure geometric distance without considering any obstacles. Because distance is symmetric—that is, the distance from node i to node j is equal to the distance from node j to node i—this matrix... Let be a real symmetric matrix, and the elements on the main diagonal are... All values ​​are 0, meaning the distance from a node to itself is 0. For example, if the coordinates of node_01 are (45020.5, 32100.8) and the coordinates of node_02 are (45220.5, 32100.8), then the difference in their x-coordinates is 200 meters, and the difference in their y-coordinates is 0 meters. The calculated distance elements... =200 meters. This calculation process is performed pairwise across all N node pairs. If N=50, a 50×50 numerical matrix is ​​generated.

[0020] Subsequently, sub-step S23 is executed. The raw sensor data stream acquired in step S1 is often asynchronous, with upload times from different devices potentially varying by several seconds to minutes. Furthermore, due to network jitter, data packets may be lost or out of order. To meet the requirements of Long Short-Term Memory (LSTM) networks for equally spaced time series inputs, this step establishes a unified time reference axis. First, a uniform time step is set across the entire campus. Based on the speed of microclimate changes and the characteristic scale of pollution diffusion in the chemical industrial park, a set... =5 minutes. This means the data will be reconstructed into a series of standard timestamps. ,in For each sensor node, its original irregular time series is mapped onto this unified time axis. Since the original sampling times often do not precisely fall on the standard timestamp (e.g., the standard time is 10:00:00, but actual sampling times are 10:00:12 and 09:59:48), linear interpolation is used to estimate the standard time value. By calculating the time proportion of the target time between two consecutive sampling times, the values ​​are synthesized proportionally to ensure data continuity and smoothness. After resampling, it may be found that some time periods are without data for extended periods due to equipment maintenance or malfunctions, forming gaps. Missing values ​​need to be filled in these cases. For short-term gaps, such as fewer than three consecutive time steps, a forward nearest neighbor filling strategy is used, assuming short-term inertia in the environmental state and filling subsequent gaps with the most recent valid observation. For longer-term gaps, they are marked as invalid periods or filled using the spatial Kriging interpolation mean of all other stations in the entire park. After the above processing, the original discrete data stream is transformed into an aligned sensor data stream. This data stream is a three-dimensional array structure with dimensions of [dimensions missing in original text]. ,in This represents the total number of time steps. For the number of nodes, The original number of feature channels is 3: pollutant concentration, wind speed, and wind direction / angle. Taking Node_01 as an example, at 10:00:00, its data has been confirmed as: concentration 0.45 mg / m³. 3 Wind speed 2.1 m / s, wind direction 270 degrees.

[0021] Finally, execute sub-step S24. For the wind direction data, the original wind direction... The wind speed is represented by angle values ​​ranging from 0 to 360 degrees. This representation presents a numerical chasm; 359 degrees and 1 degree are physically very close (only 2 degrees apart), but numerically they differ significantly (358 degrees). Directly inputting this into a neural network would cause a large gradient error when processing the northerly wind component, disrupting the continuity of the physical meaning. Therefore, an orthogonal vector decomposition method is used to decompose the wind speed scalar... Angle of wind direction The joint mapping is a wind speed component vector on a two-dimensional plane. The decomposition formula is as follows: First, convert the meteorological wind direction angle into a standard mathematical polar coordinate angle (e.g., a meteorological 270-degree westerly wind corresponds to 0 degrees or the positive X-axis in the mathematical coordinate system), or use a decomposition formula that fits the meteorological definition to ensure that the direction of the decomposed vector is consistent with the actual physical flow direction. This represents the east-west component of wind speed (positive values ​​indicate eastward-flowing winds, and negative values ​​indicate westward-flowing winds). This represents the north-south component of wind speed (positive values ​​indicate northward-flowing wind, negative values ​​indicate southward-flowing wind). Through this transformation, 359 degrees and 1 degree are mapped to similar vector coordinates, thus ensuring the smooth continuity of the input features in the numerical space. For example, with a wind speed of 2.1 m / s and a wind direction of 270 degrees, the calculated... =2.1m / s, flowing eastward. =0 m / s. (Regarding pollutant concentration data) Due to its numerical range (e.g., 0-100 mg / m³) 3 Since the components of the model (e.g., -10 to 10 m / s) have different dimensions from the wind speed components, to accelerate model convergence and prevent large numerical features from dominating the loss function, max-min normalization is used. The normalization formula is as follows: ,in and These represent the minimum and maximum observed values ​​of this pollutant in the historical dataset. After processing, the concentration values ​​are scaled to the [0,1] interval. Wind speed component. and Similar normalization or standardization processes can be applied to distribute them within a similar numerical range. At this point, each node at each time step t not only possesses location information but also a processed feature vector. The feature dimension F=3. To capture the spatiotemporal causal effects of pollutant transport (i.e., upstream emissions take time to reach downstream areas), the model cannot rely solely on current data but must input a historical observation window. The length of the historical observation window is set. For example, 12 (corresponding to data from the past 60 minutes, with a step size of 5 minutes). Using a sliding window technique, the window slides along the time axis T, and for each target time t, a segment is extracted from... arrive Continuity of time A sequence of feature vectors. All N nodes in The data within the window is stacked and encapsulated to construct a standardized input feature tensor. The dimension of this tensor is or This is used for single-run inference. Where B is the batch size and N is the number of nodes, which is 50 in this application. The time series length is 12 in this application, and F is the number of feature channels, which is 3 in this application. This four-dimensional tensor completely preserves the concentration evolution trend and wind field vector change trajectory of all monitoring points in the park over the past hour, and the data has been spatiotemporally aligned and numerically standardized.

[0022] Specifically, step S3: Using real-time wind field information from the standardized input feature tensor and building outline information from the geographic information system data, physical mask filtering and Gaussian plume weight estimation are performed on the static distance matrix to obtain a dynamic adjacency matrix sequence. It is understandable that the migration and diffusion behavior of pollutants is constrained by both micro-meteorological conditions and macro-physical structures. Relying solely on the static distance matrix in Euclidean geometry cannot accurately reflect the transport patterns of matter in fluid media, because pollutants do not diffuse uniformly in a straight line, but strictly follow an advection transport mechanism following the wind direction, and are affected by physical barriers and flow around physical structures. If the machine learning model only aggregates features based on fixed physical distances between monitoring stations, it will inevitably introduce backwind or through-wall noise that violates physical common sense, severely weakening the accuracy and physical interpretability of the source tracing results. To endow deep neural networks with the ability to perceive dynamic flow fields and physical boundaries, a graph topology that can evolve over time and contains physical constraints needs to be constructed. Therefore, step S3 is implemented.

[0023] Figure 3 This is a flowchart of step S3 in the machine learning-based pollution source tracing method for chemical industrial parks according to an embodiment of this application. Figure 3As shown, in one specific embodiment, step S3 includes: S31, extracting the global dominant wind field vector sequence from the standardized input feature tensor to obtain the global wind field vector sequence at each time step; S32, using the building vector contours in the geographic information system data to perform ray projection detection on the connections between sensor nodes, and calculating the cosine value of the wind direction deviation between sensor nodes in combination with the global wind field vector sequence to obtain the physical blocking mask matrix and the wind direction deviation factor sequence; S33, processing the static distance matrix based on the wind direction deviation factor sequence and the physical blocking mask matrix to obtain the dynamic adjacency matrix sequence.

[0024] Step S3 of this application is executed as follows: First, sub-step S31 is executed. Although each sensor node in the park records local microscopic wind field data, the wind speed readings of a single sensor are easily affected by turbulent vortices caused by surrounding local obstacles (such as trees and small equipment), making them unable to represent the macroscopic transport flow field at the park scale. In order to extract the dominant airflow direction driving the large-scale spatial migration of pollutants, spatial integration of the wind field data is required. In specific implementation, firstly, the standardized input feature tensor is used... Feature Dimensions The middle slice extracts the channel data with indices 1 and 2, that is, the data of all nodes at all time steps. and Numerical values. For each specific time step t (t iterates from 1 to...), Select the set of wind speed vectors from all N sensors within the park at that moment. The average wind speed vector for the entire park at that moment is calculated using the vector averaging algorithm. By calculating the arithmetic mean of the wind speed components at all spatially discrete nodes, random Gaussian noise and micro-turbulence components caused by local disturbances can be effectively smoothed out, thus preserving the dominant wind direction characteristics that reflect the overall motion trend of the atmospheric boundary layer. This process is performed independently for each time slice within the time window, ultimately... The global wind speed vectors calculated at each time step are arranged in temporal order to generate a global wind field vector sequence. . =12 corresponds to one hour of data, so the sequence contains 12 two-dimensional vectors, which fully describe the evolution trajectory of the wind field over time in the past hour.

[0025] Next, sub-step S32 is executed. In the physical obstruction detection stage, the set of vector polygon outlines of all buildings or enclosed structures marked as large buildings or enclosed devices is first extracted from the geographic information system data obtained in step S1. , Let N be the total number of polygons in the building. Then, iterate through the N sensor nodes determined in step S2, constructing a virtual connection segment for any two node pairs (i,j). Line segments are detected using ray casting algorithms or line segment intersection algorithms from computational geometry. Is it related to the set? any polygon in An intersection occurs. The determination logic is as follows: if the line segments... If a point intersects with the edge of any building polygon, it means that the fluid cannot be directly transmitted from node i to node j along a straight line, i.e., there is a physical barrier. In this case, a blocking mask element is set. =0; Conversely, if the line segment does not intersect any building, it indicates that there is physical visibility or airflow directness between the two points, setting... =1. After traversing and calculating all N×N node pairs, a static binary matrix—the physical blocking mask matrix—is constructed. This matrix is ​​constant over time unless the building layout changes. In the wind direction deviation calculation, the aim is to quantify the degree to which the relative positions between nodes match the current wind direction. For any two sensor nodes i (source hypothetical point) and j (target hypothetical point), their planar coordinates are first used... and Calculate the geometric azimuth angle from i to j. This angle is defined as a vector. The angle relative to true north (or the positive X-axis, consistent with the definition of wind direction) can be calculated using the inverse trigonometric function arctan2. Next, the global wind field vector sequence generated in sub-step S31 is introduced. For each time step t, the current global prevailing wind direction angle is extracted. This angle is determined by the global wind speed vector. The calculation shows that, To determine whether node j is downwind of node i, the geometric azimuth angle is calculated. with wind angle The cosine of the difference between the two values ​​is the wind direction deviation factor. The calculation formula is: Analyzing the formula: This formula calculates vectors. (Node i points to j) and vector Cosine similarity between (wind direction vectors). If Approaching 1 indicates That is, the wind is blowing from i to j, and node j is located strictly downwind of node i, indicating a very strong causal relationship; if A value close to -1 indicates that the two directions are opposite, with node j located upwind of node i. According to the principle of advection diffusion, the upwind node cannot be polluted by emissions from the downwind node (ignoring turbulent diffusion at extremely low wind speeds), indicating a very weak causal relationship. A value of 0 indicates crosswind. This calculation is performed for each time step t and for each pair of nodes (i,j), generating a three-dimensional sequence of wind direction deviation factors with dimensions of [missing information]. .

[0026] Finally, sub-step S33 is executed. This step is the synthesis stage for generating the physically-aware graph structure. Its core idea is to construct a dynamic weight function that comprehensively considers distance attenuation, wind alignment, and physical blockage, based on the physical equations of the Gaussian plume model. In traditional graph convolutional networks, the adjacency matrix is ​​determined by distance, but in this scheme, the connection weight between any two nodes i and j at time t is determined by distance. It is determined by the product of three factors: distance attenuation, wind directionality, and physical constraints. The specific calculation formula is as follows: A detailed analysis of each term in the formula: The first term is the distance attenuation term, based on the elements in the static distance matrix obtained in step S2. Calculation. Using the Gaussian kernel function. ,in The preset bandwidth parameter, such as setting it to 50% of the average sensor spacing in the park, controls the attenuation rate of spatial correlation. The closer the distance ( The smaller the value, the closer this term is to 1; the farther the distance, the more exponentially the weight decreases. This reflects the dilution effect during the diffusion process. The second term is the wind directionality term, which uses the ReLU activation function to process the wind direction deviation factor obtained in step S32. .Right now This operation has crucial physical significance: when node j is downwind (positive cosine value), its value is retained as a weight coefficient, with a higher weight corresponding to the wind direction; when node j is upwind or in the opposite direction of the vertical wind (negative cosine value), the ReLU function forces it to 0. This completely severs the upwind connection in the graph structure, transforming the undirected graph into a directed graph, forcing information flow to propagate only along the wind direction, thus conforming to the physical causality of matter transmission. The third term is the physical constraint term, directly multiplied by the physical blocking mask generated in step S32. If there is a building blocking the way between two points, that is... =0, then regardless of how close the distance or how consistent the wind direction, the final connection weight is... All will become 0. This simulates the effect of buildings on the cutoff of airflow and pollutant plumes. The synthesis process iterates along the time dimension, calculating the corresponding weight matrix for each time step t corresponding to the input tensor using the global prevailing wind direction at that time, and finally stacking them to generate a four-dimensional or three-dimensional tensor object, i.e., a dynamic adjacency matrix sequence. Its dimensions are To illustrate with a specific example: In this application =12, N=50. At time t=5, the prevailing wind direction in the entire park is westerly (270 degrees), and the airflow direction is due east. For sensor node 1 and node 2, the distance between them is 200 meters ( =200), and node 2 is located directly east of node 1 (i.e., downwind), with no buildings obstructing the view. If the distance term is calculated to be 0.6, the wind direction term is cos(0)=1, and the obstruction term is 1, then the final weight... =0.6×1×1=0.6. For the connection from node 2 to node 1, since it's a headwind, the wind direction term is... =0, therefore =0. For nodes 1 and 3, although they are only 100 meters apart and downwind, there is a tall factory building in between, so the obstruction term is 0. =0. Thus, the original static distance matrix is ​​transformed into a series of sparse adjacency matrices that dynamically evolve over time. This sequence of dynamic adjacency matrices... The topology of possible pollutant transport paths at each moment is precisely encoded.

[0027] Specifically, step S4 involves simultaneously inputting the standardized input feature tensor and the dynamic adjacency matrix sequence into a graph convolutional neural network model embedded with long short-term memory units for spatiotemporal causal feature extraction to obtain spatiotemporal latent feature vectors. It should be understood that in the extremely complex environmental flow field of a chemical industrial park, the diffusion process of pollutants exhibits highly spatiotemporal coupling characteristics. Sensor readings at a single point are not isolated events, but rather the physical superposition of historical emissions from the upwind region after advection transport and local turbulent diffusion at the current moment. Although steps S2 and S3 have constructed standardized feature tensors and dynamic adjacency matrices containing wind field physical constraints, these data are still in a shallow discrete representation state and have not yet numerically simulated the dynamic mixing process of pollutants in a non-Euclidean topological space. To automatically learn complex nonlinear diffusion laws from the data and elevate the discrete monitoring point information to a higher-order semantic space capable of representing the global pollution situation, a deep learning architecture capable of simultaneously handling graph topology and time-series dependencies must be adopted. Therefore, step S4 is implemented.

[0028] In one specific implementation, step S4 includes: S41, performing time-step graph convolution aggregation and Laplacian normalization on the standardized input feature tensor using the dynamic adjacency matrix sequence to obtain a spatial feature sequence tensor; S42, flattening the spatial feature sequence tensor and inputting it into a long short-term memory network unit for temporal state evolution and causal accumulation, and extracting the hidden state vector of the last time step to obtain a temporal hidden state sequence; S43, using a fully connected layer to perform dimensionality reduction compression and nonlinear projection on the temporal hidden state sequence to obtain a spatiotemporal hidden feature vector.

[0029] Step S4 of this application is executed as follows: First, sub-step S41 is executed, the main purpose of which is to simulate the spatial flow effect of pollutants within a single time slice by utilizing the spectral or spatial aggregation mechanism of graph convolutional networks. Since the wind field in the chemical industrial park changes dynamically over time, the connection relationships between nodes... Since the time step t is different at each time step, static graph convolution cannot be used. Instead, the tensor is unfolded along the time axis, and dynamic graph convolution is performed step-by-step. First, the four-dimensional feature tensor is... and three-dimensional adjacency matrix sequence Along the time dimension Perform slice decomposition. For each discrete time step t, where t iterates from 1 to 12, extract the node feature matrix at the current time step. and instantaneous adjacency matrix The characteristic matrix at this time It includes normalized concentration, wind speed X-component, and Y-component data from 50 sensors within the park at that moment; adjacency matrix. The magnitude of its element values ​​precisely describes the probability weight of material transport between any two nodes under the current wind direction and building obstruction conditions. Next, to ensure the numerical stability of features during propagation in the multi-layer network and prevent feature amplitude explosion or disappearance due to differences in node degree (i.e., different connection numbers), the instantaneous adjacency matrix needs to be adjusted. Perform Laplace normalization. Before the actual calculation, first... Add the identity matrix The enhanced adjacency matrix is ​​obtained. The physical significance of this step lies in adding self-loops to the graph structure, ensuring that during feature aggregation, nodes not only receive information from their neighbors but also retain their original feature information. Subsequently, the calculation... degree matrix The degree matrix is ​​a diagonal matrix, and its elements on the main diagonal are... It is equal to the sum of the elements in the i-th row of the enhanced adjacency matrix, that is... This value reflects the total connectivity strength of node i in the current wind field network. Based on this, a normalized Laplace matrix is ​​constructed. The calculation formula is as follows: This operation is a symmetric normalization. For directed graph structures, random walk normalization can also be used. To conform to directed transport logic, the weights in the adjacency matrix are scaled by the node degree by left and right multiplying by the inverse square root of the degree matrix, respectively. If a node is connected to a large number of upstream nodes, normalization appropriately reduces the weight of each connection, thus maintaining the conservation of overall input energy. After normalizing the graph structure, graph convolution aggregation is performed. This is a crucial step in feature extraction, designed to flow feature information from upstream nodes to downstream nodes. The model initializes a learnable weight matrix. ,in The input feature dimension is 3. The dimension of the output spatial latent features is set to 64. The parameters of this weight matrix are initialized using a Glorot uniform distribution and continuously adjusted during model training using the backpropagation algorithm to learn the nonlinear coefficients in the pollutant diffusion equation. The formula for calculating graph convolution is: In this formula, matrix multiplication A spatial aggregation process was physically performed: for node i, its new eigenvector is obtained by weighted summation of the eigenvectors of all its connected nodes j, i.e., nodes upwind and unobstructed. This step simulates the physical process of matter being transported by wind from its neighborhood to the current point. Subsequently, the result is compared with the weight matrix. Multiplication completes the linear transformation and dimensional mapping of the feature space. Finally, the ReLU activation function is applied. By introducing nonlinear capabilities and eliminating negative responses, we obtain the spatial characteristic matrix after spatial mixing at the current time t. Taking a specific numerical example, at time t=5, node 2 is located downwind and east of node 1, and It has a relatively high weight. After... After the operation, the high-concentration features of node 1 are weighted and transferred to node 2. If node 2 originally had a lower concentration, but is located downwind of the high-concentration node 1, after convolution, node 2's representation in the latent feature space will contain potential state information about the impending contamination. Finally, all... The spatial feature matrix calculated at each time step Reassembled according to the original time sequence, it is reconstructed into a new three-dimensional tensor, namely the spatial feature sequence tensor. Its dimensions become In this application, the tensor is (1,50,12,64). Each slice in this tensor represents the spatial distribution characteristics at a specific moment after considering the wind field morphology and physical barriers of the entire park, thus eliminating the spatial discreteness of the original data.

[0030] Next, sub-step S42 is executed. First, the spatial feature sequence tensor from S41 is received. Continuing with the previous numerical example, 50 sensor nodes (N=50) were deployed within the chemical industrial park, with a historical observation window length of 12 time steps. =12), and the spatial feature dimension after graph convolution aggregation is 64 ( =64). At this point, each time slice in the tensor The 50×64 dimension contains the local feature distribution at that moment after physical constraints and spatial mixing. To enable the subsequent Long Short-Term Memory (LSTM) network to process this complex topological data, a feature flattening operation is first performed. At each discrete time step t, the 50×64 feature matrix is ​​flattened into a one-dimensional high-dimensional vector of length 3200. This operation is not merely a simple data rearrangement; its core purpose is to construct a global pollution fingerprint that represents the pollution distribution of the entire chemical industrial park at time t. In this long vector, the first 64 elements represent the features of node 1, the next 64 elements represent the features of node 2, and so on, up to node 50. Since the order of the nodes remains strictly consistent across all time steps, this flattening operation implicitly preserves the independent identity information and topological index of the sensor nodes, which is crucial for subsequent model inversion to locate the specific source coordinates. After this processing, the original four-dimensional tensor is reconstructed into a three-dimensional temporal input tensor with dimensions (1, 12, 3200), which is a sequence of global feature vectors containing 12 consecutive time steps. The sequence is then fed into a Long Short-Term Memory (LSTM) network unit for temporal state evolution and causal accumulation. LSTM is a variant of recurrent neural networks specifically designed to solve the vanishing gradient problem in long sequence training. Its core architecture includes a forget gate, an input gate, an output gate, and a cell state that spans all time steps. The model initializes a series of learnable parameters for the LSTM unit, including the weight matrix and bias terms, which are obtained using the Glorot uniform distribution or orthogonal initialization method. The input vector at the first time step t=1 is then processed. When the global fingerprint is contaminated at time t=1, the LSTM unit first uses a forget gate to determine which information from the previous time step needs to be discarded (for t=1, the previous time step state is initialized to a zero vector). The forget gate uses the Sigmoid activation function to read the input. and the hidden state of the previous moment The output is a numerical vector between 0 and 1, controlling the proportion of cell state information retained. Next, the input gate determines how much of the new information at the current time step needs to be stored in the cell state. This involves two operations: first, using a Sigmoid layer to determine which values ​​to update; and second, using a Tanh layer to create a new candidate cell state vector. (Old cell state) Under the control of the forget gate, invalid historical information is removed, and new features filtered by the input gate are added to update the cell state to the current moment. This mechanism allows the model to remember long-term lag effects. For example, the model can remember the high concentration value that appeared in the northwest corner of the cell at time t=1 and retain it until time t=5, so as to establish a causal relationship with the concentration peak that appeared in the southeast corner at that time. Finally, the output gate is based on the current cell state. and input The hidden state vector at the current time step is calculated by processing the data through the Tanh function and gating it with the Sigmoid function. This process is executed cyclically along the time axis from t=1 to t=12, with the hidden state at each time step... All of these incorporate all historical accumulated information and spatial evolution characteristics from the initial time to the current time. Ultimately, the LSTM unit outputs a complete temporal hidden state sequence. Each of them The dimensions are set to be consistent with the input dimensions or have a sufficiently large capacity, such as 3200 dimensions, to ensure lossless transmission of information.

[0031] Finally, sub-step S43 is executed. First, the hidden state vector of the last time step is extracted from the temporal hidden state sequence generated in step S42. In the temporal modeling theory of deep learning, the last output state of a sequence model is considered a compressed summary of the entire input sequence, which theoretically encodes the entire time window. This refers to the complete spatiotemporal context information within the past 60 minutes regarding where the current pollution distribution was caused by emissions and the transmission routes taken. At this point... Although it contains rich causal logic, its dimensionality is still as high as 3200 dimensions, and it contains a large amount of redundant background noise and non-critical features. To extract the core causal features and adapt them to the subsequent source inversion task, fully connected layers are needed. Dimensionality reduction and nonlinear projection are performed. The fully connected layer contains a weight matrix. and a bias vector If the latent space dimension is preset to 128 dimensions (i.e., the target output dimension), then the weight matrix... The shape is 3200×128, and the bias vector is... The length is 128. The calculation process is as follows: First, a linear transformation is performed to transform the high-dimensional vector... With weight matrix Multiply and add bias. This linear operation maps feature information scattered across 3200 dimensions into a compact 128-dimensional subspace, forcing the model to learn the most discriminative principal components in the data and removing irrelevant interference such as minor fluctuations in wind speed or sensor white noise. Subsequently, a nonlinear activation function (such as ReLU or Tanh) is applied to process the result of the linear transformation to increase the model's nonlinear expressive power, enabling it to fit complex pollution source-receptor response relationships. Finally, the calculated 128-dimensional vector is the spatiotemporal latent feature vector. This vector is a highly condensed mathematical representation that no longer contains intuitive physical quantities (such as specific concentration or wind speed values). Instead, it precisely describes, in an abstract coded form, the current spatiotemporal response pattern of the entire park caused by potential pollution sources under specific wind field conditions.

[0032] Specifically, step S5 involves inverting the spatiotemporal latent feature vector to obtain a heatmap of the source probability distribution. Correspondingly, after the deep feature extraction in the preceding steps, the model has successfully compressed and encoded the complex micro-meteorological fluctuations, physical spatial barriers, and spatiotemporal response patterns of the sensor network within the chemical industrial park into an abstract one-dimensional spatiotemporal latent feature vector. However, while this feature vector in the high-dimensional latent space mathematically contains all the causal logic required for source tracing, its representation is unreadable to environmental management personnel responsible for emergency response and lacks intuitive geospatial orientation. To transform this abstract machine encoding into specific location information that can directly assist decision-making, a decoding mapping mechanism from the latent feature space to the physical coordinate space needs to be constructed. This mechanism reprojects the highly condensed causal features back onto the two-dimensional geographic plane and quantifies the probability of each area within the park as a potential pollution source in the form of a probability distribution. Therefore, step S5 is implemented.

[0033] In one specific implementation, step S5 includes: S51, performing high-dimensional spatial projection and feature decoding on the spatiotemporal latent feature vector to obtain the original response grid; S52, based on the Softmax function, probabilizing and structuring the original response grid to obtain the source probability distribution heatmap.

[0034] Step S5 of this application is executed as follows: First, sub-step S51 is executed. At the beginning of implementation, based on the actual physical scale of the chemical industrial park and the accuracy requirements for source tracing, a virtual grid resolution is predefined. For example, if the area of ​​the chemical industrial park is 5 km × 5 km, in order to ensure that the source tracing results can be accurate to specific factory buildings or equipment units, the grid resolution is set to G × G, where G = 50. This means that the entire continuous geographical plane of the park is discretized into 50 × 50 = 2500 independent grid units, each grid representing a potential emission area of ​​100 m × 100 m. This discretization process provides a unified spatial benchmark for subsequent classification and regression tasks. On this basis, a fully connected layer is used as a decoder to process the input 128-dimensional spatiotemporal latent feature vector. Projected onto this 2500-dimensional grid space. This fully connected layer contains a learnable weight matrix. and a bias vector Among them, the weight matrix The dimension is 2500×128, and each row of weights actually encodes the correlation strength between a specific feature combination and a specific grid coordinate; bias vector The length is 2500. The projection calculation formula is as follows: This formula performs a linear transformation operation. The input is a 128-dimensional feature vector. With weight matrix Multiplication decompresses the compressed semantic information and distributes it to 2500 output channels, generating a one-dimensional grid log-probability vector of length 2500. The value of each element in this vector initially reflects the probability that the corresponding grid cell is a pollution source; the larger the value, the higher the probability. These weight parameters... This is obtained automatically during the model training phase by minimizing the cross-entropy loss function between the predicted source and the true source using the backpropagation algorithm. This results in a one-dimensional vector. Subsequently, the processing module performs a dimensionality transformation operation, rearranging the matrix according to the preset grid width and height (G,G), i.e., (50,50), thereby restoring its two-dimensional spatial structure and obtaining a preliminary two-dimensional matrix. However, the matrix obtained by direct mapping often exhibits numerical jumps between adjacent grids, manifesting as discrete noise points. This contradicts the characteristic that real pollution sources have a certain continuity in physical space (i.e., the source is a region rather than an isolated mathematical point). To ensure the continuity of physical distribution and eliminate high-frequency noise, a small two-dimensional Gaussian convolution kernel is applied to smooth the matrix. Let the Gaussian kernel be... The smoothed matrix is The convolution operation is represented as: The original response mesh obtained after smoothing. It is a 50×50 real number matrix with smooth spatial transitions between its elements, which can better characterize the central location of the pollution source and its edge decay trend.

[0035] Then, sub-step S52 is executed. The original response mesh generated in step S51. The element values ​​in the graph belong to the real number domain (which may include negative numbers), representing only relative response intensities (Logits), lacking clear probabilistic statistical meaning, and the sum of the values ​​across the entire plane is not fixed, making it unsuitable as a direct confidence output. To satisfy the axiomatic definition of a probability distribution, i.e., the probability of all possible locations is non-negative and sums to 1, a spatial softmax function is applied. Unlike conventional image classification tasks where softmax is performed only along the channel dimension, this operation is performed globally across all elements of the entire two-dimensional grid plane (G×G). For the element with coordinates (i,j) in the original response grid... The corresponding source probability The calculation formula is as follows: The formula is analyzed as follows: The numerator uses an exponential function to map real-valued responses to non-negative numbers. Simultaneously, the exponential function's expansion property significantly widens the gap between high and low response values, making the prediction results sharper, i.e., more clearly pointing to the high-probability region. The denominator is the sum of the exponential response values ​​of all elements across the entire 50×50 grid plane, used as a normalization constant. This division operation ensures that each element in the output matrix... All values ​​fall within the interval [0,1], and the sum of the probabilities of all 2500 grids in the entire park is strictly equal to 1. The final output is a 50×50 source probability distribution heatmap. In this heatmap, the grid values ​​for most non-source areas will approach 0 (displayed in cool tones), while the grids of the actual pollution source and its adjacent areas will show higher probability values ​​(displayed in warm tones or highlighted). For example, if the probability value at coordinates (23,14) is 0.85, it means that the model considers the area covered by that grid to be the source of the pollution incident with an 85% confidence level.

[0036] Specifically, step S6: Based on the node connection logic in the dynamic adjacency matrix sequence, the coordinates of high-probability areas in the source probability distribution heatmap are corrected based on continuity constraints to obtain the final pollution source coordinates. In other words, although the source probability distribution heatmap generated in the preceding steps provides a global visual representation of potential pollution sources within the chemical industrial park, this heatmap is essentially based on statistical inference results from a deep learning model in the latent feature space. It does not fully integrate the transmission logic in the physical topology and is limited by the preset virtual grid resolution. Directly reading probability extreme points often faces the dual challenges of false positives and insufficient positioning accuracy. In the actual complex micrometeorological environment, the model may be affected by local noise or instantaneous concentration fluctuations. High-probability peaks are predicted in areas that, while statistically significant, cannot physically transmit pollutants (such as dead corners completely isolated by buildings or edge areas downstream of all sensors). These false sources, violating physical causality, need to be eliminated. Furthermore, since gridding discretizes continuous geographic space, the coordinates of a single maximum grid center can only provide coarse positioning at the meter or ten-meter level, which cannot meet the urgent need for sub-pixel accuracy in emergency response. Therefore, a post-processing mechanism is introduced to perform physical consistency verification of the probability distribution using the topological connection logic of the sensor network under dynamic wind fields, and to transform the discrete grid probabilities into precise geographic coordinates through continuity constraints. To this end, step S6 is implemented.

[0037] In one specific implementation, step S6 includes: S61, selecting several peak points from the source probability distribution heat map based on an adaptive threshold, and calculating the sum of the out-degree intensities of the sensor nodes corresponding to each peak point according to the dynamic adjacency matrix sequence to obtain a corrected candidate list; S62, performing pollution source inversion on the grid coordinates of each candidate point based on the corrected candidate list to obtain the final pollution source coordinates.

[0038] Step S6 of this application is executed as follows: First, sub-step S61 is executed. The core of this step is to use the transmission logic of the physical world to clean the statistical noise of the data world. First, a 50×50 source probability distribution heatmap is read, and all 2500 elements in the matrix are traversed to find the maximum probability value. For example, in a certain source tracing calculation, The threshold is 0.85. To eliminate widespread low-probability noise in the background, caused by sensor baseline drift or environmental background, an adaptive dynamic threshold is set. The value is taken as 50% of the maximum probability value, that is A binary mask matrix with the same dimensions as the heatmap is generated using this threshold. All grid cells with probability values ​​below 0.425 are set to 0, while those above retain their original values, thus identifying high-confidence candidate regions. Within these high-probability regions after masking, instead of directly selecting a single maximum value, the K highest-probability local peaks are extracted to form a candidate point set. Setting K=5 means the model will simultaneously consider the five most likely potential leakage locations. The candidate point set is organized as a list of five tuples, each recording the grid coordinates of that point. For example, (23,14) and its corresponding original probability value. A key issue at this point is that these grid points do not directly exist in the graph neural network topology constructed in step S3; the nodes in the graph structure are actual existing sensors. To utilize the physical connectivity information in the dynamic adjacency matrix sequence for verification, a mapping from the virtual grid to the physical nodes needs to be established. The calculation of each candidate grid point... Find the Euclidean distance between the candidate point and the planar coordinates of all N sensor nodes within the park, and map the candidate point to the index of its spatially nearest sensor node. For example, candidate point 1 is closest to sensor Node_05, so its proxy node index is 5. The physical assumption of this mapping is that if the grid point is a real pollution source, then the nearest sensor should capture a high concentration signal, and this sensor should have the ability to transmit pollutants to other nodes in the wind farm network. Next, physical consistency verification is performed. For each candidate point k, its proxy sensor node is retrieved. The dynamic adjacency matrix sequence generated in step S3 The connection state in the dynamic adjacency matrix. The network is divided into slices at time steps, where each slice A(t) describes the influence weights between nodes under wind direction and obstruction conditions at time t. If a candidate point is a real source, its proxy node should be located upwind within the past time window and have strong outward transmission capabilities, meaning the sum of edge weights (out-degree) pointing to downstream nodes in the dynamic graph should be large. Conversely, if the proxy node is an isolated point (surrounded by buildings) or only acts as a receiver (always downwind, with a large in-degree but a zero out-degree), it indicates that pollutants cannot effectively diffuse from that location to other parts of the network, and the candidate point is very likely a misjudgment of the downwind accumulation effect by the model. Combining the westerly scenario of S3 mentioned above, if a candidate point is mapped to node 1 (located upwind), then node 1 should have a high-strength connection (large out-degree) pointing to downstream node 2 in the dynamic graph, thus verifying the physical rationality of the candidate point; conversely, if the candidate point is mapped to node 2 (located at the very downstream), its out-degree is 0, and it will be suppressed by the physical filter. To quantify this physical transmission capability, the sum of degree strengths is calculated as a physical correction coefficient. Specifically, this is applied to candidate point k and its proxy nodes. The algorithm iterates through the entire time window (e.g., 12 time steps) and all target nodes j (from 1 to N), accumulating the element values ​​in the dynamic adjacency matrix. Then, it uses the hyperbolic tangent function to map the accumulated result to the (0,1) interval, serving as the confidence decay factor for that candidate point. (Probability correction) The calculation formula is as follows: The double summation term in the formula represents the total output energy of the proxy node throughout the entire historical observation window. The out-degree at time t reflects the instantaneous range of the node's influence on other sensors in the entire park; summing over time t reflects the duration of this influence. Used to calculate time averages. It is a preset scaling hyperparameter, such as 1.0, used to adjust the sensitivity of the correction. The function plays a crucial role: when the sum of the out-degrees is very large (i.e., the point is upwind and unobstructed, making it very easy to spread), the diffusion is very easy. The value approaches 1, the corrected probability Keep close to the original probability When the sum of out-degrees approaches 0 (i.e., the point is blocked or in a dead zone). The value approaches 0, and the corrected probability It will be forcibly decayed to 0. This operation is equivalent to applying a physical filter to the purely data-driven probability result, eliminating artifacts that could not possibly be the source in fluid dynamics. After calculating all K candidate points, a corrected candidate list is generated, which contains a set of candidate points that have been physically and logically cleaned and reweighted.

[0039] Subsequently, sub-step S62 is executed. After obtaining the physically verified candidate list, if only the center coordinates of the grid point with the highest correction probability are simply selected as the result, the tracing accuracy will be limited by the grid resolution, such as 100 meters. To overcome this discretization limitation and achieve sub-pixel-level precise positioning, the processing module no longer relies on a single maximum point, but instead uses a weighted centroid method, comprehensively utilizing the information of all non-zero probability candidate points in the list for spatial interpolation. In specific implementation, all corrected probabilities are selected from the corrected candidate list. Valid candidate points. For example, after filtering, M points remain ( Using the corrected probabilities of these points as weights, a weighted average is calculated for their grid center coordinates. The final pollution source coordinates are then obtained. The calculation formula is as follows: Through this weighted averaging, the final predicted coordinates are no longer confined to the geometric center of the grid, but rather shift towards directions with higher probability density based on the probability distribution, thus falling within the continuous space between grid cells. For example, if candidate point A has a probability of 0.8 at (23,14) and candidate point B has a probability of 0.2 at (23,15), the final coordinates will fall at (23,14.2), which is closer to the actual physical source location than simply outputting (23,14). This effectively eliminates the quantization error caused by a single grid resolution, transforming discrete classification results into continuous regression results, and achieving precise location of pollution sources in chemical industrial parks. The final output... This is the final conclusion of this source tracing mission, which can be directly used to guide drone patrols or manual investigations.

[0040] Regarding the weighted centroid localization mechanism mentioned in the above embodiments, when dealing with pollution source tracing problems in complex wind field environments of chemical industrial parks, the traditional simple weighted average algorithm implicitly assumes that the spatial positioning error is uniformly distributed across all dimensions, thus ignoring the extremely special wind field constraints during atmospheric transport. In actual fluid dynamics scenarios, the diffusion scale of pollutants in the downwind direction is much more affected by advection than by turbulent diffusion in the crosswind direction, exhibiting strong anisotropy. This leads to the predicted coordinates being easily affected by excessive interference from lateral noise points under strong wind and stable conditions, resulting in lateral drift. The generated source tracing results often deviate from the actual physical transport path. Furthermore, the original formula only uses the magnitude of the probability values ​​for weighting, ignoring the local geometric features of the probability distribution surface. This causes the algorithm to be unable to effectively distinguish between two forms with the same probability integral but completely different physical meanings: the sharp isolated peaks representing deterministic point source leakage and the flat plateaus representing the blurred region after diffusion. To address the aforementioned shortcomings and enhance the model's ability to accurately focus on high-confidence leakage point sources amidst complex background noise and achieve sub-pixel accuracy, this scheme constructs a source localization mechanism based on the coupling of wind field turbulence and curvature tensor in step S62.

[0041] Based on this, in a specific preferred embodiment, sub-step S62, based on the revised candidate list, performs pollution source inversion on the grid coordinates of each candidate point to obtain the final pollution source coordinates, including: Based on the global wind field vector sequence derived from the dynamic adjacency matrix sequence, the wind field consistency factor and the wind field principal axis rotation matrix are determined. This step aims to establish a reference coordinate system that conforms to hydrodynamic characteristics, distinguishing between laminar and turbulent states by quantifying the temporal consistency of airflow, thereby providing an adaptive physical benchmark for subsequent anisotropic constraints. Specifically, the global wind field vector sequence generated in step S3 is first extracted. This sequence contains the data within the time window, i.e. =12 Global average wind speed vector at each time step These vectors are used to calculate the arithmetic mean wind vector within the time window. and wind field consistency factor The calculation formula is as follows: In the above formula, This represents the average wind vector within the time window; This indicates the length of the time window, which is set to 12 in this embodiment; The instantaneous wind speed vector at time t; The Euclidean norm (modulus) of a vector. The value represents the wind field consistency factor, which ranges from [0,1]. The closer the value is to 1, the more stable the wind direction (laminar flow dominates). The closer the value is to 0, the more turbulent the wind direction (turbulent flow dominates). This represents the wind field principal axis rotation matrix, used for coordinate system transformation; Represents the average wind vector The corresponding azimuth (relative to due east or the positive X-axis). For example, continuing the scenario of westerly winds (airflow eastward), if the wind direction has been extremely stable over the past hour, with the vector pointing due east for most of the time, the calculated average wind vector... =(2.1,0), its modulus is 2.1. If the average value of the instantaneous wind speed modulus is also approximately 2.1, then ≈1.0 indicates a strong steady-state wind field, where lateral diffusion should be strictly suppressed. The azimuth angle at this point... =0 degrees (mathematical coordinate system), rotation matrix It is an identity matrix.

[0042] Hessian matrix curvature analysis is performed on the probability distribution in the corrected candidate list to obtain a sharpness enhancement weight list. This step aims to utilize the second-order differential features of the probability surface to identify and enhance peak regions with point source characteristics, while suppressing flat background noise. Specifically, for each candidate point k in the corrected candidate list, its local curvature on the probability surface is extracted using the discrete Laplacian operator, and sharpness enhancement weights are constructed. The calculation formula is as follows: In the above formula, This represents the sharpness enhancement weight for the k-th candidate point; This represents the corrected original probability value output in step S61; This represents the sharpness sensitivity coefficient, a preset hyperparameter used to control the intensity of the reward for peak features. In one specific implementation, It is set to 10.0; Indicates the coordinates of the candidate point The discrete Laplacian operator value at a given point (approximately the trace of the Hessian matrix) is used to characterize the local concavity / convexity of the probability surface. Since the second derivative (curvature) at a local maximum point is negative, The operation ensures that only convex regions (with negative curvature) that exhibit local maxima are enhanced, while concave regions remain unchanged. The negative sign in the function makes the curvature more negative (i.e., the peak more sharp), and the exponent term larger, thus giving that point a higher weight. For example, if candidate point A has a probability of 0.8 and is located at a sharp peak ( =-0.5), the probability of candidate point B is also 0.8, but it is located on a gentle plateau ( =-0.01), in When the value is 10, the weight of point A will be amplified. The magnification is multiplied by a factor of 1, while point B is only magnified by a factor of 1. This forces the model to focus on point A.

[0043] Based on the wind field principal axis rotation matrix, the points in the corrected candidate list are projected onto the wind field coordinate system. A compression constraint is applied to the crosswind component based on the wind field consistency factor, and an anisotropic manifold projection integral is performed using a sharpness enhancement weight list to obtain the final pollution source coordinates. This step maps the spatial candidate points to a solution space that dynamically changes with the wind field. Anisotropic compression penalties are applied to the crosswind component using the wind field consistency factor to eliminate lateral drift and output accurate coordinates. Specifically, the implementation process is as follows: using the highest probability point as the reference origin, all candidate points are projected onto the wind field coordinate system, and the stability of the wind field is considered. The crosswind coordinates are nonlinearly compressed (the more stable the wind, the higher the lateral compression rate), while the downwind coordinates remain unchanged. Finally, a weighted integral and inverse transformation are performed on the compressed coordinates using sharpness enhancement weights. The calculation formula is as follows: In the above formula, This represents the relative coordinate vector of the k-th candidate point in the wind field coordinate system; This represents the transpose (i.e., inverse rotation) of the wind field principal axis rotation matrix. This represents the original absolute coordinate vector of the k-th candidate point; This represents the coordinate vector of the reference point with the highest probability. Represents a coordinate vector after anisotropic compression; matrix Let be an anisotropic compression tensor, where The term determines the compressibility in the crosswind direction. When In a strong, stable wind field, this coefficient approaches 0, meaning the lateral deviation will be forcibly zeroed out; when In turbulent flow, this coefficient approaches 1, degenerating into isotropic treatment. This represents the final calculated absolute coordinate vector of the pollution source; This represents summing over all candidate points. Specifically, as shown below... =0.9 (wind direction stable, eastward), highest probability point =(45100, 32100). There exists a high-probability interference point. =(45100,32150), which is located 50 meters due north of the reference point (i.e., a crosswind deviation of 50 meters). After projection =(0,50). After compression, the x-coordinate becomes 50×(1-0.9)^250=0.5 meters. This means that in the final weighted calculation, the contribution of this point to the Y-axis coordinate is compressed from 50 meters to 0.5 meters, thereby eliminating the interference of lateral noise on the final positioning result and improving the prediction result. It strictly conforms to the physical transmission path. By constructing this solution space that deforms with the wind, the technical solution eliminates errors that violate the laws of fluid mechanics at the physical level, achieving sub-pixel-level high-precision traceability.

[0044] In summary, the machine learning-based pollution source tracing method for chemical industrial parks, based on the embodiments of this application, has been clarified. Addressing the problem of traditional models failing due to the sparse, unstructured distribution of sensors and unsteady micro-meteorological interference in chemical industrial parks, this scheme abandons traditional static Euclidean distance modeling. Instead, it utilizes building outline information and real-time wind field data from GIS to perform physical mask filtering and Gaussian plume weight reconstruction on the static distance matrix, generating a dynamic adjacency matrix sequence that characterizes airflow around-the-wall effects and dynamic causal relationships. This effectively avoids the wall-penetrating fallacy and captures the transient topological structure with the wind. Based on this, a graph convolutional neural network embedding long short-term memory units is constructed, synchronously inputting the standardized sensor spatiotemporal tensor and the dynamic adjacency matrix to extract spatiotemporal causal features in a non-Euclidean space. Finally, by inverting the source probability distribution heatmap and combining it with continuity constraints for coordinate correction, accurate location of pollution sources is achieved under complex interference environments.

[0045] Figure 4 The diagram below shows a machine learning-based pollution source tracing system for chemical industrial parks according to an embodiment of this application. Figure 4As shown, the machine learning-based pollution tracing system 100 for chemical industrial parks includes: a chemical industrial park data acquisition module 110, used to acquire raw sensor data streams and geographic information system (GIS) data from the chemical industrial park; a data alignment and tensor quantization module 120, used to perform spatiotemporal alignment and sliding window tensor quantization on the raw sensor data streams based on GIS data to obtain a standardized input feature tensor and a static distance matrix; and a dynamic adjacency matrix generation module 130, used to utilize real-time wind field information in the standardized input feature tensor and building outline information in the GIS data to perform physical mask filtering and Gaussian plume weight estimation on the static distance matrix. The system obtains a dynamic adjacency matrix sequence; a spatiotemporal latent feature extraction module 140, which synchronously inputs the standardized input feature tensor and the dynamic adjacency matrix sequence into a graph convolutional neural network model embedded with long short-term memory units to extract spatiotemporal causal features to obtain spatiotemporal latent feature vectors; a source probability distribution heatmap analysis module 150, which performs pollution source probability heatmap inversion on the spatiotemporal latent feature vectors to obtain a source probability distribution heatmap; and a pollution source coordinate determination module 160, which performs coordinate correction based on continuity constraints on high-probability regions in the source probability distribution heatmap based on the node connection logic in the dynamic adjacency matrix sequence to obtain the final pollution source coordinates.

[0046] Here, those skilled in the art will understand that the specific operations of each step in the above-mentioned machine learning-based pollution source tracing system for chemical industrial parks have been referenced above. Figures 1 to 3 The description of the machine learning-based pollution source tracing method for chemical industrial parks is detailed here, and therefore, its repeated description will be omitted.

[0047] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for tracing pollution sources in chemical industrial parks based on machine learning, characterized in that, include: S1: Acquire raw sensor data streams and geographic information system data from the chemical industrial park; S2: Based on geographic information system data, perform spatiotemporal alignment and sliding window tensor quantization on the raw sensor data stream to obtain a standardized input feature tensor and a static distance matrix; S3: Using real-time wind field information in the standardized input feature tensor and building outline information in geographic information system data, physical mask filtering and Gaussian plume weight estimation are performed on the static distance matrix to obtain a dynamic adjacency matrix sequence. S4: Simultaneously input the standardized input feature tensor and the dynamic adjacency matrix sequence into the graph convolutional neural network model embedded with long short-term memory units to extract spatiotemporal causal features and obtain spatiotemporal latent feature vectors; S5: Perform pollution source probability heatmap inversion on the spatiotemporal latent feature vectors to obtain the source probability distribution heatmap; S6: Based on the node connection logic in the dynamic adjacency matrix sequence, the coordinates of high-probability regions in the source probability distribution heatmap are corrected based on continuity constraints to obtain the final pollution source coordinates.

2. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 1, characterized in that, The raw sensor data stream includes pollutant concentrations, wind speeds, and wind directions at each monitoring station; the geographic information system data includes metadata of the monitoring stations and geographic information of the park.

3. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 2, characterized in that, Step S2 includes: The planar coordinates of each monitoring station are analyzed from the geographic information system data; A static distance matrix is ​​constructed based on the planar coordinates of each monitoring station; The original sensor data stream is resampled along a unified time axis and missing values ​​are imputed to obtain an aligned sensor data stream; The wind direction angle in the aligned sensor data stream is orthogonally decomposed and the pollutant concentration is normalized by minimax. The processed multidimensional data is then sliced ​​and stacked based on a preset time window length to obtain a standardized input feature tensor.

4. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 1, characterized in that, Step S3 includes: The global dominant wind field vector sequence is extracted from the standardized input feature tensor to obtain the global wind field vector sequence at each time step; The connection between sensor nodes is detected by ray projection using building vector contours in geographic information system data, and the cosine value of wind direction deviation between sensor nodes is calculated by combining the global wind field vector sequence to obtain the physical blocking mask matrix and wind direction deviation factor sequence. Based on the wind direction deviation factor sequence and the physical blocking mask matrix, the static distance matrix is ​​processed to obtain the dynamic adjacency matrix sequence.

5. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 1, characterized in that, Step S4 includes: The spatial feature sequence tensor is obtained by performing time-step graph convolution aggregation and Laplacian normalization on the standardized input feature tensor using the dynamic adjacency matrix sequence. The spatial feature sequence tensor is flattened and input into a long short-term memory network unit to perform temporal state evolution and causal accumulation, and the hidden state vector of the last time step is extracted to obtain the temporal hidden state sequence. A fully connected layer is used to perform dimensionality reduction compression and nonlinear projection on the temporal hidden state sequence to obtain the spatiotemporal hidden feature vector.

6. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 1, characterized in that, Step S5 includes: High-dimensional spatial projection and feature decoding are performed on the spatiotemporal latent feature vectors to obtain the original response grid. Based on the Softmax function, the original response grid is probabilized and structured to obtain a heatmap of the source probability distribution.

7. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 1, characterized in that, Step S6 includes: Several peak points are selected from the source probability distribution heatmap based on an adaptive threshold, and the sum of the out-degree intensities of the sensor nodes corresponding to each peak point is calculated according to the dynamic adjacency matrix sequence to obtain a corrected candidate list. Based on the revised candidate list, pollution source inversion is performed on the grid coordinates of each candidate point to obtain the final pollution source coordinates.

8. The machine learning-based pollution source tracing method for chemical industrial parks according to claim 7, characterized in that, Based on the revised candidate list, pollution source inversion is performed on the grid coordinates of each candidate point to obtain the final pollution source coordinates, including: Based on the global wind field vector sequence derived from the dynamic adjacency matrix sequence, the wind field consistency factor and the wind field principal axis rotation matrix are determined. Hessian matrix curvature analysis was performed on the probability distribution of the modified candidate list to obtain the sharpness enhancement weight list; The points in the corrected candidate list are projected onto the wind field coordinate system based on the wind field principal axis rotation matrix. Compression constraints are applied to the crosswind component based on the wind field consistency factor, and anisotropic manifold projection integration is performed in combination with the sharpness enhancement weight list to obtain the final pollution source coordinates.

9. A pollution source tracing system for chemical industrial parks based on machine learning, characterized in that, include: The chemical industrial park data acquisition module is used to acquire raw sensor data streams and geographic information system data from the chemical industrial park. The data alignment and tensor quantization module is used to perform spatiotemporal alignment and sliding window tensor quantization on the raw sensor data stream based on geographic information system data to obtain a standardized input feature tensor and a static distance matrix. The dynamic adjacency matrix generation module is used to perform physical mask filtering and Gaussian plume weight estimation on the static distance matrix by utilizing real-time wind field information in the standardized input feature tensor and building outline information in the geographic information system data to obtain a dynamic adjacency matrix sequence. The spatiotemporal latent feature extraction module is used to synchronously input the standardized input feature tensor and the dynamic adjacency matrix sequence into the graph convolutional neural network model embedded with long short-term memory units to extract spatiotemporal causal features and obtain spatiotemporal latent feature vectors. The source probability distribution thermal analysis module is used to invert the pollution source probability thermal map from the spatiotemporal latent feature vector to obtain the source probability distribution thermal map. The pollution source coordinate determination module is used to perform coordinate correction based on continuity constraints on high-probability regions in the source probability distribution heatmap, based on the node connection logic in the dynamic adjacency matrix sequence, to obtain the final pollution source coordinates.

Images