Groundwater pollution concentration prediction method of time-space diagram neural network fused with physical constraint

Abstract

The invention discloses a groundwater pollution concentration prediction method of a space-time diagram neural network fused with physical constraints, and relates to the technical field of groundwater pollutant prediction.The method comprises the steps that hydrology and pollutant data in a modeling area are preprocessed, and a three-dimensional parameter sampling point set is constructed; determining spatial distribution of hydrodynamic force and solute transport parameters through an optimization algorithm; the method comprises the following steps: firstly predicting a water level field based on a cascaded graph neural network architecture, and then performing pollutant concentration prediction by taking a water level prediction result as a feature input; a physical constraint loss function is adopted for optimization, and a prediction result is synchronously updated through a space-time attention mechanism. Therefore, by the adoption of the method, collaborative prediction of the water level and the pollutant concentration is achieved, prediction reasonability is guaranteed through physical constraints, a prediction model can be built only through a small amount of data, the calculation burden of a traditional method for building a complex mechanism model is avoided, and prediction accuracy under the condition of data scarcity is remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of groundwater pollutant prediction technology, and in particular to a groundwater pollutant concentration prediction method that incorporates a spatiotemporal graph neural network with physical constraints. Background Technology

[0002] Groundwater pollution prediction is a critical scientific problem that urgently needs to be solved in the field of water resource management. With rapid industrialization and urbanization, groundwater pollution is becoming increasingly serious. Traditional prediction methods have shown significant limitations in addressing this challenge and are unable to meet the needs of modern water resource management.

[0003] Current mainstream prediction technologies are mainly based on two approaches: physical mechanism models and data-driven models. While physical mechanism models have a clear physical basis, they require extremely high accuracy in hydrogeological parameters and boundary conditions, exhibiting significant uncertainty in data-scarce regions. Data-driven models, while reducing the requirements for data quality, often struggle to accurately depict the complex physical processes of pollutant migration. In particular, emerging methods such as deep learning still have significant shortcomings in areas such as embedding physical constraints and handling spatiotemporal heterogeneity, performing poorly in addressing the complex interactions between point source pollution (e.g., from pumping wells) and regional pollutant diffusion.

[0004] More importantly, existing forecasting technologies generally adopt a step-by-step forecasting model for water level and pollutant concentration. This fragmented forecasting approach weakens the coupling relationship between hydrodynamic processes and solute migration processes. In practical applications, this technical deficiency often manifests as forecast results violating fundamental physical laws such as the law of conservation of mass, or as a slow response in handling sudden pollution events, failing to meet the timeliness requirements of emergency decision-making.

[0005] Therefore, given the increasingly severe groundwater pollution situation, developing next-generation prediction technologies has become an urgent task. An ideal technological solution needs to overcome the limitations of traditional methods, organically integrating physical mechanisms with data-driven approaches, and establishing a collaborative prediction framework for water levels and pollutant concentrations. This is not only crucial for ensuring drinking water safety but will also provide strong technical support for sustainable water resource management. Summary of the Invention

[0006] The purpose of this invention is to provide a groundwater pollution concentration prediction method that integrates a spatiotemporal graph neural network with physical constraints, in order to solve the problems of traditional groundwater pollution prediction models, such as high data dependence, low computational efficiency, and insufficient representation of spatial heterogeneity, so as to achieve more efficient and accurate groundwater pollution prediction and provide strong support for water resource management.

[0007] To achieve the above objectives, this invention provides a groundwater pollution concentration prediction method incorporating a spatiotemporal graph neural network with physical constraints, comprising the following steps:

[0008] S1. Extract all hydrological data within the modeling area, including monitoring water level, pollutant concentration, pumping well flow rate, permeability coefficient, and boundary condition data. Organize these data into a structured dataset according to a standard format, construct a high-dimensional parameter space database dominated by physical boundaries, and use hierarchical Latin hypercube sampling to generate an initial set of parameter sampling points.

[0009] S2. Based on the hydrogeological conditions of the modeling area, the parameter value range of the hydraulic condition boundary in the modeling area is set. The spatial distance between sampling points is quantified by the spatial energy function. The parameter distribution uniformity is evaluated by the maximum and minimum distance criterion. The parameter spatial distribution is optimized by the simulated annealing algorithm to obtain the optimized high-dimensional parameter sampling point set.

[0010] S3. Based on the optimized parameter spatial distribution, a numerical model of the groundwater system is constructed. Based on the spatial distribution characteristics of hydrogeological parameters, a dynamic hydrodynamic field boundary condition system is established, and a structure consistency-driven hydrodynamic field configuration is designed. Through the boundary disturbance and response test mechanism, the grid data of the numerical simulation is output.

[0011] S4. Convert the grid data output from the numerical model of the groundwater system into graph structure data, with nodes representing hydrological units and edges representing spatial adjacency relationships. Construct a cascaded prediction network architecture that integrates physical priors, and use a spatiotemporal attention gating mechanism to collaboratively update and predict groundwater levels and pollution concentrations.

[0012] Furthermore, S1 includes the hydrogeological conditions of the comprehensive modeling area, extracts the main control parameters affecting the dynamics of pollutants affecting groundwater levels, and constructs a high-dimensional parameter spatial database of three dimensions and above. Each sampling point has a unique index ID and is bound to its spatial location and parameter combination.

[0013] Furthermore, in S2, the spatial energy function is defined as the sum of the minimum distances between all pairs of points or their reciprocals, used to evaluate the distribution of the current sampling point in three-dimensional space.

[0014] Furthermore, S2 includes a multi-modal neighborhood search strategy that uses simulated annealing to iteratively optimize the parameter space distribution. The strategy is as follows:

[0015] For the parameter sampling point set, a variety of neighborhood perturbation operations are designed, including parameter perturbation, spatial location perturbation and structural perturbation; in each iteration, different neighborhood perturbation operations are selected in rotation according to the set probability.

[0016] A dual objective function is set based on the parameter uniformity of the parameter sampling point set and the spatial energy function, and the two objective functions are fused by normalized weighting. Among them, the optional indicators for parameter uniformity include variance and entropy.

[0017] Establish adaptive sampling constraints guided by the fusion boundary structure: When there are significant changes in boundary conditions in the modeling region, set a spatial guiding function to perform weighted perturbation on the sampling points near the boundary;

[0018] In the simulated annealing process, an initial temperature is set and a nonlinear annealing strategy is adopted. The acceptance criterion adopts the Metropolis standard, and inferior solutions are accepted with a certain probability when the energy function increases.

[0019] Furthermore, in S3, a numerical model of the groundwater system is constructed, including:

[0020] The optimized high-dimensional parameter sampling point set is mapped to the model spatial domain to form a parameter field consistent with the geological stratigraphic structure;

[0021] A spatial location-related boundary parameter expression mechanism is introduced to dynamically couple spatiotemporal boundary conditions and construct a dynamic boundary condition system.

[0022] Based on the actual layered structure of the groundwater system and the connectivity of aquifers within the modeling area, vertical seepage channels are set between multiple aquifer units, and interlayer discharge coefficients and storage parameters are configured to reconstruct the complete hydrodynamic field structure.

[0023] Furthermore, in S4, the attributes of nodes include the water level, pollutant concentration, hydraulic gradient, and parameter values ​​of the corresponding hydrological unit at different times; the edge connection strategy includes both horizontal adjacent unit connections and vertical interlayer connectivity, reflecting the hydraulic connection structure of the groundwater system; the edge attributes include directionality, hydraulic gradient difference, or boundary state identifier; at the same time, the graph structure at each time moment is organized into graph sequence data to form a graph time series, and the graph structure at different times shares the topology, but the node attributes evolve and update over time.

[0024] Furthermore, S4 includes: for nodes or edges with specific boundary types, adding a boundary type encoding vector, and combining it with the attributes of the node or edge to form input features.

[0025] Furthermore, in S4, the cascaded prediction network architecture that integrates physical priors includes a water level prediction module and a pollutant prediction module;

[0026] In the water level prediction module, the input features include the three-dimensional spatial coordinates of hydrological units, key hydrogeological parameters, and historical water level observation sequences. The processing includes using a multilayer perceptron to perform deep encoding of node features and processing them through residual map convolutional blocks. Each residual map convolutional block contains two parts: neighborhood information aggregation and skip connections. At the same time, a boundary condition masking mechanism is designed for various hydrological boundary conditions.

[0027] In the pollutant prediction module, the input features include water level prediction results, which are concatenated with pollutant concentration features to form an extended feature vector. Through a gated graph convolutional network structure, the migration and diffusion patterns of pollutants in the aquifer are captured.

[0028] The water level prediction module and the pollutant prediction module share a spatiotemporal attention mechanism. The spatial attention weight quantifies the hydraulic connection strength between hydrological units and automatically identifies key spatial dependencies, while the temporal attention weight models the evolution of the state of the hydrological unit itself to capture temporal features.

[0029] Furthermore, in S4, a loss function for multi-task fusion of physical constraints is designed in the cascaded prediction network architecture that integrates physical priors, including water level prediction loss term, pollutant prediction loss term and physical constraints.

[0030] Water level prediction loss item An improved relative Huber loss function is adopted, defined as:

[0031]

[0032] In the formula, δ is the adaptive threshold parameter, and σ h h represents the standard deviation of water level observations. i Let i be the actual water level value of the i-th node. To predict water level values, N is the number of nodes;

[0033] Pollutant Predicted Loss Items The mean square error, which preserves the dimensions of the physical quantity, is defined as:

[0034]

[0035] In the formula, c i This represents the actual pollutant concentration. Predicted concentration values;

[0036] Physical constraints include mass conservation constraints. Boundary condition constraints Source and sink constraints

[0037]

[0038] In the formula, ρ is density, v is Darcy velocity, Q is source-sink term, B is the set of boundary nodes, |B| is the number of boundary nodes, and h bc Given the boundary values, n is the normal vector, and q i This is the measured flux;

[0039] The total loss function is:

[0040]

[0041] In the formula, the weighting coefficients are adopted. For k = 1, 2, 3, adaptive adjustment is performed, σ k is the standard deviation parameter, and ∈ is the smoothing coefficient.

[0042] Furthermore, S4 includes updating the node state through a spatiotemporal attention gating mechanism, as follows:

[0043]

[0044]

[0045] in,

[0046]

[0047] In the formula, || represents feature splicing. Represents spatial neighborhood information, Represents time-series historical information, z i (t) is the update gate, r i (t) represents the reset gate. As a candidate state, h i (t) represents the final state output; W z 、W r 、W h For a trainable parameter matrix, Let i be the set of spatial neighbors of node i. K is the weight matrix for spatial attention values. w For the size of the time attention window, This is the weight matrix for time attention values.

[0048] Therefore, the groundwater pollution concentration prediction method using the above-mentioned spatiotemporal graph neural network that integrates physical constraints has the following technical advantages:

[0049] (1) This invention constructs a multidimensional parameter space dominated by physical boundaries, supports the construction of an integrated graph structure of parameter-node-space, and combines hierarchical Latin hypercube sampling and simulated annealing algorithm to optimize the parameter space distribution using a dual objective function; at the same time, it designs an adaptive sampling constraint guided by the fusion boundary structure to enhance the model’s ability to learn boundary effects.

[0050] (2) This invention constructs a dynamic boundary condition system, realizes the dynamic coupling expression of spatiotemporal boundary conditions, which is more in line with the needs of graph structure models to learn the change law in the time dimension, and helps to improve the accuracy of spatiotemporal prediction.

[0051] (3) The present invention designs a cascaded graph neural network prediction architecture to ensure that the input data is highly compatible with the graph model design in terms of spatial structure, attribute dimension and temporal organization, and explicitly encodes the boundary conditions to guide the graph neural network to explicitly pay attention to the differences in physical boundaries, thereby improving the model fitting ability in the prediction of boundary regions. At the same time, a spatiotemporal collaborative state update and prediction mechanism is adopted to realize the dynamic evolution of node states and prediction output.

[0052] (4) The present invention designs a joint optimization mechanism that integrates data-driven and physical constraints, namely a joint optimization strategy that integrates multiple tasks and physical constraints, which can improve the prediction accuracy of the model while ensuring the rationality of the results in a physical sense.

[0053] 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

[0054] Figure 1 This is a flowchart of a groundwater pollution concentration prediction method that integrates a spatiotemporal graph neural network with physical constraints;

[0055] Figure 2 This is a schematic diagram of the cascaded graph neural network architecture in an embodiment of a groundwater pollution concentration prediction method that integrates physical constraints and spatiotemporal graph neural networks;

[0056] Figure 3 This is a schematic diagram of a fixed head boundary in an embodiment of a groundwater pollution concentration prediction method that integrates a spatiotemporal graph neural network with physical constraints.

[0057] Figure 4 This is a sampling location map of pumping wells in an embodiment of a groundwater pollution concentration prediction method that integrates a spatiotemporal graph neural network with physical constraints. Detailed Implementation

[0058] The present invention will be explained in more detail through the following embodiments. The purpose of disclosing the present invention is to protect all changes and modifications within the scope of the present invention. The present invention is not limited to the following embodiments.

[0059] Example 1

[0060] like Figure 1 As shown, this invention provides a method for predicting groundwater pollution concentration using a spatiotemporal graph neural network that integrates physical constraints, comprising the following steps:

[0061] S1. Data preprocessing and parameter sampling.

[0062] S1.1 Extract all hydrological data within the modeling area, including monitoring water level, pollutant concentration, pumping well flow rate, permeability coefficient, and boundary condition data, and organize them into a structured dataset according to the standard format.

[0063] The processing of multi-source heterogeneous data includes adopting a multi-source data fusion strategy to uniformly project heterogeneous data such as monitoring water level, pumping volume, and formation permeability into a unified spatial grid system, thereby avoiding interference from uneven data distribution on the modeling results.

[0064] S1.2. Based on the Latin hypercube design, an initial parameter sampling point set is generated, and a multidimensional parameter space database containing pumping volume (Q) and water level boundary (R) is established. Each sampling point has a unique index in the parameter space.

[0065] To enhance the ability of graph neural networks to perceive boundary changes and spatial heterogeneity in groundwater pollutant modeling tasks, the design of initial parameter sampling points must not only cover key hydrodynamic variables such as pumping volume (Q) and water level boundary (R), but also fully consider spatial structural features and the requirements of model graph transformation. Therefore, the following optimized sampling strategy is introduced:

[0066] 1. Construct a multidimensional parameter space dominated by physical boundaries:

[0067] Based on the comprehensive modeling of the hydrogeological conditions of the area, the main control parameters affecting the dynamics of pollutants in groundwater level are extracted, including the pumping intensity of pumping wells, the elevation of water level control boundaries, and the permeability coefficient of key geological zones, forming a high-dimensional parameter space of QRK three dimensions or above, providing physical constraints for subsequent graph nodes.

[0068] 2. Introduce a three-dimensional spatial coordinate system to assist sampling:

[0069] Each set of parameter sampling points is bound to a specific monitoring location, and its three-dimensional coordinate information (x, y, z) is obtained to ensure that each sample has a mapping relationship with a spatial node, thereby supporting the subsequent construction of an integrated "parameter-node-space" graph structure.

[0070] 3. Employing Stratified Latin Hypercube Sampling (LSS) to improve boundary coverage:

[0071] The entire parameter space is divided into regions (e.g., by boundary type, hydrological zone, etc.), and Latin hypercube sampling is performed on each subspace to ensure that the point density in areas sensitive to changes in boundary conditions is relatively higher.

[0072] 4. Incorporate spatial distance uniformity constraints to prevent overfitting:

[0073] By incorporating a minimum spatial distance criterion into the selection of sampling points, the distance between any two sample points in three-dimensional space is greater than a set threshold. This avoids local clustering that could lead to overly dense connections among some nodes in the graph, thereby enhancing the robustness of the graph structure.

[0074] 5. A unique indexing mechanism for sampling points supports graph node identification:

[0075] A unique index ID is generated for each parameter point and bound to its spatial location and parameter combination. This facilitates the correspondence between the static feature input, dynamic state update and supervision signal (such as water level) of the nodes in the subsequent graph neural network, enabling tracking and optimization throughout the modeling process.

[0076] S2, Simulated annealing parameter optimization.

[0077] S2.1. Based on the hydrogeological conditions of the modeling area, set the parameter range of the hydraulic condition boundary existing in the modeling area, such as water level boundary and pumping volume.

[0078] S2.2. The spatial distance between sampling points is quantified using an energy function, and the uniformity of parameter distribution is evaluated by maximizing the minimum distance criterion.

[0079] To improve the spatial distribution rationality of initial parameter samples during graph neural network modeling and ensure the simulation results have global representativeness of the hydrological response within the region, it is necessary to further optimize the spatial uniformity of sampling points based on the initial stratified Latin hypercube sampling. This embodiment adopts a sampling point uniformity optimization strategy that integrates spatial structure perception, as follows:

[0080] 1. Construct a spatial energy function to quantify differences in sample distribution:

[0081] The three-dimensional Euclidean distance between each sampling point is set as the basic metric, and the global energy function is defined as the sum of the minimum distances between all pairs of points or their reciprocals. This function is used to evaluate whether the current sampling distribution has local density or sparseness problems in three-dimensional space, thus providing an objective function for simulated annealing optimization.

[0082] 2. Introduce the "Maximin" criterion to improve distribution coverage:

[0083] In each sampling point perturbation iteration, we attempt to maximize the minimum distance (i.e., Euclidean distance) between any two sampling points to prevent the sampling points from being highly concentrated in space, especially in areas where the model is sensitive to local spatial perturbations (such as pumping zones and boundary transition zones), thereby improving the model's responsiveness to node state changes in graph structures.

[0084] 3. Balancing both parameter space and geographic space homogeneity:

[0085] The optimization process not only considers the uniformity of the distribution of parameter dimensions (such as Q, R, K, etc.), but also integrates the geographic spatial distribution to ensure that the hydrological units represented by the nodes in the figure have good spatial coverage, thereby avoiding the problem of some regional nodes lacking training data support.

[0086] 4. Iterative perturbation and adaptive cooling mechanism:

[0087] Simulated Annealing (SA) is employed, which involves small perturbations to the position or corresponding parameter values ​​of parameter points in each round. By controlling the perturbation amplitude and acceptance probability, the algorithm approximates the global optimum from a local search. The cooling function is set based on the energy function's rate of decline to ensure stable optimization results during the convergence phase.

[0088] 5. Spatial sensitive area weighting mechanism guides sampling adjustment direction:

[0089] For areas with large model prediction errors, higher weights are given during the sampling adjustment process (such as strengthening the sampling point spacing constraint in areas with significant errors), guiding the sampling distribution to focus on areas with modeling difficulties, thereby improving the final graph neural network's ability to predict water level changes in these areas.

[0090] S2.3. Based on the multi-mode neighborhood search strategy, the simulated annealing algorithm is used to iteratively optimize the parameter space distribution, further improving the adaptability and sensitivity of the parameter space samples in the subsequent graph neural network modeling process.

[0091] The simulated annealing algorithm optimizes the structure of the parameter point set. The optimization goal is not only uniform coverage in the parameter space, but also consistency of spatial structure and full expression of the model's sensitive regions. The main strategies are as follows:

[0092] 1. Construct a multimodal perturbation mechanism to enhance search diversity:

[0093] Multiple neighborhood perturbation operations are designed for the sampling point set, including parameter perturbation (such as adjusting the pumping rate Q, water level boundary R, ​​etc.), spatial location perturbation (fine-tuning the geographical coordinates of the sampling points), and structural perturbation (local resampling). Different perturbation modes are selected alternately with a set probability in each round of optimization, which improves the diversity of the search and avoids getting trapped in local optima.

[0094] 2. Joint optimization based on multiple objective functions:

[0095] A dual objective function is introduced to simultaneously optimize the parameter uniformity (such as variance or entropy indices) and spatial energy function (such as the maximum-minimum distance criterion) of the sample point set. The two objectives are fused through a normalized weighting method, balancing the model's ability to identify parameters in the model space with the coherence of the node spatial distribution in the graph structure during optimization.

[0096] 3. Adaptive sampling constraints guided by fusion boundary structure:

[0097] When there are significant changes in boundary conditions in the modeling region (such as river boundaries, underground barriers, or pumping area boundaries), a spatial guiding function is set to weight the sampling points near the boundary, so that they are preferentially distributed in the region where the boundary gradient changes significantly, thereby enhancing the model's ability to learn from boundary effects.

[0098] 4. Dynamic temperature control and acceptance criterion design:

[0099] In the simulated annealing process, an initial temperature is set and a nonlinear annealing strategy (such as exponential or adaptive annealing) is adopted to increase the search range in the early stage and gradually converge in the later stage. The acceptance criterion adopts the Metropolis standard, which accepts inferior solutions with a certain probability when the energy function increases, thereby improving the ability to escape local extrema.

[0100] S3. Numerical modeling of groundwater flow.

[0101] S3.1 Based on the optimized parameter space distribution, construct the MODFLOW 6 mesh model and set appropriate mesh cell sizes.

[0102] S3.2 Construct a numerical model of the groundwater system with integrated structural priors. Based on the spatial distribution characteristics of hydrogeological parameters, establish a complete hydrodynamic field boundary condition system to realize the dynamic coupling and parameterized characterization of various boundary conditions.

[0103] Based on the completion of basic mesh construction and initial boundary condition setting, the focus is on the fine construction of the numerical model of the groundwater system within the modeling area, aiming to provide high-fidelity numerical simulation output that can be mapped to graph structures for graph neural networks. This includes the following key steps:

[0104] 1. Mapping from parameter space to spatial domain:

[0105] The optimized multidimensional parameter sampling point set (including pumping volume Q, water level boundary R, ​​permeability coefficient K, etc.) is mapped into the model spatial domain to form a parameter field consistent with the geological stratigraphic structure. This mapping follows the geological structure control logic and enhances the local parameter resolution in boundary-sensitive areas, thereby ensuring the response accuracy of subsequent simulations in key areas.

[0106] 2. Construction of dynamic boundary condition system:

[0107] Unlike traditional static boundary setting methods, this approach introduces spatially relevant boundary parameter expression mechanisms (such as dynamic water level functions for river boundaries and time-varying flow sequences for pumping wells) to achieve a dynamic coupling expression of spatiotemporal boundary conditions. This method better aligns with the needs of graph structure models to learn change patterns in the time dimension, thus helping to improve the accuracy of spatiotemporal predictions.

[0108] 3. Structural consistency-driven hydrodynamic field configuration:

[0109] Based on the actual layered structure and aquifer connectivity of the groundwater system within the modeling area, vertical seepage channels are established between multiple aquifer units. Interlayer discharge coefficients and storage parameters are rationally configured to reconstruct the complete hydrodynamic field structure within the area. This configuration not only enhances the physical realism of the numerical simulation but also provides accurate structural priors for cross-layer adjacency relationships in graph neural networks.

[0110] 4. Integration of boundary disturbance and response testing mechanisms:

[0111] A control perturbation mechanism (such as boundary pumping abrupt change and hydraulic head abrupt change) is introduced into the model to simulate the response behavior of the groundwater system under unsteady conditions, and a set of perturbation response curves is output. This response set serves as an important component of the training data for the graph neural network and can be used to enhance the model's predictive robustness in complex boundary scenarios.

[0112] 5. Consistent integration with the node structure of graph neural networks:

[0113] During the model output phase, the numerical simulation results (such as time-series water level and pollutant concentration, hydraulic gradient, etc.) are reorganized according to the predetermined grid node structure to ensure that their spatial distribution characteristics and adjacency relationships are consistent with the node definitions in the graph neural network, thereby achieving seamless connection from the numerical model to the graph structure model and providing structurally complete supervisory data support for subsequent graph neural network learning.

[0114] S4, Cascaded Graph Neural Network Prediction.

[0115] S4.1. Convert the grid data output from the groundwater system numerical model into graph-structured data, where nodes represent hydrological units and edges represent spatial adjacency relationships. Ensure that the input data is highly compatible with the graph model design in terms of spatial structure, attribute dimensions, and temporal organization, as detailed below:

[0116] 1. Node definition and attribute assignment:

[0117] Each hydrological unit center in the 3D mesh model is used as a node in the graph structure. The node attributes include the water level and pollutant concentration, hydraulic gradient, and parameter values ​​(such as permeability coefficient K and pumping volume Q) of the corresponding hydrological unit at different times. Auxiliary attributes such as layer number and burial depth can also be added to provide multi-scale input features for the subsequent network coding module.

[0118] 2. Edge connection construction:

[0119] Based on the adjacency information of the grid in the numerical model, an edge structure is constructed in the graph. The edge connection strategy includes connections between horizontally adjacent cells as well as vertical interlayer connectivity, reflecting the hydraulic connection structure of the regional groundwater system. Edge attributes can include directionality, hydraulic gradient difference, or boundary state identifiers, enhancing the graph structure's ability to represent real hydrodynamic processes.

[0120] 3. Graph sequence organization in the time dimension:

[0121] According to the time steps of the model simulation, the graph structure at each time step is organized into graph sequence data, forming a spatio-temporal graph sequence. The graph at each time step shares the topological structure, but the node attributes evolve and update over time, which meets the requirements of spatio-temporal modeling in graph neural networks.

[0122] 4. Explicit coding of boundary conditions:

[0123] For nodes or edges with specific boundary types (such as constant head boundary, pumping boundary, zero flux boundary, etc.), a boundary type encoding vector is added, and this vector is combined with the node or edge attributes to form the input features. This approach guides the graph neural network to explicitly focus on the differences in physical boundaries, improving the model's fitting ability in boundary region predictions.

[0124] 5. Standardization and encapsulation of graph structure input:

[0125] The final generated graph structure data is encapsulated in a standard graph data format (such as graph data objects supported by PyTorchGeometric or DGL), and preprocessed by normalization, missing data completion, and anomaly removal to provide high-quality, structurally consistent, and temporally continuous input data for subsequent graph neural networks.

[0126] S4.2 Construct a cascaded prediction network architecture that integrates physical priors, including a water level prediction module and a pollutant prediction module, as detailed below:

[0127] Collaborative modeling of two tasks is achieved through a shared spatiotemporal attention mechanism.

[0128] 1. In the design of the water level prediction module, a multilayer perceptron is first used to deeply encode the node features. The input features include the three-dimensional spatial coordinates (x, y, z) of the hydrological unit, key hydrogeological parameters (permeability coefficient K, specific yield Sy, etc.), and historical water level observation sequences. The encoded features are processed through residual map convolutional blocks. Each residual map convolutional block contains two parts: neighborhood information aggregation and skip connections. Its mathematical expression is as follows: in, W represents the hidden state of node i in the l-th layer network. (l)The trainable parameter matrix of the l-th layer network; α ij σ represents the attention weight, used to describe the normalized association strength between nodes i and j; σ is the Swish activation function. Specifically, a boundary condition masking mechanism is designed for various hydrological boundary conditions (such as constant head boundaries, pumping well boundaries, etc.), ensuring that the prediction of boundary nodes conforms to physical laws by applying strong constraints.

[0129] The pollutant prediction module uses water level prediction results as an important input feature, which is then concatenated with pollutant concentration features to form an extended feature vector. This module employs a gated graph convolutional network (GRU) structure, and its state update formula is as follows: By selectively retaining important neighborhood information through gating mechanisms, the migration and diffusion patterns of pollutants in aquifers can be effectively captured.

[0130] The two prediction modules share a core spatiotemporal attention mechanism. Spatial attention weights. It quantifies the hydraulic connection strength between hydrological units and can automatically identify key spatial dependencies such as the influence range of pumping wells; time attention weighting. This model encodes the evolution of the hydrological unit's own state, accurately capturing temporal features such as seasonal fluctuations. The entire network uses LayerNorm for layer normalization and introduces DropEdge regularization to prevent overfitting. Gradient clipping ensures numerical stability during training.

[0131] S4.3 Design a loss function for multi-task fusion physical constraints, including the relative Huber loss for water level prediction and the dimensional preservation loss for pollutant prediction.

[0132] The joint optimization strategy for multi-task fusion of physical constraints consists of three parts: the water level prediction loss term adopts an improved relative Huber loss function, defined as follows: Where δ is the adaptive threshold parameter, σ h h represents the standard deviation of water level observations. i Let i be the actual water level value of the i-th node. To predict water level values, N represents the number of nodes; the pollutant prediction loss term uses the mean square error that preserves the physical dimensions. Where c i This represents the actual pollutant concentration. To predict concentration values, the physical constraints section includes three key components: mass conservation constraints. in The gradient symbol is ρ, the density is v, the Darcy velocity is Q, and the source and sink terms are Q, ensuring that the prediction results satisfy the continuity equation; boundary conditions are constrained. Where B is the set of boundary nodes, |B| is the number of boundary nodes, and h bc Given known boundary values, force the predicted values ​​of boundary nodes to conform to the known boundary conditions; and impose source and sink constraints. Where n is the wellbore normal vector, q i To ensure the measured flux and the reasonableness of the predicted gradient near the injection / pumping wells, the total loss function is calculated using... To achieve multi-objective optimization, the weight coefficients are adopted. 1, 2, 3 are adaptively adjusted, where σ k Here, is the standard deviation parameter, ∈ is the smoothing coefficient, and it is used in the backpropagation process. To achieve dynamic equilibrium, where The total loss function is applied to the dynamic weight parameters w k The partial derivative (i.e., gradient) of w k η is the dynamic weighting parameter, balancing the optimization intensity of different tasks, and the learning rate. This system ensures, through rigorous mathematical constraints, that the prediction results both fit the observed data and conform to the fundamental principles of groundwater dynamics.

[0133] S4.4 To achieve effective simulation and prediction of groundwater dynamic processes, a joint optimization mechanism integrating data-driven and physical constraints is designed, namely a multi-task joint optimization strategy that integrates physical constraints. This aims to improve the model's prediction accuracy while ensuring the physical rationality of the results.

[0134] A spatiotemporal collaborative state update and prediction mechanism includes constructing a spatiotemporal attention gating mechanism to achieve dynamic evolution and predictive output of node states. See also... Figure 2 This embodiment uses a multi-head attention architecture to handle spatiotemporal dependencies, where the spatial attention weights are calculated as follows:

[0135]

[0136] The time attention weights are calculated as follows:

[0137]

[0138] This mechanism can automatically identify key spatiotemporal patterns, among which Let d be the learnable parameter matrix, and d be the feature dimension.

[0139] Node state updates use physically booted GRU units:

[0140]

[0141] Among them, spatial neighborhood information The set of spatial neighbors of node i The weight matrix for spatial attention values; temporal historical information. K w The size of the time attention window represents the number of historical time steps considered. This is the weight matrix for the temporal attention values. This design retains the advantages of traditional GRU in processing time-series data, while enhancing the ability to extract spatial features through the attention mechanism.

[0142] To ensure the physical rationality of the prediction results, this embodiment also applies post-processing constraints to the output results, including forcing the predicted values ​​of boundary nodes to conform to known boundary conditions.

[0143]

[0144] Example 2

[0145] This invention employs a fusion of improved Latin hypercube sampling and adaptive simulated annealing algorithms to systematically construct 200 sets of training condition datasets with strict physical constraints, providing high-quality input data for spatiotemporal graph neural networks. In the parameter space construction stage, a Stratified Latin Hypercube Sampling (LSS) method is used, with parameters determined by pumping rates Q (50-100 m³ / s). 3 In a three-dimensional parameter space comprised of the water level boundary R (initial value 90.1-92.5m, final value 81.25-82.5m) and the permeability coefficient K (0.0001-0.0007m / s, stratified sampling), each parameter dimension was divided into 200 equally probable intervals and optimized to generate an initial parameter sample set. Spatial correlation constraints were specifically introduced during the sampling process. By maximizing the minimum distance criterion and using a parameter correlation optimization algorithm, the 200 generated initial parameter combinations simultaneously met the technical requirements of a spatial coverage index h = 0.18 (better than the traditional LHS 0.25), a Pearson correlation coefficient |ρ| < 0.1 between parameters, and a 40% increase in sampling density in the boundary region.

[0146] In the parameter optimization stage, a multi-objective simulated annealing algorithm based on physical constraints was designed. A composite objective function was constructed, comprising the sum of the reciprocals of the distances between sampling points (E_dist), the parameter correlation index (E_corr), and the boundary matching penalty term (E_bound), and adaptive weight coefficients of 0.5, 0.3, and 0.2 were assigned, respectively. The optimization process adopted a temperature adaptive control strategy: T_k = T_0·exp(-α·k / N), where the initial temperature T_0 = 1000℃, the cooling coefficient α = 0.95, and the maximum number of iterations N = 1000. The specific implementation is divided into three optimization phases: In the global exploration phase (1-300 iterations), a large neighborhood search with a 15% amplitude is implemented, focusing on optimizing the global parameter distribution; in the local optimization phase (301-700 iterations), a medium perturbation with an 8% amplitude is used, simultaneously adjusting the well location coordinates and boundary parameters, and introducing an elite retention mechanism; in the fine-tuning phase (701-1000 iterations), a small perturbation with a 3% amplitude is used, focusing on optimizing boundary-sensitive areas, especially implementing enhanced optimization for the river boundary area (model columns 35, rows 10-55). The optimization results are as follows: Figure 4 As shown.

[0147] After rigorous physical consistency checks, the optimized parameter set was used for automated modeling via the Flopy library and ultimately integrated into the MODFLOW6 model input parameters. This model system contains 5400 active elements (75 rows × 125 columns of mesh), with basic parameters set at a top elevation of 100m, a bottom elevation of 70m, and an initial head of 40m. A stable flow gradient is formed through constant head boundaries on the northwest (100m) and southeast (80m) sides. Figure 3 As shown. Specific parameter configurations include: spatial coordinates and dynamic pumping parameters for the 10 optimized well locations in the Well package; river boundary water level parameters using a linear gradient function in the RIV module; the permeability field maintaining a five-level zoning structure (K1 = 0.0004 to K5 = 0.0007 m / s) in the NPF module; and the vertical recharge rate (0.8-1.8 × 10⁻⁶ m / s) based on permeability zoning in the RCH / EVT module. -9 Key hydrogeological parameters such as m / s.

[0148] The 75×125 grid data output by MODFLOW6 was converted into a spatiotemporal graph structure. Each node contains spatial coordinates (x, y, z), permeability coefficient K, water level h, pollutant concentration c, and state features (prev_h, prev_c) for five historical time steps. Two types of edge connections were constructed: spatial adjacency edges used a 4-neighborhood + vertical connection method, and time series edges associated with the same node in consecutive time steps. Boundary conditions were explicitly encoded using a 1D mask vector, including constant head boundary identifiers, river boundary water level values, and pumping well flow parameters, ensuring complete preservation of physical boundary information.

[0149] like Figure 2 As shown, the network employs a two-level prediction architecture: the water level prediction module (HeadGNN) contains three layers of residual graph convolutions, each with a hidden dimension of 128, using the ReLU activation function; the pollutant prediction module (ConcGNN) uses two layers of gated graph convolutions (GRU units), with a hidden dimension of 64, and takes as input the stitched water level prediction results and the original node features (total dimension 192). Both modules output the final predicted value as linear layers, achieving collaborative modeling of hydrodynamics and pollutant migration through cascading connections.

[0150] The training employs a composite loss function: water level prediction uses a relative Huber loss (δ = 0.1, σ_h = 1.2m), while pollutant prediction uses a dimensionless MSE loss (concentration normalized to [0,1]). Physical constraints such as flux conservation and boundary matching are introduced, and the two prediction tasks are balanced by λ = 0.3, while the strength of the physical constraints is controlled by β = 0.5, ensuring that the prediction results simultaneously satisfy data fitting and physical laws.

[0151] A dynamic update mechanism with 4-head spatial attention and 5-step temporal attention windows is employed. Training uses the Adam optimizer (learning rate 1e-4) with a batch size of 16, reaching its optimum after 200 training rounds (early stopping strategy). The final model results on the training, validation, and test sets are shown in Table 1, where the mean square error of water level (MSE_Head) and the mean square error of concentration (MSE_Conc) are also included.

[0152] Table 1. Prediction accuracy of cascaded graph neural networks

[0153]

[0154] Therefore, the present invention adopts the above-mentioned groundwater pollution concentration prediction method based on spatiotemporal graph neural network with physical constraints, which realizes the coordinated prediction of water level and pollutant concentration. The physical constraints ensure the rationality of the prediction. Only a small amount of data is needed to build the prediction model, avoiding the computational burden of building complex mechanism models by traditional methods, and significantly improving the prediction accuracy under data-scarce conditions.

[0155] 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 groundwater pollution concentration prediction method incorporating a spatiotemporal graph neural network with physical constraints, characterized in that, Includes the following steps: S1. Extract all hydrological data within the modeling area, including monitoring water level, pollutant concentration, pumping well flow rate, permeability coefficient, and boundary condition data. Organize these data into a structured dataset according to a standard format, construct a high-dimensional parameter space database dominated by physical boundaries, and use hierarchical Latin hypercube sampling to generate an initial set of parameter sampling points. S2. Based on the hydrogeological conditions of the modeling area, the parameter value range of the hydraulic condition boundary in the modeling area is set. The spatial distance between sampling points is quantified by the spatial energy function. The parameter distribution uniformity is evaluated by the maximum and minimum distance criterion. The parameter spatial distribution is optimized by the simulated annealing algorithm to obtain the optimized high-dimensional parameter sampling point set. S3. Based on the optimized parameter spatial distribution, a numerical model of the groundwater system is constructed. Based on the spatial distribution characteristics of hydrogeological parameters, a dynamic hydrodynamic field boundary condition system is established, and a structure consistency-driven hydrodynamic field configuration is designed. Through the boundary disturbance and response test mechanism, the grid data of the numerical simulation is output. S4. Convert the grid data output from the numerical model of the groundwater system into graph structure data, with nodes representing hydrological units and edges representing spatial adjacency relationships. Construct a cascaded prediction network architecture that integrates physical priors, and use a spatiotemporal attention gating mechanism to collaboratively update and predict groundwater levels and pollution concentrations.

2. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, S1 includes the hydrogeological conditions of the comprehensive modeling area, extracts the main control parameters that affect the dynamics of pollutants affecting the groundwater level, and constructs a high-dimensional parameter spatial database of three dimensions and above. Each sampling point has a unique index ID and is bound to its spatial location and parameter combination.

3. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints as described in claim 1, characterized in that, In S2, the spatial energy function is defined as the sum of the minimum distances between all pairs of points or their reciprocals, and is used to evaluate the distribution of the current sampling point in three-dimensional space.

4. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, S2 includes a multi-modal neighborhood search strategy that uses simulated annealing to iteratively optimize the parameter space distribution. The strategy is as follows: For the parameter sampling point set, a variety of neighborhood perturbation operations are designed, including parameter perturbation, spatial location perturbation and structural perturbation; in each iteration, different neighborhood perturbation operations are selected in rotation according to the set probability. A dual objective function is set based on the parameter uniformity of the parameter sampling point set and the spatial energy function, and the two objective functions are fused by normalized weighting. Among them, the optional indicators for parameter uniformity include variance and entropy. Establish adaptive sampling constraints guided by the fusion boundary structure: When there are significant changes in boundary conditions in the modeling region, set a spatial guiding function to perform weighted perturbation on the sampling points near the boundary; In the simulated annealing process, an initial temperature is set and a nonlinear annealing strategy is adopted. The acceptance criterion adopts the Metropolis standard, and inferior solutions are accepted with a certain probability when the energy function increases.

5. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, In S3, a numerical model of the groundwater system is constructed, including: The optimized high-dimensional parameter sampling point set is mapped to the model spatial domain to form a parameter field consistent with the geological stratigraphic structure; A spatial location-related boundary parameter expression mechanism is introduced to dynamically couple spatiotemporal boundary conditions and construct a dynamic boundary condition system. Based on the actual layered structure of the groundwater system and the connectivity of aquifers within the modeling area, vertical seepage channels are set between multiple aquifer units, and interlayer discharge coefficients and storage parameters are configured to reconstruct the complete hydrodynamic field structure.

6. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, In S4, the attributes of a node include the water level, pollutant concentration, hydraulic gradient, and parameter values ​​of the corresponding hydrological unit at different times; the edge connection strategy includes both horizontal adjacent unit connections and vertical interlayer connectivity, reflecting the hydraulic connection structure of the groundwater system; the edge attributes include directionality, hydraulic gradient difference, or boundary state identifier; at the same time, the graph structure at each time moment is organized into graph sequence data to form a graph time series, and the graph structure at different times shares the topology, but the node attributes evolve and update over time.

7. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, S4 includes: for nodes or edges with specific boundary types, adding a boundary type encoding vector, and combining it with the attributes of the node or edge to form the input feature.

8. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, In S4, the cascaded prediction network architecture that integrates physical priors includes a water level prediction module and a pollutant prediction module. In the water level prediction module, the input features include the three-dimensional spatial coordinates of hydrological units, key hydrogeological parameters, and historical water level observation sequences. The processing includes using a multilayer perceptron to perform deep encoding of node features and processing them through residual map convolutional blocks. Each residual map convolutional block contains two parts: neighborhood information aggregation and skip connections. At the same time, a boundary condition masking mechanism is designed for various hydrological boundary conditions. In the pollutant prediction module, the input features include water level prediction results, which are concatenated with pollutant concentration features to form an extended feature vector. Through a gated graph convolutional network structure, the migration and diffusion patterns of pollutants in the aquifer are captured. The water level prediction module and the pollutant prediction module share a spatiotemporal attention mechanism. The spatial attention weight quantifies the hydraulic connection strength between hydrological units and automatically identifies key spatial dependencies, while the temporal attention weight models the evolution of the state of the hydrological unit itself to capture temporal features.

9. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, In S4, a loss function for multi-task fusion of physical priors is designed in the cascaded prediction network architecture, including water level prediction loss term, pollutant prediction loss term and physical constraints. Water level prediction loss item An improved relative Huber loss function is adopted, defined as: In the formula, δ is the adaptive threshold parameter, and σ h h represents the standard deviation of water level observations. i Let i be the actual water level value of the i-th node. To predict water level values, N is the number of nodes; Pollutant Predicted Loss Items The mean square error, which preserves the dimensions of the physical quantity, is defined as: In the formula, c i This represents the actual pollutant concentration. Predicted concentration values; Physical constraints include mass conservation constraints. Boundary condition constraints Source and sink constraints In the formula, ρ is density, v is Darcy velocity, Q is source-sink term, B is the set of boundary nodes, |B| is the number of boundary nodes, and h bc Given the boundary values, n is the normal vector, and q i This is the measured flux; The total loss function is: In the formula, the weighting coefficients are adopted. Perform adaptive adjustment, σ k is the standard deviation parameter, and ∈ is the smoothing coefficient.

10. The groundwater pollution concentration prediction method based on a spatiotemporal graph neural network incorporating physical constraints according to claim 1, characterized in that, S4 includes updating node states through a spatiotemporal attention gating mechanism, as follows: in, In the formula, || represents feature splicing. Represents spatial neighborhood information, Represents time-series historical information, z i (t) is the update gate, r i (t) represents the reset gate. As a candidate state, h i (t) represents the final state output; W z 、W r 、W h For a trainable parameter matrix, Let i be the set of spatial neighbors of node i. K is the weight matrix for spatial attention values. w For the size of the time attention window, This is the weight matrix for time attention values.

Images