Industrial park-oriented groundwater pollution source top-n probability tracing method

Abstract

The application discloses an industrial park-oriented groundwater pollution source Top-N probability tracing method, belongs to the technical field of groundwater pollution prevention and hydrogeological exploration, and comprises the following steps: preprocessing groundwater monitoring data and enterprise discharge information, constructing an enterprise-level candidate pollution source unit set and calculating initial prior probability; constructing a groundwater propagation graph and a candidate pollution source correlation graph to form a double-graph coupling structure; completing node feature initialization, abnormal information back propagation and posterior probability calculation based on a graph neural network to obtain Top-N sorting; and finally constructing a three-dimensional path probability cloud through reverse particle tracking, correcting the posterior probability by spatial overlap degree, and outputting the final result. The application integrates the advantages of physical mechanism and data driving, adapts to a multi-source pollution scene, improves tracing accuracy and practicability, and provides support for park pollution control.

Description

Technical Field

[0001] This invention relates to the field of groundwater pollution prevention and control and hydrogeological exploration technology, and in particular to a Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks. Background Technology

[0002] With the rapid advancement of my country's industrialization, industrial parks have become important carriers of regional economic development, but also high-risk areas for groundwater pollution. Characteristic pollutants generated by enterprises in these parks during production, storage, and discharge processes can easily enter groundwater aquifers through seepage and overflow. Groundwater pollutants are characterized by strong concealment, complex migration paths, and significant environmental lag. Especially in industrial park settings, the large number of enterprises, diverse emission behaviors, and the obvious superposition of multiple pollution sources result in significant three-dimensional heterogeneity in the diffusion of pollution within groundwater aquifers. This makes accurate identification and tracing of groundwater pollution sources a core prerequisite and key technical challenge for water environment management and pollution control in industrial parks.

[0003] Existing groundwater pollution source tracing technologies mainly fall into two categories: One category is numerical simulation inversion methods based on hydrogeological mechanisms. These methods can physically characterize groundwater flow and pollutant migration processes through groundwater dynamics and solute transport equations. However, they are highly sensitive to hydrogeological parameters such as aquifer permeability, boundary conditions, and source-sink terms. In complex hydrogeological conditions within industrial parks, parameter acquisition is difficult and uncertain, easily leading to significant deviations in inversion results. Furthermore, the inversion calculation process is complex and inefficient, making it difficult to adapt to the rapid source tracing and emergency response needs of multi-source pollution scenarios in industrial parks. The other category is data-driven analysis methods based on machine learning. These methods construct a mapping relationship between pollution characteristics and pollution sources through monitoring data, alleviating the reliance on precise hydrogeological parameters. However, existing data-driven methods generally lack effective physical constraints on groundwater dynamic processes, easily leading to source tracing results that contradict actual hydrogeological laws, failing to guarantee the physical rationality of the results, and making them difficult to directly apply to pollution supervision and enforcement scenarios in industrial parks.

[0004] Meanwhile, existing source tracing technologies are mostly developed based on the assumption of single-source pollution, making them difficult to adapt to the complex pollution scenarios of multiple enterprises and multiple emission sources overlapping in industrial parks, and unable to achieve precise pollution source location down to the enterprise level. Furthermore, existing methods are mostly based on deterministic inversion results, lacking probabilistic quantitative assessments and prioritization of pollution source contributions, and thus failing to provide a feasible decision-making basis for graded management and precise enforcement of groundwater pollution in industrial parks. Therefore, there is an urgent need to develop a technical method that combines the constraints of hydrogeological and physical mechanisms with the advantages of data-driven intelligent analysis, adapts to the multi-source pollution scenarios of industrial parks, and enables enterprise-level probabilistic and precise source tracing. Summary of the Invention

[0005] The purpose of this invention is to propose a Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks. This method integrates the advantages of groundwater hydrogeological and physical mechanisms with graph neural network data-driven approaches. It addresses the technical problems of inaccurate source location, lack of physical rationality in source tracing results, and inability to achieve enterprise-level probabilistic ranking under conditions of multiple pollution sources superimposed in industrial parks and uncertain hydrogeological parameters. This significantly improves the accuracy, stability, and engineering applicability of groundwater pollution source tracing, providing reliable technical support for the control and treatment of groundwater pollution in industrial parks.

[0006] To achieve the above objectives, this invention provides a Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks, comprising the following steps: Step S1: Construction of candidate pollution source units and initialization of prior probabilities: Extract groundwater monitoring data and enterprise emission information from the industrial park, complete standardized preprocessing and structured organization, construct a set of enterprise-level candidate pollution source units based on the spatial distribution of enterprises, hydrogeological zoning and pollutant exceedance characteristics of monitoring points, and calculate the initial prior probability of each candidate pollution source unit. Step S2: Construction of the dual-graph coupling structure: Based on the hydrogeological structure and groundwater dynamic conditions of the study area, a groundwater propagation map characterizing the groundwater migration process of pollutants is constructed; simultaneously, a candidate pollution source association map characterizing the correlation characteristics between pollution sources is constructed; and a spatial mapping relationship between the two types of graph structures is established to form a dual-graph coupling structure. Step S3: Dual-graph coupling probability inference and Top-N ranking: Construct a graph neural network model based on the dual-graph coupling structure, complete the initialization of node features in both graphs, propagate the monitoring anomaly information in the groundwater propagation graph in reverse along the groundwater flow path, complete the aggregation and update of pollution source features in the candidate pollution source association graph, calculate the response score of each candidate pollution source unit by fusing the two types of graph structure information, obtain the posterior probability after normalization, and complete the Top-N ranking of candidate pollution source units according to the posterior probability. Step S4: Physical consistency correction based on three-dimensional path constraints: Based on the groundwater velocity field in the study area, reverse particle tracking is carried out from the abnormal monitoring points to construct a three-dimensional path probability cloud of pollution migration; the spatial overlap between each candidate pollution source unit and the three-dimensional path probability cloud is calculated, and the posterior probability is weighted and corrected with the spatial overlap as a physical constraint, and the final Top-N probability ranking result of pollution sources is output.

[0007] Preferably, step S1 includes the following specific steps: Step S11: Standardize and preprocess the groundwater monitoring data of the industrial park. The monitoring data includes the spatial coordinates of the monitoring well, the observation time, the water level value and the pollutant concentration value. After data cleaning, outlier removal and missing value filling are completed, all data are uniformly mapped to the same spatial coordinate system and time series framework to establish a spatial-temporal integrated data index. Step S12: Structure and organize enterprise emission information, including enterprise site boundaries, emission locations, characteristic pollutant types, and historical emission records. Combined with hydrogeological zoning information of the study area, match the spatial location of enterprises with aquifer structure and groundwater flow field distribution. Based on the enterprise site range, divide the area into enterprise-level candidate pollution source units, forming a set of candidate pollution source units. ; in, For the set of candidate pollution source units, For the first One candidate pollution source unit, The total number of candidate pollution source units; construct a corresponding pollutant feature vector for each candidate pollution source unit, and establish the spatial correlation between the candidate pollution source units and the groundwater system; Step S13: Based on the spatial distance relationship between candidate pollution source units and abnormal monitoring points, the consistency of groundwater flow direction, and the degree of matching of pollutant characteristics, a prior probability calculation model is constructed to quantitatively evaluate each unit in the candidate pollution source unit set. After normalization, the initial prior probability distribution of each candidate pollution source unit is obtained.

[0008] Preferably, in step S13, the prior probability calculation model is a multi-factor fusion model, and the specific expression is as follows: ; in, For the first The original prior scores of each candidate pollution source unit, As an abnormal monitoring point, , , These are the weighting coefficients of each influencing factor, and ; The distance decay factor, expressed in exponential decay form, is as follows: ,in The spatial distance between the candidate pollution source unit and the abnormal monitoring point. The preset distance attenuation scale; The groundwater flow direction consistency factor is expressed as follows: ,in The angle between the direction of the line connecting the candidate pollution source unit to the abnormal monitoring point and the mainstream direction of groundwater flow; The pollutant characteristic matching factor is assigned a value based on the degree of matching between the characteristic pollutants of the candidate pollution source and the pollutants detected at the abnormal monitoring point. The value is 1.0 when there is a complete match, 0.6 when there is a partial match, and 0.1 when there is basically no match. The original prior scores are normalized to obtain the initial prior probability, expressed as: ; in, For the first The initial prior probability of each candidate pollution source unit.

[0009] Preferably, step S2 includes the following specific steps: Step S21: Discretize the study area into several hydrological unit nodes and construct a groundwater propagation map. , where the set of nodes For the location of hydrological units and monitoring wells, the edge is set. The directed edges between nodes represent the direction of groundwater flow. The weights of the directed edges comprehensively represent the spatial adjacency relationship, hydraulic connectivity relationship, hydraulic gradient, permeability coefficient and interlayer connectivity parameters between nodes, forming a directed weighted graph structure that reflects the groundwater flow process. Step S22: Construct a candidate pollution source association graph , where the set of nodes Set of candidate pollution source units Each unit in the set of edges The association edges between candidate pollution source units represent the spatial proximity relationship, drainage system relationship, or belonging relationship of the same hydrogeological unit between units; the characteristics of each candidate pollution source unit are encoded to form a structured feature representation; Step S23: Establish a mapping relationship between the groundwater propagation map and the candidate pollution source association map, associate each candidate pollution source unit with the hydrological unit nodes within its influence range, and realize the coupling of groundwater propagation information and candidate pollution source information through the mapping relationship to form a dual-map coupling structure.

[0010] Preferably, step S3 includes the following specific steps: Step S31: Initialization of node features in dual-graph coupling structure: Based on the dual-graph coupling structure, feature extraction and encoding are performed on the hydrological unit nodes of the groundwater propagation graph and the pollution source nodes of the candidate pollution source association graph, respectively, to complete the feature initialization of the two types of nodes and build standardized input features for the graph neural network model; Step S32: Graph Neural Network Feature Propagation and Update: Based on the constructed graph neural network model, the node features in the dual-graph coupling structure are aggregated and propagated and updated through the topological connection relationship and weighted edges between nodes, which characterizes the spatial propagation process of pollutants in the groundwater system and updates the hidden feature representation of each node. Step S33: Backpropagation of abnormal information and calculation of posterior probability: Using the pollution anomaly information of the monitoring well as the initial input signal, the backpropagation of multi-layer anomaly information is completed along the reverse path of groundwater flow in the groundwater propagation map. Through the spatial mapping relationship of the dual-map coupling structure, the anomaly response information is aggregated to the candidate pollution source node. Combined with the initial prior probability, the posterior probability of each candidate pollution source unit is calculated. Step S34: Top-N ranking output of candidate pollution sources: Sort all candidate pollution source units in descending order of posterior probability. Based on the preset output scale N, select the top N candidate pollution source units and output the Top-N ranking result of pollution sources.

[0011] Preferably, in step S31, the specific method for initializing node features is as follows: For the hydrological unit nodes in the groundwater propagation map, construct the initial feature vector. : ; in, For node spatial coordinates, The node water level value. The value represents the concentration of pollutants at the node. The hydraulic gradient at the node, The permeability coefficient of the aquifer corresponding to the node. For effective porosity, For the thickness of the aquifer, This is a vector of node attribute masks. For the pollution source nodes in the candidate pollution source association graph, construct an initial feature vector. : ; in, Spatial location of candidate pollution source units For pollutant feature vectors, As a measure of emission intensity, This refers to a priori probability values ​​or risk level parameters. By using linear transformation to map the features of the two types of nodes to a latent space of the same dimension, the standardized construction of the input features of the graph neural network is completed.

[0012] Preferably, in step S33, during the posterior probability calculation, the fusion formula between the response score and the initial prior probability is: ; in, For the first The posterior score of each candidate pollution source unit. Let be the initial prior probability. For the first The normalized response scores of each candidate pollution source unit are calculated; the posterior scores are then normalized to obtain the final posterior probability. .

[0013] Preferably, step S4 includes the following specific steps: Step S41: Based on the hydrogeological parameters and water level field data of the study area, a three-dimensional groundwater velocity field is constructed using Darcy's law to provide a hydrodynamic basis for reverse particle tracing. Step S42: Starting from the anomaly monitoring point, release virtual particles in the reverse direction in the three-dimensional flow field and carry out reverse particle tracking. Obtain the set of reverse motion trajectories of all particles through discrete iterative calculation. Step S43: Discretize the study area into a three-dimensional voxel grid, count the cumulative occurrence and residence time of particles in each voxel, construct the probability density distribution function of the pollution migration path, and obtain the three-dimensional probability cloud of the pollution migration path after spatial smoothing. Step S44: Calculate the spatial overlap between each candidate pollution source unit region and the three-dimensional path probability cloud. Use the spatial overlap as a physical consistency constraint to weight and correct the posterior probability. Re-sort according to the corrected probability and output the final Top-N probability ranking result of pollution sources.

[0014] Preferably, in step S42, the particle motion control equation for reverse particle tracking is: ; in, For the particle at time Spatial location, This represents the groundwater flow velocity at that location; the negative sign indicates backtracking in the opposite direction of the groundwater flow. The discrete formula for particle position update is: ; in, For particles in Spatial location at a given time This represents the spatial position of the particle in the previous time step. For time step, The random perturbation term, which follows a zero-mean Gaussian distribution, is used to characterize the influence of aquifer heterogeneity on particle motion. Several virtual particles are released at different vertical aquifer levels at each anomaly monitoring point, and the reverse trajectory calculation of all particles is completed by setting the maximum number of tracking steps.

[0015] Preferably, in step S43, the probability density distribution function expression of the pollution migration path is: ; in, For the particle's position in space The cumulative number of occurrences in the corresponding voxel, with the denominator being the sum of the statistical values ​​of all voxels; In step S44, the formula for calculating spatial overlap is: ; in, For the first Spatial overlap of candidate pollution source units The scope of the study area; The weighted correction formula for the posterior probability is: ; in, For the first The final probability after correction for each candidate pollution source unit For posterior probability, This represents the total number of candidate pollution source units.

[0016] Therefore, the present invention employs the above-described Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks, which has the following advantages: (1) Adapt to the needs of park supervision and achieve precise source tracing at the enterprise level: This invention constructs enterprise-level candidate pollution source units, which directly fits the engineering reality of groundwater pollution supervision in industrial parks. It can achieve precise pollution source location down to the enterprise level, solving the problem that traditional methods are based on single-source assumptions and cannot adapt to the superposition of multiple enterprises and multiple emission sources in the park. The source tracing results can directly support the park's pollution classification control and law enforcement supervision. (2) Dual-graph coupling and fusion of physical mechanism and data-driven approach to ensure the rationality of results: This invention designs a dual-graph coupling structure of groundwater propagation map and candidate pollution source association map, while taking into account the hydrogeological and physical mechanism of pollutant groundwater migration and the association characteristics between pollution sources. It embeds physical prior constraints for graph neural network modeling, solves the defects of pure data-driven method that lacks physical rationality and is prone to results that contradict hydrogeological laws, and significantly improves the reliability of source tracing results. (3) Probabilistic quantitative assessment and Top-N ranking provide clear control priorities: This invention can effectively characterize the pollution contribution of each pollution source under the superposition of multiple pollution sources through the reverse propagation of abnormal information and the aggregation of multi-dimensional features, realize the probabilistic quantitative assessment of pollution sources, and provide clear priority basis for pollution control in the park through Top-N ranking, which is suitable for the dual needs of emergency response and routine supervision. (4) Three-dimensional path physical constraint correction to improve the stability and accuracy of source tracing: This invention introduces a physical consistency correction link based on reverse particle tracking. The three-dimensional path probability cloud of pollution migration is used as a constraint to correct the posterior probability obtained by data driving, which further improves the fit between the source tracing results and the actual flow law of groundwater, reduces the source tracing error caused by the uncertainty of hydrogeological parameters and insufficient monitoring data, and improves the environmental adaptability of the method. (5) Integrated design of the whole process, taking into account both efficiency and accuracy: This invention realizes the integrated source tracing of the whole process of "prior probability initialization - dual-graph coupling probability inference - physical constraint correction", which can complete stable source tracing without a large number of high-precision hydrogeological parameters. It has high computational efficiency and strong adaptability, and can be widely used in the daily supervision of groundwater pollution and emergency source tracing of pollution incidents in various industrial parks.

[0017] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0018] Figure 1 This is a flowchart of a Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks, as described in an embodiment of the present invention. Figure 2 This is a diagram showing the park layout and water head field in an embodiment of the present invention; Figure 3 This is a comparison chart of the prior probabilities of candidate pollution sources in embodiments of the present invention; Figure 4 This is a backpropagation response intensity field diagram in an embodiment of the present invention; Figure 5 This is a posterior probability ranking diagram of TopN pollution sources in an embodiment of the present invention; Figure 6 This is a diagram showing the distribution of groundwater velocity field and permeability coefficient in an embodiment of the present invention; Figure 7 This is a reverse particle tracing path diagram of the monitoring well in an embodiment of the present invention; Figure 8 This is a three-dimensional path probability cloud plane projection diagram in an embodiment of the present invention; Figure 9 This is a probability comparison diagram before and after physical consistency correction in an embodiment of the present invention; Figure 10This is a schematic diagram of the reverse tracing three-dimensional path in an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0020] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0021] Example like Figure 1-10 As shown in the figure, this embodiment proposes a Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks. The specific implementation steps are as follows: S1: Candidate pollution source construction and prior probability initialization: S11: Preprocess the groundwater monitoring data of the industrial park. The monitoring data includes information such as the coordinates of the monitoring well points, observation time, water level and pollutant concentration. Clean the raw data, identify outliers and fill in missing values, and organize it into a structured dataset according to the standard format.

[0022] Data standardization and spatiotemporal alignment are crucial for groundwater monitoring data processing. The coordinates of monitoring well locations, water level observation sequences, and pollutant concentration observation sequences are standardized using a unified format. Data from different sources and time scales are also standardized. Outliers are identified and removed using statistical thresholds, time series fluctuation characteristics, or comparisons with neighboring wells. Missing values ​​are filled using interpolation between adjacent times, compensation from neighboring wells, or fitting to historical trends. After processing, all monitoring data is mapped to the same coordinate reference system and time step, constructing a spatial-temporal integrated data index structure that includes the correspondence between "monitoring well—time—observation value" to improve data consistency and usability for subsequent source tracing analysis.

[0023] S12: Structure and organize the emission information of enterprises in the industrial park, and combine it with the hydrogeological zoning information of the modeling area to conduct correlation analysis between the spatial location of enterprises and the aquifer structure, groundwater flow field distribution and monitoring well spatial layout; divide candidate pollution source units based on the enterprise's land area, emission location or potential impact area to form a set of candidate pollution source units.

[0024] To improve the ability of candidate pollution source sets to represent the actual pollution discharge patterns of industrial parks, the construction of candidate pollution source units is based not only on enterprise directory information but also on detailed subdivision considering the enterprise site scope, discharge location characteristics, and hydrogeological conditions of the modeling area. Therefore, the following construction method is introduced: 1. Delineate candidate pollution source units based on enterprise site boundaries: Use the site boundaries of enterprises within the industrial park as the basic spatial units. For enterprises with clearly defined emission areas, candidate pollution source units can be directly delineated based on enterprise sites. For enterprises with large land areas and multiple potential emission areas within them, further local subdivisions can be made within the enterprise sites to avoid simplifying potentially polluting locations with significant spatial differences into a single source point.

[0025] 2. Spatial matching by combining aquifer structure and groundwater flow field: After the candidate pollution source units are divided, their spatial location is matched with the aquifer layer structure, groundwater flow direction and hydraulic gradient distribution, so that the spatial expression of the candidate sources not only reflects the distribution characteristics of surface enterprises, but also takes into account the actual hydrodynamic connection in the groundwater system.

[0026] 3. Construct pollutant feature vectors for candidate pollution source units: For each candidate pollution source unit, extract the main pollutant types, emission characteristics and spatial location attributes of the corresponding enterprise to form a pollutant feature vector, which is used to characterize the potential impact of the unit on groundwater monitoring anomalies and to provide a basis for subsequent correlation analysis between candidate sources and monitoring points.

[0027] 4. Establish the relationship between the candidate pollution source unit set and spatial index: Assign a unique identifier to each candidate pollution source unit after segmentation, and bind it to its spatial location, the enterprise to which it belongs, and the characteristics of the pollutants to form a candidate pollution source unit set: ; This facilitates unified use in subsequent probability modeling, sorting output, and path analysis processes.

[0028] S13: Based on the spatial distance relationship between candidate pollution source units and abnormal monitoring points, the consistency of groundwater flow direction, and the degree of matching of pollutant type or concentration characteristics, a prior probability calculation model for candidate pollution source units is constructed. Each unit in the candidate pollution source unit set R is quantitatively evaluated to obtain the initial prior probability distribution of each candidate pollution source unit, providing initial conditions for subsequent posterior probability updates and Top-N ranking analysis of pollution sources.

[0029] Prior probability modeling and multi-factor constraint mechanism: To ensure that the initial prior probability of candidate pollution source units better reflects the groundwater pollution migration patterns in industrial parks, the setting of prior probabilities not only considers the spatial proximity between candidate pollution source units and abnormal monitoring points, but also needs to be comprehensively evaluated in conjunction with groundwater flow characteristics and the consistency of pollutant properties. Therefore, the following modeling method is introduced: 1. Probability decay modeling based on spatial distance: The spatial distance between candidate pollution source units and abnormal monitoring points is used as an important factor in the prior probability calculation. The closer the candidate source is to the monitoring point, the higher its probability of contributing to abnormal pollution. By introducing the distance decay relationship, the weight of distant candidate sources is reduced, thereby reflecting the spatial decay effect of pollutants during the migration process.

[0030] 2. Integrating hydrodynamic constraints of groundwater flow direction: Combining groundwater flow field direction and hydraulic gradient information, the hydrodynamic connection between candidate pollution source units and abnormal monitoring points is determined; for candidate sources located upstream of the main groundwater flow direction and with hydraulic connectivity, their prior weight is increased; for candidate sources that do not meet the groundwater propagation path, their initial probability is reduced to enhance the physical rationality of the prior distribution.

[0031] 3. Matching constraints based on pollutant characteristics: The pollutant types emitted by enterprises corresponding to candidate pollution source units are compared and analyzed with the pollutant characteristics detected by abnormal monitoring points; when the two have a high degree of consistency in pollutant types, combination characteristics or dominant pollutant components, the prior probability of the corresponding candidate source is increased, thereby enhancing the pollution source identification capability from the attribute level.

[0032] 4. Prior probability calculation based on multi-factor fusion: Based on factors such as spatial distance, hydrodynamic constraints, and pollutant characteristic matching, the prior probability calculation of the candidate pollution source unit set is performed. Each unit is quantitatively evaluated to form a multi-factor fusion prior probability distribution, which serves as the initial input for subsequent posterior probability updates and Top-N ranking analysis.

[0033] 5. Prior probability normalization and stability processing: The initial probability of each candidate pollution source unit is normalized to meet the probability distribution constraints, and extreme probability values ​​are smoothed to avoid numerical bias problems in subsequent probability propagation and improve the overall computational stability.

[0034] S2: Construction of Groundwater Propagation Map and Candidate Source Association Map: To improve the structural representation of groundwater pollution propagation processes and the accuracy of pollution source identification, the graph structure construction process not only needs to depict the spatial propagation relationships within the groundwater system but also needs to explicitly establish the association mechanism between candidate pollution sources and monitoring points. Therefore, the following construction method is introduced: S21: Based on groundwater monitoring data and hydrogeological information, the modeling area is discretized into multiple hydrological unit nodes, and the connection relationship between the nodes is constructed according to the groundwater flow direction, hydraulic gradient and spatial adjacency relationship to form a groundwater propagation map, which is used to characterize the migration path of pollutants in the groundwater system.

[0035] Construction of Groundwater Propagation Map: 1. Node Construction Based on Hydrological Units: The modeling region is divided into regular grids or irregular hydrological units, with each unit serving as a graph node. Node attributes include water level, pollutant concentration, and relevant hydrogeological parameters, thus achieving a discretized representation of the groundwater system. 2. Edge Connection Mechanism Based on Hydrodynamic Connections: Connection edges between nodes are constructed based on groundwater flow direction, hydraulic gradient, and spatial adjacency, enabling the graph structure to reflect the main migration paths of pollutants in groundwater. Connections in the mainstream direction are given higher weights to enhance the ability to represent actual hydrodynamic processes. 3. Propagation Modeling Integrating Spatial Adjacency and Directionality: During node connection, geometric adjacency relationships and water flow direction constraints are considered simultaneously, ensuring that the graph structure possesses both spatial continuity and satisfies the physical laws of groundwater pollution migration.

[0036] S22: Based on the spatial location relationship and hydrodynamic connection between the candidate pollution source unit set R and the groundwater monitoring point, construct the association relationship between the candidate pollution source unit and the monitoring node to form a candidate pollution source association diagram, which is used to describe the influence relationship between potential pollution sources and observed anomalies.

[0037] Construction of the candidate pollution source association map: 1. Establishment of source-monitoring point spatial association: Based on the spatial distance relationship between candidate pollution source units and groundwater monitoring points, an initial connection relationship is established between candidate source nodes and monitoring nodes to characterize the potential impact range. 2. Association screening based on hydrodynamic path constraints: Combining groundwater flow direction information, the connection relationship between candidate sources and monitoring points is screened, retaining only connections with hydrodynamic accessibility to avoid establishing invalid associations that do not conform to actual propagation paths. 3. Weighted expression of association strength: Based on spatial distance, flow direction consistency, and pollutant matching degree, different weights are assigned to the connection relationship between candidate sources and monitoring points to reflect the potential impact strength of different candidate sources on monitoring anomalies.

[0038] S23: The groundwater propagation map is fused with the candidate pollution source association map. By establishing the connection relationship between the candidate pollution source nodes and the groundwater nodes, a unified dual-graph coupling structure is formed, which provides a foundation for subsequent probability propagation and source tracing analysis based on graph neural networks.

[0039] The construction method of the dual-graph coupling structure is as follows: 1. Coupling connection between source nodes and groundwater nodes: Based on the groundwater propagation graph, candidate pollution source nodes are introduced and connected to adjacent or hydraulically associated groundwater nodes to realize the input representation of pollution sources to the groundwater system. 2. Unified graph structure representation mechanism: The groundwater propagation graph and the candidate pollution source association graph are integrated into a unified graph structure, so that the graph simultaneously contains groundwater nodes and pollution source nodes, and information interaction is realized through edge relationships. 3. Probability propagation-oriented structural design: In the dual-graph coupling structure, an information propagation path from monitoring nodes to candidate pollution source nodes is constructed, providing a structural foundation for subsequent backpropagation of abnormal information and posterior probability calculation. 4. Differentiated representation of multiple types of nodes and edges: Groundwater nodes and pollution source nodes are distinguished in the graph structure, and different types of edges (such as spatial adjacency edges and source-receptor association edges) are labeled to enhance the graph neural network's ability to model different relationships.

[0040] S3: Probabilistic source tracing calculation based on graph neural networks: S31: Based on the dual-graph coupling structure, feature initialization is performed on groundwater nodes and candidate pollution source nodes. The features include water level, pollutant concentration, spatial location and pollution source attribute information, which provide input features for graph neural networks.

[0041] To ensure that the dual-graph coupling structure can accurately represent the state of the groundwater system and the attributes of candidate pollution sources, during the graph structure feature initialization process, it is necessary not only to assign attribute features to groundwater nodes and candidate pollution source nodes respectively, but also to establish a unified feature representation space to support subsequent propagation, aggregation, and probability mapping calculations in the graph neural network. The specific explanation is as follows: Graph structure feature initialization: 1. Groundwater node feature construction: For the first node in the groundwater propagation map... groundwater nodes Construct its initial feature vector for: ; in, For node spatial coordinates, This is the water level value. This represents the pollutant concentration value. The hydraulic gradient at the node, Permeability coefficient, For effective porosity, For the thickness of the aquifer, This is a boundary type or monitoring attribute mask vector. These features simultaneously characterize the spatial location, hydrodynamic state, and pollution status of a node.

[0042] For groundwater nodes with time-series observation data, historical time window features can be further introduced to form an extended input. : ; in, and They represent the preceding Historical sequences of water levels and pollutant concentrations at each time step are used to enhance the model's ability to dynamically characterize the pollution evolution process.

[0043] 2. Construction of candidate pollution source node features: For the set of candidate pollution source units... The Middle Candidate pollution source nodes Construct its source node feature vector for: ; in, The spatial location of the candidate pollution source unit. This is a pollutant feature vector used to characterize the type, combination, or dominant component of emitted pollutants. As a measure of emission intensity, This refers to the prior probability value or risk level parameter.

[0044] Prior probabilities By directly embedding the source node features, we obtain: ; This allows the subsequent posterior probability update process to directly absorb prior information during the graph propagation stage.

[0045] 3. Unified Mapping of Node Features: Since groundwater nodes and candidate pollution source nodes have different origins and attribute dimensions, to facilitate unified graph neural network calculations, they are mapped to a unified latent space through linear transformation or multilayer perceptron: ; ; in, Indicates groundwater nodes The initial hidden representation in a graph neural network. Indicates candidate pollution source nodes The initial latent feature vector; and The feature encoding functions for groundwater nodes and candidate pollution source nodes are represented respectively, with the output dimension unified as follows: Hidden vectors are used to ensure that nodes of different types can exchange messages and interact with features within the same graph structure.

[0046] 4. Edge Feature Initialization: In addition to node features, edge features are assigned to the edges in the graph structure to enhance the expression of the strength of physical connections during the propagation process. For nodes... With nodes The edge between Its edge characteristics can be represented as: ; in, The distance between nodes. For water level difference, The angle between the side direction and the main flow direction of groundwater. These are the basic edge weights. By introducing edge features, the subsequent propagation process can be made more consistent with the migration mechanism of groundwater pollution.

[0047] S32: Based on graph neural networks, node features are propagated and aggregated. Spatial information is transmitted through the connection relationship between nodes, and the feature representation of each node is updated to characterize the propagation process of pollutants in the groundwater system.

[0048] Graph Structure Information Propagation and Feature Update: To characterize the propagation process of pollutants in the groundwater system, a neighborhood aggregation-based graph neural network is used to perform multi-layer propagation and update of node features in a coupled bi-graph structure. During propagation, not only are the topological connections between nodes considered, but edge weights, hydrodynamic direction, and source-receiver relationships are also introduced to modulate the information transmission intensity. Details are as follows: 1. Neighborhood message passing mechanism: Let the first... Layer Time Node The hidden state is Its neighborhood set is Then the node Aggregated messages at the next level It can be represented as: ; in, For the first The message transformation matrix of the layer, For the edge The corresponding propagation weight is used to measure the node's propagation weight. For nodes The intensity of the impact.

[0049] Node feature updates can be represented as: ; in, This is the node's own state transformation matrix. It is a non-linear activation function.

[0050] 2. Calculation of edge weights incorporating hydrodynamic constraints: To ensure that the information propagation process conforms to the migration patterns of groundwater pollution, propagation weights are constructed based on the hydrodynamic connections between nodes. One form is: ; in, Let be the distance decay function. It is a function of water level difference. For the flow consistency function, These are learnable or preset weight coefficients.

[0051] The above functions can be taken as follows: ; Therefore, when the distance between nodes is relatively short, the water level difference conforms to the upstream-to-downstream propagation trend and the direction is highly consistent, the edge weight is greater.

[0052] 3. Heterogeneous Propagation Update of Multiple Edge Types: For different types of edges in a dual-graph coupling structure, such as spatial adjacency edges between groundwater nodes and source-receptor association edges between candidate pollution source nodes and groundwater nodes, feature transformation is performed separately using heterogeneous graph propagation, i.e.: ; in, Represents a set of edge types. Represents a node In relation types The neighborhood set below, This represents the transformation matrix corresponding to the relationship type. By distinguishing different types of edges, the model's ability to model complex source-receptor relationships can be enhanced.

[0053] 4. Multi-layered propagation and long-range impact characterization: After multi-layered propagation, candidate pollution source nodes can receive information from multiple monitoring nodes and intermediate groundwater nodes, thereby characterizing the long-range diffusion impact of pollutants in the groundwater system. If the total number of propagation layers is... Then the first Layer node status Included - The comprehensive information within the skip neighborhood can be used for subsequent posterior probability calculations.

[0054] To avoid oversmoothing due to deep propagation, residual connections are introduced during the propagation process: ; In order to preserve the original feature information of the node itself.

[0055] S33: Taking the abnormal pollution information of the monitoring nodes as input, the graph structure is used to propagate the information backward, and the observed anomalies are gradually transmitted to the candidate pollution source nodes. The posterior probability is calculated by combining the node characteristics.

[0056] Anomaly Information Backpropagation and Posterior Probability Calculation: To achieve probabilistic source tracing by inferring potential pollution sources from monitored anomalies, an anomaly information backpropagation mechanism is constructed on a dual-graph coupling structure. This mechanism uses anomaly monitoring nodes as information sources, propagating anomalous pollution signals upstream and to candidate pollution source nodes along the groundwater propagation path and source-receptor association edges. The posterior probability of candidate pollution sources is then calculated using a graph neural network. Details are as follows: 1. Definition of abnormal monitoring signals: For monitoring nodes Anomaly signals are constructed based on the degree of deviation between the observed concentration and the background concentration. ; in, To monitor the actual pollutant concentration at the monitoring nodes, Background concentration, The standard deviation of concentration, To prevent the use of tiny constants with a denominator of zero. If If the threshold is exceeded, the node can be... It has been identified as an abnormal monitoring node.

[0057] abnormal node set Recorded as: ; in, This is the threshold for anomaly detection.

[0058] 2. The backpropagation process of abnormal information: targeting the set of abnormal monitoring nodes. The abnormal intensity is used as the backpropagation source signal and propagated back to the candidate pollution source node along the reachable path in the figure. For the candidate pollution source node... The reverse contribution it receives from the abnormal node It can be represented as: ; in, Indicates from abnormal nodes To candidate pollution source nodes The set of reachable paths For the path The backpropagation weights.

[0059] It can be determined by the distance, consistency of water flow direction, and edge type: ; in, Let be the side length. The angle between the edge direction and the reverse tracing direction. For edge-type modulation coefficients, To adjust the parameters.

[0060] The above expression shows that the shorter the distance, the higher the directional consistency, and the more the edge type conforms to the pollution migration inversion logic, the greater its backpropagation contribution.

[0061] 3. Fusion of prior probabilities and backpropagation results: After obtaining the backpropagation contributions of candidate pollution source nodes... Then, it is compared with the initial prior probability. The data is then fused to form a posterior score: ; in, is the balance coefficient between prior information and graph propagation information. When When the value is large, more emphasis is placed on prior probability; when When the size is smaller, more emphasis is placed on the graph propagation contribution obtained from the inversion of monitored anomalies.

[0062] Multiplicative fusion can also be used: ; Both of the above methods can be used to achieve joint updates of prior and posterior information.

[0063] S34: Sort the candidate pollution source nodes according to their posterior probability values ​​to obtain the Top-N pollution source identification results.

[0064] Top-N Pollution Source Ranking Output: After calculating the posterior probability of candidate pollution sources, to facilitate practical engineering applications and pollution source tracing decisions, the candidate pollution source units are ranked and output to form the Top-N pollution source identification results. Specific explanations are as follows: 1. Posterior probability ranking: Ranking the set of candidate pollution source units. The posterior probabilities of each node Arrange in descending order to obtain the sorted sequence: ; in, Indicates the sorted order of the first... The posterior probability of each candidate pollution source. The corresponding ranking result is denoted as: ; 2. Top-N Candidate Pollution Source Extraction: Based on the preset output scale Select the first from the sorted sequence These candidate pollution sources constitute the Top-N pollution source set: ; in, It can be set according to the size of the industrial park, regulatory needs, or traceability accuracy requirements.

[0065] 3. Cumulative Probability Discrimination Mechanism: A cumulative probability discrimination method can also be introduced to enhance the interpretability of Top-N results. (Before definition) The cumulative probability of each candidate pollution source for: ; when Reaching the preset threshold In this case, the Top-N set can be considered to cover most potential pollution source information. This method helps to adaptively determine the number of candidate pollution sources to be output based on actual needs.

[0066] 4. Results Output and Subsequent Path Analysis Integration: The ranking output not only includes the candidate pollution source unit number and its posterior probability value, but also simultaneously outputs its spatial location, pollutant type, and affiliated enterprise information, facilitating subsequent 3D path backtracking analysis. The Top-N pollution source set... As the key target for the next step of reverse particle tracing and 3D path probability cloud analysis, this will achieve an integrated connection between the probability tracing results and the path analysis process.

[0067] S4: 3D path analysis and probability correction based on reverse particle tracking: S41: Groundwater velocity field construction: Based on the output results of the groundwater numerical model, the three-dimensional velocity field of each grid node in the study area is calculated to provide a dynamic basis for the analysis of pollutant migration paths.

[0068] To accurately characterize the migration direction and velocity distribution of pollutants in the groundwater system during path inversion, a three-dimensional groundwater velocity field needs to be constructed based on the numerical model output. This velocity field is jointly determined by the water level field and hydrogeological parameters, and serves as the fundamental dynamic condition for subsequent inverse particle tracing.

[0069] In this embodiment, based on the water level distribution results calculated by the groundwater numerical model, the hydraulic gradient of each grid node is solved, and the groundwater flow velocity is calculated in conjunction with the aquifer permeability characteristics. The flow of groundwater in porous media follows Darcy's law, and its flow velocity is expressed as: ; in, This represents the three-dimensional velocity vector of groundwater. Permeability coefficient, For effective porosity, This represents the hydraulic gradient. This expression reflects the fundamental physical law governing the flow of groundwater along the direction of decreasing water level.

[0070] In the specific implementation, the finite difference method is used to spatially discretize the water level field, obtaining the gradient components of each grid cell in three directions, thereby constructing a continuous three-dimensional velocity field. Simultaneously, in regions with distinct aquifer stratification, the velocity is calculated separately for different strata and coupled to ensure the physical consistency of the velocity field in the vertical direction.

[0071] Through the above processing, a three-dimensional groundwater velocity field covering the entire modeling area can be obtained, providing a stable and reliable physical basis for pollution path inversion.

[0072] S42: Reverse particle tracing path calculation: Starting from the anomaly monitoring node, reverse particle tracing is performed in the velocity field to generate possible source paths of pollutants.

[0073] After obtaining the groundwater velocity field, in order to infer the pollution source path from the monitoring anomalies, the inverse particle tracing method is used to perform inverse calculations on the pollutant migration path in the flow field.

[0074] Specifically, the monitoring node detecting abnormal pollution is used as the starting point for particle release. Virtual particles are placed at the corresponding spatial locations and tracked in the opposite direction of groundwater flow. The particle motion process can be represented as follows: ; in, For the particle at time Spatial location, This represents the groundwater flow velocity at that location. The negative sign indicates backtracking in the opposite direction of the groundwater flow.

[0075] In numerical implementation, discretizing the above continuity equation yields the particle position update formula: ; in, The time step is defined as the time step. By iteratively updating the particle positions, the trajectory of pollutants migrating upstream from the monitoring point can be obtained.

[0076] To improve the stability and coverage of path inversion, multiple particles are released at each anomaly monitoring node, and small perturbations are introduced to the initial position or velocity to simulate the impact of groundwater heterogeneity and local uncertainties on migration paths, thereby forming a set of multiple possible pollution source paths.

[0077] S43: Construction of 3D Path Probability Cloud: Based on multiple reverse particle paths and statistical spatial distribution characteristics, a 3D path probability cloud for pollutant migration is constructed.

[0078] After obtaining multiple reverse particle paths, in order to perform spatial probability characterization of the possible source regions of pollutants, statistical analysis of particle trajectories is conducted to construct a three-dimensional path probability cloud.

[0079] Specifically, the trajectory points of all particles in space are aggregated, and the frequency or density distribution of particle passages is statistically analyzed on a spatial grid to form a spatial probability field. In this embodiment, a kernel density estimation method is used to smooth the particle distribution, and its probability density is expressed as: ; in, Indicates position Path probability density at , The total number of particles, For sampling points on the particle trajectory, For kernel function, This is the bandwidth parameter.

[0080] The three-dimensional probability cloud constructed in the above manner can reflect the spatial distribution range and concentration of possible sources of pollutants, with areas of higher probability density corresponding to areas of high probability of pollution sources.

[0081] S44: Path consistency constraint and probability correction: Combining the path probability cloud and the spatial distribution of candidate pollution sources, the pollution source probabilities obtained in S3 are physically consistent and corrected to obtain the final probability result.

[0082] After obtaining the three-dimensional path probability cloud, in order to further improve the physical rationality of the pollution source identification results, the path information is fused with the candidate pollution source probability results obtained in step S3, and the original probability is corrected.

[0083] Specifically, the spatial consistency of candidate sources is quantified by calculating the degree of overlap between the spatial regions of candidate pollution sources and the path probability cloud. For the first... The consistency index of a candidate pollution source unit can be expressed as: ; in, This indicates the degree of spatial overlap between candidate pollution source units and path probability clouds. This is an indicator function for candidate pollution source regions. The study area.

[0084] Based on this, the consistency index is fused with the posterior probability obtained in step S3 to obtain the corrected pollution source probability: ; in, This represents the posterior probability of candidate pollution sources calculated by a graph neural network. For spatial path consistency weights, This is the corrected final probability.

[0085] This correction process can effectively reduce the probability of unreasonable candidate pollution sources along the groundwater flow path, and improve the spatial consistency and engineering reliability of the source tracing results.

[0086] A specific implementation process is as follows: like Figure 1-10 As shown, this embodiment takes a typical industrial park as the research object, and the park's planar area is set as follows: The study area contains a main unconfined aquifer with a burial depth ranging from [missing information]. To facilitate probabilistic source tracing analysis of groundwater pollution, the study area was discretized horizontally. There are 1 grid cell, with a cell size of 1 The vertical discretization is divided into 3 layers, therefore the entire modeling area forms a total of 3 layers. Each hydrological unit node corresponds one-to-one with a node in the groundwater propagation map.

[0087] Five candidate pollution source units from enterprises were identified within the study area, denoted as follows: The central coordinates and pollutant characteristics of the five enterprise sites are as follows: Candidate Sources Electroplating company A, coordinates Characteristic pollutants are ; Candidate Sources Chemical Company B, coordinates The characteristic pollutants are TCE and chlorinated hydrocarbons; Candidate Sources Electronics company C, coordinates Characteristic pollutants are ; Candidate Sources Chemical cleaning company D, coordinates The characteristic pollutants are TCE and PCE; Candidate Sources Coating company E, coordinates The characteristic pollutants are benzene series compounds.

[0088] Eight groundwater monitoring wells were deployed in the study area, including monitoring wells lie in An abnormally high concentration of TCE was detected in the latest monitoring, with the measured concentration being [missing information]. The regional background concentration was Therefore, the monitoring well These are considered as abnormal monitoring points and serve as the starting point for subsequent reverse source tracing analysis.

[0089] S1: Candidate pollution source unit construction and prior probability initialization: First, the monitoring well data was standardized and organized, unifying the spatial coordinate system, the observation time step to one month, and establishing a "monitoring well-time-concentration / water level" index table. Then, the enterprise emission information was structured and organized to form a set of candidate pollution source units. ; For abnormal monitoring wells The prior probabilities of candidate sources are constructed using three factors: spatial distance, consistency of hydrodynamic flow direction, and pollutant matching degree. For the first... The original prior score of each candidate source is defined as: ; in, For distance attenuation factor, As a groundwater flow direction consistency factor, As a pollutant matching factor, For the weighting coefficients, take... .

[0090] The distance decay factor adopts an exponential decay form: ; in, The planar distance from the candidate source center to the anomalous well. For distance attenuation scale, the example uses... .

[0091] The flow consistency factor is defined as: ; in, This is the angle between the direction of the line connecting the candidate source to the monitoring well and the direction of the mainstream groundwater flow. If the candidate source is located upstream of the anomalous well and in the same direction, this value is larger; if it deviates significantly from the mainstream direction, this value is smaller.

[0092] The pollutant matching factor is defined as: ; The original prior scores of the five candidate sources were calculated and are shown in Table 1 below: Table 1. Original Prior Scores of Candidate Sources

[0093] Normalizing the above scores yields the initial prior probabilities: ; The calculation results are shown in Table 2 below: Table 2. Calculation results of the initial prior probabilities of candidate sources

[0094] This shows that before graph inference is performed, candidate sources and It has already demonstrated a high degree of prior probability. This process corresponds perfectly with the logic of "constructing a prior probability model based on spatial distance relationships, consistency of groundwater flow direction, and degree of matching of pollutant characteristics".

[0095] S2: Construction of groundwater propagation map and candidate source association map: 1. Based on 6000 discretized hydrological units of the study area, a groundwater propagation map was constructed. Among them, the node set Represents the set of all hydrological unit nodes and edges. This represents the hydraulic connection between nodes. For adjacent hydrological units... and The edge weights are defined as follows: ; in, The equivalent permeability coefficient of adjacent units. For water level difference, The distance between nodes. The inter-layer connectivity coefficient. For the weight parameters, take respectively .

[0096] Simultaneously construct a candidate pollution source association map. ,in If two candidate sources are located in the same drainage system, the same hydrogeological zone, or are less than 250 m apart, then an association edge is established between them.

[0097] Further establish the mapping relationship between candidate source nodes and groundwater map nodes. For each candidate source... Select all hydrological nodes within the enterprise's land parcel and the three water-bearing units below the land parcel to form a mapping node set. .For example: Mapped to 18 groundwater nodes; Mapped to 21 groundwater nodes.

[0098] This forms a dual-graph coupling structure, which is used for subsequent backpropagation of anomaly information and calculation of candidate source response scores.

[0099] For abnormal monitoring wells Construction anomaly strength: ; in,- , - , - , - .

[0100] Then we have: ; This indicates that the well has significant anomalies.

[0101] 2. Reverse propagation in groundwater propagation maps: In groundwater propagation maps, starting from the node where the anomalous well is located, propagation proceeds in the opposite direction to the groundwater flow. Backpropagation. (Round 1) Round propagation node The response intensity is denoted as Its update formula is: ; in, For nodes The upstream neighbor set, For normalized propagation weights.

[0102] During initialization, only the abnormal well nodes have response values: ; After four rounds of backpropagation, an abnormal response field can be formed in the groundwater node space.

[0103] 3. Candidate Source Response Score Aggregation: Aggregate the node responses in the groundwater propagation graph to candidate source nodes. For candidate sources... The response score is defined as: ; in The mapping weights are assigned inversely proportional to the distance between the node and the candidate source center, and satisfy the following: ; The response scores of the five candidate sources obtained in this embodiment are shown in Table 3 below: Table 3 Response Scores of Candidate Sources

[0104] 4. Posterior probability calculation: The initial prior probability is fused with the response score to obtain the posterior score. ; in .

[0105] Table 4 is obtained based on the above formula: Table 4. Posterior scores of candidate sources

[0106] Further normalization: ; In this example Since the probability has been approximately normalized, the posterior probability can be directly used as the ranking criterion. The final Top-N ranking result is as follows: ; If take The top three candidate pollution sources are: ; in As the most likely source of pollution, This is a secondary potential source of pollution.

[0107] S4: 3D Path Probability Cloud Generation and Physical Consistency Correction: 1. Construction of the three-dimensional velocity field: Based on the groundwater level distribution and permeability field of the study area, Darcy's law was used to calculate the three-dimensional velocity field. ; in, This is the groundwater velocity vector. Permeability coefficient, For effective porosity, this embodiment takes... , This represents the hydraulic gradient.

[0108] Numerical calculations show that the main flow direction in the study area is from northwest to southeast, and the average planar velocity is approximately... The vertical flow velocity amplitude is .

[0109] 2. Reverse Particle Tracing: Using Anomaly Monitoring Wells The corresponding location is the particle release point. Particles are released in three vertical layers, 100 particles per layer, for a total of 300 reverse-direction particles. The particle position update formula is: ; in,- , The term represents a random perturbation, used to indicate the heterogeneity of the medium, and follows a zero-mean Gaussian distribution. The maximum number of tracking steps is set to 120.

[0110] After reverse tracing, the particle paths mainly converge at and It is located near the upstream area.

[0111] 3. Three-dimensional path probabilistic cloud computing: The entire study area is discretized into three-dimensional voxels consistent with the groundwater grid. The number of times a particle path passes through each voxel is statistically analyzed, and the path probability density is defined. ; in For particles in volume element The estimated number of passes or dwell time within the path is then used. A three-dimensional Gaussian kernel is further employed for smoothing to obtain a continuous distribution, forming a three-dimensional path probability cloud.

[0112] 4. Path overlap and physical consistency correction: For each candidate source Calculate its land area Spatial overlap with 3D path probability cloud: ; Table 5 is obtained from this embodiment: Table 5. Calculation results of spatial overlap of candidate sources

[0113] Spatial overlap is used as a physical consistency constraint to modify the posterior probability: ; in These are the physical constraint weights.

[0114] The final probability results after correction are shown in Table 6 below: Table 6. Probabilities of candidate sources after spatial overlap correction

[0115] The final sorting result is still: ; Among them, candidate sources The highest final probability indicates that chemical cleaning company D is the most likely source of pollution for this abnormal monitoring well contamination event. Meanwhile, the 3D path probability cloud shows that the pollution migration path mainly propagates along the northwest-southeast mainstream groundwater channel in the study area, forming a high-probability convergence zone downstream of the abnormal monitoring well. This result corroborates the results of the dual-graph coupling inference, demonstrating that the method of this invention can achieve integrated output of enterprise-level probabilistic source tracing and 3D path analysis.

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

Claims

1. A Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks, characterized in that, Includes the following steps: Step S1: Construction of candidate pollution source units and initialization of prior probabilities: Extract groundwater monitoring data and enterprise emission information from the industrial park, complete standardized preprocessing and structured organization, construct a set of enterprise-level candidate pollution source units based on the spatial distribution of enterprises, hydrogeological zoning and pollutant exceedance characteristics of monitoring points, and calculate the initial prior probability of each candidate pollution source unit. Step S2: Construction of the dual-graph coupling structure: Based on the hydrogeological structure and groundwater dynamic conditions of the study area, a groundwater propagation map characterizing the groundwater migration process of pollutants is constructed; simultaneously, a candidate pollution source association map characterizing the correlation characteristics between pollution sources is constructed; and a spatial mapping relationship between the two types of graph structures is established to form a dual-graph coupling structure. Step S3: Dual-graph coupling probability inference and Top-N ranking: Construct a graph neural network model based on the dual-graph coupling structure, complete the initialization of node features in both graphs, propagate the monitoring anomaly information in the groundwater propagation graph in reverse along the groundwater flow path, complete the aggregation and update of pollution source features in the candidate pollution source association graph, calculate the response score of each candidate pollution source unit by fusing the two types of graph structure information, obtain the posterior probability after normalization, and complete the Top-N ranking of candidate pollution source units according to the posterior probability. Step S4: Physical consistency correction based on three-dimensional path constraints: Based on the groundwater velocity field in the study area, reverse particle tracking is carried out from the abnormal monitoring points to construct a three-dimensional path probability cloud of pollution migration; the spatial overlap between each candidate pollution source unit and the three-dimensional path probability cloud is calculated, and the posterior probability is weighted and corrected with the spatial overlap as a physical constraint, and the final Top-N probability ranking result of pollution sources is output.

2. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 1, characterized in that: Step S1 includes the following specific steps: Step S11: Standardize and preprocess the groundwater monitoring data of the industrial park. The monitoring data includes the spatial coordinates of the monitoring well, the observation time, the water level value and the pollutant concentration value. After data cleaning, outlier removal and missing value filling are completed, all data are uniformly mapped to the same spatial coordinate system and time series framework to establish a spatial-temporal integrated data index. Step S12: Structure and organize enterprise emission information, including enterprise site boundaries, emission locations, characteristic pollutant types, and historical emission records. Combined with hydrogeological zoning information of the study area, match the spatial location of enterprises with aquifer structure and groundwater flow field distribution. Based on the enterprise site range, divide the area into enterprise-level candidate pollution source units, forming a set of candidate pollution source units. ； in, For the set of candidate pollution source units, For the first One candidate pollution source unit, The total number of candidate pollution source units; construct a corresponding pollutant feature vector for each candidate pollution source unit, and establish the spatial correlation between the candidate pollution source units and the groundwater system; Step S13: Based on the spatial distance relationship between candidate pollution source units and abnormal monitoring points, the consistency of groundwater flow direction, and the degree of matching of pollutant characteristics, a prior probability calculation model is constructed to quantitatively evaluate each unit in the candidate pollution source unit set. After normalization, the initial prior probability distribution of each candidate pollution source unit is obtained.

3. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 2, characterized in that: In step S13, the prior probability calculation model is a multi-factor fusion model, and the specific expression is as follows: ； in, For the first The original prior scores of each candidate pollution source unit, As an abnormal monitoring point, , , These are the weighting coefficients of each influencing factor, and ; The distance decay factor, expressed in exponential decay form, is as follows: ,in The spatial distance between the candidate pollution source unit and the abnormal monitoring point. The preset distance attenuation scale; The groundwater flow direction consistency factor is expressed as follows: ,in The angle between the direction of the line connecting the candidate pollution source unit to the abnormal monitoring point and the mainstream direction of groundwater flow; The pollutant characteristic matching factor is assigned a value based on the degree of matching between the characteristic pollutants of the candidate pollution source and the pollutants detected at the abnormal monitoring point. The value is 1.0 when there is a complete match, 0.6 when there is a partial match, and 0.1 when there is basically no match. The original prior scores are normalized to obtain the initial prior probability, expressed as: ； in, For the first The initial prior probability of each candidate pollution source unit.

4. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 1, characterized in that: Step S2 includes the following specific steps: Step S21: Discretize the study area into several hydrological unit nodes and construct a groundwater propagation map. , where the set of nodes For the location of hydrological units and monitoring wells, the edge is set. The directed edges between nodes represent the direction of groundwater flow. The weights of the directed edges comprehensively represent the spatial adjacency relationship, hydraulic connectivity relationship, hydraulic gradient, permeability coefficient and interlayer connectivity parameters between nodes, forming a directed weighted graph structure that reflects the groundwater flow process. Step S22: Construct a candidate pollution source association graph , where the set of nodes Set of candidate pollution source units Each unit in the set of edges The association edges between candidate pollution source units represent the spatial proximity relationship, drainage system relationship, or belonging relationship of the same hydrogeological unit between units; the characteristics of each candidate pollution source unit are encoded to form a structured feature representation; Step S23: Establish a mapping relationship between the groundwater propagation map and the candidate pollution source association map, associate each candidate pollution source unit with the hydrological unit nodes within its influence range, and realize the coupling of groundwater propagation information and candidate pollution source information through the mapping relationship to form a dual-map coupling structure.

5. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 1, characterized in that: Step S3 includes the following specific steps: Step S31: Initialization of node features in dual-graph coupling structure: Based on the dual-graph coupling structure, feature extraction and encoding are performed on the hydrological unit nodes of the groundwater propagation graph and the pollution source nodes of the candidate pollution source association graph, respectively, to complete the feature initialization of the two types of nodes and build standardized input features for the graph neural network model; Step S32: Graph Neural Network Feature Propagation and Update: Based on the constructed graph neural network model, the node features in the dual-graph coupling structure are aggregated and propagated and updated through the topological connection relationship and weighted edges between nodes, which characterizes the spatial propagation process of pollutants in the groundwater system and updates the hidden feature representation of each node. Step S33: Backpropagation of abnormal information and calculation of posterior probability: Using the pollution anomaly information of the monitoring well as the initial input signal, the backpropagation of multi-layer anomaly information is completed along the reverse path of groundwater flow in the groundwater propagation map. Through the spatial mapping relationship of the dual-map coupling structure, the anomaly response information is aggregated to the candidate pollution source node. Combined with the initial prior probability, the posterior probability of each candidate pollution source unit is calculated. Step S34: Top-N ranking output of candidate pollution sources: Sort all candidate pollution source units in descending order of posterior probability. Based on the preset output scale N, select the top N candidate pollution source units and output the Top-N ranking result of pollution sources.

6. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 5, characterized in that: In step S31, the specific method for initializing node features is as follows: For the hydrological unit nodes in the groundwater propagation map, construct the initial feature vector. : ； in, For node spatial coordinates, The node water level value. The value represents the concentration of pollutants at the node. The hydraulic gradient at the node, The permeability coefficient of the aquifer corresponding to the node. For effective porosity, For the thickness of the aquifer, This is a vector of node attribute masks. For the pollution source nodes in the candidate pollution source association graph, construct an initial feature vector. : ； in, Spatial location of candidate pollution source units For pollutant feature vectors, As a measure of emission intensity, This refers to a priori probability values ​​or risk level parameters. By using linear transformation to map the features of the two types of nodes to a latent space of the same dimension, the standardized construction of the input features of the graph neural network is completed.

7. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 5, characterized in that: In step S33, during the posterior probability calculation, the fusion formula between the response score and the initial prior probability is: ； in, For the first The posterior score of each candidate pollution source unit. Let be the initial prior probability. For the first The normalized response scores of each candidate pollution source unit are calculated; the posterior scores are then normalized to obtain the final posterior probability. .

8. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 1, characterized in that: Step S4 includes the following specific steps: Step S41: Based on the hydrogeological parameters and water level field data of the study area, a three-dimensional groundwater velocity field is constructed using Darcy's law to provide a hydrodynamic basis for reverse particle tracing. Step S42: Starting from the anomaly monitoring point, release virtual particles in the reverse direction in the three-dimensional flow field and carry out reverse particle tracking. Obtain the set of reverse motion trajectories of all particles through discrete iterative calculation. Step S43: Discretize the study area into a three-dimensional voxel grid, count the cumulative occurrence and residence time of particles in each voxel, construct the probability density distribution function of the pollution migration path, and obtain the three-dimensional probability cloud of the pollution migration path after spatial smoothing. Step S44: Calculate the spatial overlap between each candidate pollution source unit region and the three-dimensional path probability cloud. Use the spatial overlap as a physical consistency constraint to weight and correct the posterior probability. Re-sort according to the corrected probability and output the final Top-N probability ranking result of pollution sources.

9. The Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 8, characterized in that: In step S42, the particle motion control equations for reverse particle tracking are: ； in, For the particle at time Spatial location, This represents the groundwater flow velocity at that location; the negative sign indicates backtracking in the opposite direction of the groundwater flow. The discrete formula for particle position update is: ； in, For particles in Spatial location at a given time This represents the spatial position of the particle in the previous time step. For time step, The random perturbation term, which follows a zero-mean Gaussian distribution, is used to characterize the influence of aquifer heterogeneity on particle motion. Several virtual particles are released at different vertical aquifer levels at each anomaly monitoring point, and the reverse trajectory calculation of all particles is completed by setting the maximum number of tracking steps.

10. A Top-N probabilistic source tracing method for groundwater pollution sources in industrial parks according to claim 8, characterized in that: In step S43, the probability density distribution function of the pollution migration path is expressed as follows: ； in, For the particle's position in space The cumulative number of occurrences in the corresponding voxel, with the denominator being the sum of the statistical values ​​of all voxels; In step S44, the formula for calculating spatial overlap is: ； in, For the first Spatial overlap of candidate pollution source units The scope of the study area; The weighted correction formula for the posterior probability is: ； in, For the first The final probability after correction for each candidate pollution source unit For posterior probability, This represents the total number of candidate pollution source units.

Images