Groundwater pollution concentration prediction method of spatio-temporal graph neural network fused with physical constraints

Abstract

The application discloses a groundwater pollution concentration prediction method of a space-time graph neural network fusing physical constraints, relates to the technical field of groundwater pollution prediction, and comprises the following steps: preprocessing hydrological and pollutant data in a modeling area and constructing a three-dimensional parameter sampling point set; determining the spatial distribution of hydrodynamic and solute transport parameters through an optimization algorithm; based on a cascaded graph neural network architecture, first predicting a water level field, and then inputting the water level prediction result as a feature to predict the pollutant concentration; and adopting a physical constraint loss function to optimize and synchronously updating the prediction result through a space-time attention mechanism. Therefore, the above method realizes the collaborative prediction of the water level and the pollutant concentration, the physical constraint guarantees the prediction rationality, the prediction model can be established only by using a small amount of data, the calculation burden of constructing a complex mechanism model in the traditional method is avoided, and the prediction accuracy under the condition of data scarcity is significantly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of groundwater pollution prediction, in particular to a groundwater pollution concentration prediction method based on a spatio-temporal graph neural network with physical constraints. BACKGROUND

[0002] Groundwater pollution prediction is a key scientific problem in water resource management that needs to be solved urgently. With the rapid development of industrialization and urbanization, groundwater pollution problems are becoming increasingly serious. Traditional prediction methods have obvious limitations in dealing with this challenge and are difficult to meet the needs of modern water resource management.

[0003] The current mainstream prediction technology is mainly based on two technical routes: physical mechanism models and data-driven models. Although the physical mechanism model has a clear physical basis, it requires high precision for hydrogeological parameters and boundary conditions, and has significant uncertainty in data-scarce areas. Although the data-driven model reduces the requirement for data quality, it is often difficult to accurately describe the complex physical process of pollutant migration. In particular, emerging methods such as deep learning still have obvious shortcomings in embedding physical constraints and handling spatio-temporal heterogeneity, and perform poorly in dealing with the complex interaction between point source pollution such as pumping wells and regional pollutant diffusion.

[0004] More importantly, existing prediction technologies generally use a step-by-step prediction mode for water level and pollutant concentration, which weakens the coupling relationship between hydrodynamic processes and solute transport processes. In practical applications, this technical defect often manifests itself in predicted results that violate basic physical laws such as conservation of mass, or in slow response when dealing with sudden pollution incidents, which cannot meet the timeliness requirements of emergency decision-making.

[0005] Therefore, in the face of the increasingly serious situation of groundwater pollution, it is urgent to develop a new generation of prediction technology. An ideal technical solution needs to break through the limitations of traditional methods, achieve the organic integration of physical mechanisms and data-driven methods, and establish a collaborative prediction framework for water level and pollutant concentration. This not only has important significance for ensuring drinking water safety, but also provides strong technical support for sustainable management of water resources. SUMMARY

[0006] The purpose of the present application is to provide a groundwater pollution concentration prediction method based on a spatio-temporal graph neural network with physical constraints, to solve the problems such as high data dependency, low computational efficiency, and insufficient representation of spatial heterogeneity in traditional groundwater pollution prediction models, and to achieve more efficient and accurate groundwater pollution prediction, providing strong support for water resource management.

[0007] To achieve the above purpose, the present application provides a groundwater pollution concentration prediction method based on a spatio-temporal 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] Adaptive sampling constraint guided by fusion of boundary structure: when there is a significant change in boundary condition in the modeling area, set a spatial guide function to weight the disturbance of the sampling points near the boundary;

[0018] In the simulated annealing process, set the initial temperature and adopt the nonlinear annealing strategy, and accept the Metropolis standard as the acceptance criterion, and accept the poor solution with a certain probability when the energy function rises.

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

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

[0021] Introduce a boundary parameter expression mechanism related to the spatial position to dynamically couple and express the spatiotemporal boundary conditions, and construct a dynamic boundary condition system;

[0022] According to the actual layered structure and aquifer connectivity of the groundwater system in the modeling area, set the vertical seepage channels between the multi-layer aquifer units, and configure the interlayer drainage coefficient and storage parameters, and reconstruct the complete hydrodynamic field structure.

[0023] Further, in S4, the attributes of the nodes include the water level values, pollutant concentrations, hydraulic gradients, and parameter values of the corresponding hydrological units at different times; the connection strategy of the edges includes both the connection between horizontally adjacent units and the connectivity between vertically adjacent layers, reflecting the hydraulic connection structure in the groundwater system; the attributes of the edges include directionality, hydraulic gradient difference, or boundary state identifier; at the same time, the graph structure at each time is organized into a graph sequence data to form a graph time sequence, and the graph structures at different times share the topological structure, but the attributes of the nodes evolve and update over time.

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

[0025] Further, 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 the hydrological units, key hydrogeological parameters, and historical water level observation sequences, and the processing process includes deep encoding of node features using a multilayer perceptron and processing through residual graph convolution blocks, each residual graph convolution block containing neighborhood information aggregation and skip connection; at the same time, a boundary condition mask mechanism is designed for various types of hydrological boundary conditions;

[0027] In the pollutant prediction module, the input features include water level prediction results, and the pollutant concentration features are spliced to form an extended feature vector, and the diffusion and migration law of the pollutant in the aquifer is captured through a gated graph convolution network structure;

[0028] The water level prediction module and the pollutant prediction module share the spatio-temporal attention mechanism, wherein the spatial attention weight quantifies the hydraulic connection strength between hydrological units, automatically identifies the key spatial dependence relationship, and the temporal attention weight models the evolution law of the state of the hydrological unit itself, and is used to capture the time sequence features.

[0029] Further, in the cascade prediction network architecture fusing physical priors in S4, a loss function fusing physical constraints for multiple tasks is designed, including a water level prediction loss term, a pollutant prediction loss term and a physical constraint;

[0030] The water level prediction loss term An improved relative Huber loss function is used, which is defined as:

[0031]

[0032] In the formula, δ is an adaptive threshold parameter, σ h is the standard deviation of the water level observation, h i is the true water level value of the i th node, is the predicted water level value, and N is the number of nodes;

[0033] The pollutant prediction loss term The mean square error of the physical dimension is used, which is defined as:

[0034]

[0035] In the formula, c i is the true pollutant concentration, and the predicted concentration value;

[0036] The physical constraint includes mass conservation constraint boundary condition constraint and source and sink term constraint

[0037]

[0038] In the formula, ρ is the density, v is the Darcy flow rate, Q is the source and sink term, B is the set of boundary nodes, |B| is the number of boundary points, h bc is the known boundary value, n is the normal vector, q i is the measured flux;

[0039] The total loss function is:

[0040]

[0041] In the formula, the weight coefficient adopts k=1, 2, 3 is adaptively adjusted, and sigma k is a standard deviation parameter, and epsilon is a smoothing coefficient.

[0042] Further, S4 comprises updating the node state through a spatio-temporal attention gating mechanism, as follows:

[0043]

[0044]

[0045] wherein,

[0046]

[0047] In the formula, || represents feature splicing, represents spatial neighborhood information, represents time series historical information, and z i (t) is an update gate, and r i (t) is a reset gate, is a candidate state, and h i (t) is a final state output; W z , W r , W h are trainable parameter matrices, is a spatial neighbor set of node i, is a weight matrix of spatial attention values, and K w is a time attention window size, is a weight matrix of time attention values.

[0048] Therefore, the groundwater pollution concentration prediction method adopting the above-mentioned spatio-temporal graph neural network fusing physical constraints has the following technical effects:

[0049] (1) The application constructs a multi-dimensional parameter space dominated by physical boundaries, supports the construction of a parameter-node-space integrated graph structure, and combines hierarchical Latin hypercube sampling and simulated annealing algorithm to optimize the parameter space distribution by using a double-objective function; meanwhile, an adaptive sampling constraint fusing boundary structure guidance is designed, which enhances the learning ability of the model to the boundary effect.

[0050] (2) The application constructs a dynamic boundary condition system, realizes dynamic coupling expression of spatio-temporal boundary conditions, is more in line with the demand of learning change law in time dimension of the graph structure model, and is helpful to improve the spatio-temporal prediction accuracy.

[0051] (3) The application designs a cascaded graph neural network prediction architecture, ensures that the input data is highly compatible with the graph model design in terms of spatial structure, attribute dimension and time sequence organization, and explicitly encodes the boundary conditions, guides the graph neural network to explicitly pay attention to the physical boundary difference, and improves the model fitting capability in the prediction of the boundary area; at the same time, a space-time collaborative state updating and prediction mechanism is adopted to realize the dynamic evolution and prediction output of the node state.

[0052] (4) The application designs a set of joint optimization mechanism fusing data driving and physical constraints, that is, a multi-task fusion physical constraint joint optimization strategy, which can improve the prediction accuracy of the model while ensuring the rationality of the results in the physical sense.

[0053] The technical solutions of the application will be further described in detail below with the help of the drawings and examples. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 is a flowchart of the groundwater pollution concentration prediction method of the space-time graph neural network fusing physical constraints;

[0055] Figure 2 is a schematic diagram of the cascaded graph neural network architecture in the groundwater pollution concentration prediction method embodiment of the space-time graph neural network fusing physical constraints;

[0056] Figure 3 is a schematic diagram of the constant head boundary in the groundwater pollution concentration prediction method embodiment of the space-time graph neural network fusing physical constraints;

[0057] Figure 4 is a pumping well sampling position map in the groundwater pollution concentration prediction method embodiment of the space-time graph neural network fusing physical constraints. DETAILED DESCRIPTION

[0058] The application can be explained in more detail through the following examples, and the purpose of disclosing the application is to protect all changes and improvements within the scope of the application, and the application is not limited to the following examples.

[0059] Example 1

[0060] As shown in Figure 1 , the application provides a groundwater pollution concentration prediction method of a space-time graph neural network fusing physical constraints, including the following steps:

[0061] S1, data preprocessing and parameter sampling.

[0062] S1.1, extract all hydrological data in the modeling area, including monitoring water level, pollutant concentration, pumping well flow, permeability coefficient and boundary condition data, and arrange them into a structured data set according to the standard format.

[0063] The processing of the multi-source heterogeneous data includes projecting monitoring water level, pumping volume, and stratum permeability and other heterogeneous data into a unified spatial grid system by using a multi-source data fusion strategy, so as to avoid the interference of uneven data distribution on the modeling result.

[0064] S1.2, generating an initial parameter sampling point set based on a Latin hypercube design, establishing a multi-dimensional parameter space database containing pumping volume (Q) and water level boundary (R), and each sampling point having a unique index in the parameter space.

[0065] In order to improve the perception ability of the graph neural network to boundary changes and spatial heterogeneity in the modeling task of groundwater pollutants, the design of the initial parameter sampling point not only needs to cover the water power key variables such as pumping volume (Q) and water level boundary (R), but also needs to fully consider the requirements of spatial structure characteristics and model graph conversion. Therefore, the following preferred sampling strategy is introduced:

[0066] 1. Constructing a multi-dimensional parameter space dominated by physical boundaries:

[0067] Based on the hydrogeological conditions of the comprehensive modeling area, the main control parameters affecting the dynamic of groundwater level pollutants are extracted, including the pumping strength of the pumping well, the elevation value of the water level control boundary, the permeability coefficient of the key geological zone, and the like, to form a high-dimensional parameter space of Q-R-K three-dimensional or more, providing physical constraint conditions for subsequent graph nodes.

[0068] 2. Introducing a three-dimensional space coordinate system to assist sampling:

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

[0070] 3. Using stratified Latin hypercube sampling (Stratified LHS) to improve boundary coverage:

[0071] The entire parameter space is divided into regions (such as by boundary type, hydrological partition, etc.), and Latin hypercube sampling is performed on each subspace to ensure that the point density in the boundary condition change sensitive area is relatively higher.

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

[0073] The spatial minimum distance criterion is added in the sampling point selection, so that the distance between any two sample points in the three-dimensional space is greater than a set threshold, avoiding local aggregation leading to excessive connection of some nodes in the graph, and enhancing the robustness of the graph structure.

[0074] 5. Sampling point unique index mechanism supports graph node identification:

[0075] A unique index ID is generated for each parameter point, which is bound with its spatial location and parameter combination, facilitating the static feature input of nodes in the subsequent graph neural network, dynamic state update, and correspondence with supervision signals (such as water level), realizing tracking and optimization in the whole modeling process.

[0076] S2, simulated annealing parameter optimization.

[0077] S2.1, according to the hydrogeological conditions of the modeling area, set the parameter value range of the existing hydraulic condition boundary in the modeling area, such as water level boundary and pumping capacity, etc.

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

[0079] In order to improve the spatial distribution rationality of initial parameter samples in the process of graph neural network modeling, and ensure that the simulation results have global representativeness of hydrological response in the region, the spatial uniformity of sampling points needs to be further optimized on the basis of 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 the difference in sample distribution:

[0081] Set the three-dimensional Euclidean distance between each sampling point as the basic measurement index, and define the global energy function as the sum or reciprocal sum of the minimum distance values between all point pairs. This function is used to evaluate whether there is a problem of local concentration or dispersion in the three-dimensional space of the current sampling distribution, thereby providing an objective function for simulated annealing optimization.

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

[0083] In each iteration of sampling point disturbance, try to maximize the minimum distance (i.e. Euclidean distance) between any two sampling points, prevent the concentration of sample points in space, especially in areas sensitive to local spatial disturbance (such as pumping area, boundary transition zone), and improve the response ability of the model to changes in node state in the graph structure.

[0084] 3. Consider both parameter space and geographical space uniformity:

[0085] In the optimization process, not only the distribution uniformity of parameter dimensions (such as Q, R, K, etc.) is considered, but also the geographical space distribution is integrated to ensure that the hydrological units represented by the nodes in the graph have good spatial coverage, thereby avoiding the problem of lack of training data support for nodes in some areas.

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

[0087] In each round, the position of the parameter point or the corresponding parameter value is slightly disturbed by using the simulated annealing algorithm (SA). By controlling the disturbance amplitude and acceptance probability, the approximation from local search to global optimum is realized. The cooling function is set according to the energy function descent rate to ensure the stability of the optimization result in the convergence stage.

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

[0089] For areas with large model prediction errors, higher weights are given in the sampling adjustment process (such as strengthening the sampling point spacing constraint in the error significant area), guiding the sampling distribution to focus on the modeling difficulty area, thereby improving the prediction ability of the final graph neural network for 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 sample in the subsequent graph neural network modeling process.

[0091] Among them, the simulated annealing algorithm optimizes the structure of the parameter point set, and the optimization target is not only the uniform coverage in the parameter space, but also the consistency of the space structure and the sufficient expression of the model sensitive area. The main strategies are as follows:

[0092] 1. Constructing a multi-modal disturbance mechanism to improve search diversity:

[0093] A variety of neighborhood disturbance operations are designed for the sampling point set, including parameter disturbance (such as adjusting the pumping quantity Q and water level boundary R value), spatial position disturbance (fine-tuning the geographic coordinates of the sampling point), and structure disturbance (local resampling). Different disturbance modes are selected with a set probability in each round of optimization to improve search diversity and avoid falling into local optimum.

[0094] 2. Joint optimization based on multi-objective function:

[0095] A double-objective function is introduced to optimize the parameter uniformity (such as variance or entropy index) and spatial energy function (such as maximum and minimum distance criterion) of the sampling point set. By normalizing and weighting the two objectives, the identification ability of the model parameter space and the continuity of the node spatial distribution in the graph structure are considered in the optimization.

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

[0097] When there are significant boundary condition changes in the modeling area (such as river boundaries, underground barriers, or pumping area boundaries), setting up a spatial guiding function to weight the disturbance of sampling points near the boundary, so that they are preferentially distributed in areas with significant boundary gradient changes, thereby enhancing the model's learning ability for boundary effects.

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

[0099] In the simulated annealing process, set the initial temperature and use a nonlinear annealing strategy (such as exponential or adaptive annealing), which can expand the search range in the early stage and gradually converge in the later stage. The acceptance criterion uses the Metropolis standard, which accepts inferior solutions with a certain probability when the energy function rises, to improve 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 grid model, and set appropriate grid cell size.

[0102] S3.2, construct a groundwater system numerical model that integrates structural priori, based on the spatial distribution characteristics of hydrogeological parameters, establish a complete hydrodynamic field boundary condition system, realize the dynamic coupling and parameterization of various boundary conditions.

[0103] After completing the basic grid construction and initial boundary condition setting, focus on the fine construction of the numerical model of the groundwater system in the modeling area, aiming to provide high-fidelity numerical simulation output that can be mapped to the graph structure for the graph neural network, including the following key links:

[0104] 1. Parameter space mapping to spatial domain:

[0105] Map the optimized multi-dimensional parameter sampling point set (including pumping volume Q, water level boundary R, permeability coefficient K, etc.) to the model spatial domain to form a parameter field consistent with the geological horizon structure. This mapping follows the geological structure control logic and strengthens the local parameter resolution in the boundary sensitive area, thereby ensuring the response accuracy of subsequent simulation in key areas.

[0106] 2. Dynamic boundary condition system construction:

[0107] Unlike traditional static boundary setting methods, introduce spatial position-dependent boundary parameter expression mechanisms (such as river boundary dynamic water level function, pumping well time-varying flow sequence), realize the dynamic coupling expression of spatio-temporal boundary conditions. This method is more in line with the needs of graph structure models to learn the change rules in the time dimension, and helps to improve the spatio-temporal prediction accuracy.

[0108] 3. Structure-consistent hydrodynamic field configuration:

[0109] According to the actual layered structure and aquifer connectivity of the groundwater system in the modeling area, the vertical seepage channels between multi-layer aquifer units are set, and the interlayer drainage coefficient and storage parameter are reasonably configured to reconstruct the complete hydrodynamic field structure in the area. This configuration not only improves the physical reality of numerical simulation, but also provides accurate structural priors for the cross-layer adjacency relationship in the graph neural network.

[0110] 4. Integration of boundary disturbance and response test mechanism:

[0111] The control disturbance mechanism (such as boundary pumping mutation, water head mutation) is introduced into the model to simulate the response behavior of the groundwater system under unsteady state conditions, and the disturbance response curve set is output. This response set is an important part of the graph neural network training data, which can be used to enhance the prediction robustness of the model under complex boundary scenarios.

[0112] 5. Consistency with the node structure of the graph neural network:

[0113] In the model output stage, the numerical simulation results (such as time series water level and pollutant concentration, hydraulic gradient, etc.) are reorganized according to the established grid node structure to ensure that their spatial distribution characteristics, adjacency relationships and node definitions in the graph neural network are consistent, realizing seamless docking from numerical model to graph structure model, providing structural complete supervision data support for subsequent graph neural network learning.

[0114] S4, cascaded graph neural network prediction.

[0115] S4.1, convert the grid data output by the groundwater system numerical model into graph structure data, with nodes representing hydrological units and edges representing spatial adjacency relationships, to ensure that the input data is highly compatible with the graph model design in terms of spatial structure, attribute dimension and time series organization, as follows:

[0116] 1. Node definition and attribute assignment:

[0117] The center of each hydrological unit in the three-dimensional grid model is taken as the node of the graph structure, and the node attributes include the water level and pollutant concentration of the corresponding hydrological unit at different times, hydraulic gradient, parameter value (such as permeability coefficient K, pumping quantity Q), etc., and auxiliary attributes such as layer number and burial depth can be added to provide multi-scale input features for the subsequent network coding module.

[0118] 2. Edge connection relationship construction:

[0119] Based on the adjacency information of the grid in the numerical model, the edge structure in the graph is constructed. The connection strategy of the edge includes both the connection between the horizontally adjacent units and the connectivity between the vertically adjacent layers, reflecting the hydraulic connection structure in the regional groundwater system. The attributes of the edge can include directionality, hydraulic gradient difference or boundary state identification, which can strengthen the representation ability of the graph structure to the real hydrodynamic process.

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

[0121] According to the time step of the model simulation, the graph structure at each time is organized into a graph sequence data, forming a graph time sequence (Spatio-temporal Graph Sequence). The graph at each time shares the topological structure, but the node attribute evolves and updates with time, meeting the requirements of spatio-temporal modeling in graph neural networks.

[0122] 4. Explicit encoding 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, which is combined with the node or edge attribute to form the input feature. This way can guide the graph neural network to explicitly pay attention to the difference of physical boundaries and improve the model fitting ability in the prediction of boundary areas.

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

[0125] The finally generated graph structure data is packaged into a standard graph data format (such as the graph data object supported by PyTorchGeometric or DGL), and normalized, missing data completed and abnormal data removed, etc. preprocessing operations are performed to provide high-quality, structurally consistent and time-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 follows:

[0127] Through the shared spatio-temporal attention mechanism, the collaborative modeling of the two tasks is realized.

[0128] 1. In the design of the water level prediction module, first, a multi-layer perceptron is used to encode the node features, and the input features include the three-dimensional spatial coordinates (x, y, z) of the hydrological unit, the key hydrogeological parameters (permeability coefficient K, specific yield Sy, etc.) and the historical water level observation sequence. The encoded features are processed through residual graph convolution blocks, each of which contains neighborhood information aggregation and jump connection, and its mathematical expression is where, is the hidden state of node i in the l-th layer network; W (l)Trainable parameter matrix of the first layer network; a ij is the normalized correlation strength between nodes i and j, and a is the Swish activation function. In particular, a boundary condition mask mechanism is designed for various types of hydrological boundary conditions (such as constant head boundary, pumping well boundary, etc.), which ensures that the prediction of the boundary nodes conforms to the physical law by imposing strong constraints.

[0129] The pollutant prediction module takes the water level prediction result as an important input feature, and splices it with the pollutant concentration feature to form an extended feature vector The module adopts a Gated Graph Convolutional Network (GRU) structure, and its state update formula is Through the gating mechanism, important neighborhood information is selectively retained, effectively capturing the migration and diffusion law of pollutants in the aquifer.

[0130] The two prediction modules share the core spatio-temporal attention mechanism. The spatial attention weight quantifies the hydraulic connection strength between hydrological units, and can automatically identify the influence range of pumping wells and other key spatial dependencies; the temporal attention weight models the evolution law of the state of the hydrological unit itself, and can accurately capture the seasonal fluctuations and other time series features. The entire network adopts LayerNorm for layer normalization, and introduces DropEdge regularization strategy to prevent overfitting, and gradient clipping to ensure numerical stability during training.

[0131] S4.3, design a loss function that integrates multi-task and physical constraints, including the relative Huber loss of water level prediction and the dimension preserving loss of pollutant prediction.

[0132] The joint optimization strategy of multi-task integration and physical constraints consists of three parts: the water level prediction loss term adopts an improved relative Huber loss function, defined as where δ is an adaptive threshold parameter, σ h is the standard deviation of water level observation, h i is the true water level value of the i-th node, is the predicted water level value, and N is the number of nodes; the pollutant prediction loss term adopts the mean square error that preserves the dimension of the physical quantity where c i is the true pollutant concentration, is the predicted concentration value. The physical constraint part contains three key contents: mass conservation constraint where is the gradient sign, p is the density, v is the Darcy velocity, and Q is the source and sink term, to ensure that the prediction result satisfies the continuity equation; boundary condition constraint where B is the set of boundary nodes, |B| is the number of boundary points, h bc is the known boundary value, the boundary node prediction value is forced to meet the known boundary condition; and the source-sink term constraint where n is the well wall normal vector, q i is the measured flux, ensuring that the predicted gradient near the pumping and injection wells is reasonable. The total loss function is obtained by Multi-objective optimization is achieved, where the weight coefficient is adjusted adaptively by 1, 2, 3, where σ k is the standard deviation parameter, ∈ is the smoothing coefficient, and is adjusted during the back propagation process by Dynamic balance is achieved, where is the partial derivative (i.e., the gradient) of the total loss function with respect to the dynamic weight parameter w k , w k is the dynamic weight parameter, and η is the learning rate. This system ensures that the prediction results not only fit the observed data, but also meet the basic principles of groundwater dynamics through strict mathematical constraints.

[0133] S4.4, to achieve effective simulation and prediction of groundwater dynamics, a joint optimization mechanism that combines data-driven and physical constraints is designed, i.e., a multi-task joint optimization strategy that integrates physical constraints, aiming to improve the prediction accuracy of the model while ensuring the reasonableness of the results in the physical sense.

[0134] The spatiotemporal collaborative state update and prediction mechanism includes the construction of a spatiotemporal attention gating mechanism to achieve dynamic evolution of node states and prediction output. Please refer to Figure 2 This embodiment uses a multi-head attention architecture to process spatiotemporal dependencies, where the spatial attention weight is calculated as:

[0135]

[0136] The temporal attention weight is calculated as:

[0137]

[0138] This mechanism can automatically identify key spatiotemporal patterns, where is a learnable parameter matrix, and d is the feature dimension.

[0139] The node state update uses a physically guided GRU unit:

[0140]

[0141] where the spatial neighborhood information is the spatial neighbor set of node i, a weight matrix of spatial attention values; temporal history information K w is the size of the temporal attention window, indicating the number of history time steps considered, is a weight matrix of temporal attention values. This design not only retains the advantages of traditional GRU processing time series data, but also enhances the extraction ability of spatial features through attention mechanism.

[0142] To ensure the physical rationality of the prediction results, the embodiment also imposes post-processing constraints on the output results, including forcing the predicted values of boundary nodes to comply with known boundary conditions.

[0143]

[0144] Embodiment two

[0145] The improved Latin hypercube sampling and adaptive simulated annealing algorithm are adopted to systematically construct 200 sets of training working condition data sets with strict physical constraints, providing high-quality input data for the spatio-temporal graph neural network. In the parameter space construction stage, the Stratified LHS method is adopted to divide each parameter dimension into 200 equal probability intervals and optimize the combination in the three-dimensional parameter space composed of pumping capacity Q (50-100 m 3 / d, logarithmic scale), water level boundary R (initial value 90.1-92.5 m, final value 81.25-82.5 m) and permeability coefficient K (0.0001-0.0007 m / s, stratified sampling). The initial parameter sample set is generated. Spatial correlation constraints are particularly introduced in the sampling process. Through the Maximin Criterion and parameter correlation optimization algorithm, the 200 initial parameter combinations generated satisfy the spatial coverage index h=0.18 (better than the traditional LHS of 0.25), the Pearson correlation coefficient between parameters |ρ|<0.1 and the boundary region sampling density is improved by 40%.

[0146] In the parameter optimization stage, a multi-objective simulated annealing algorithm based on physical constraints was designed. A composite objective function was constructed, including the sum of the reciprocal of the distance between sampling points (E_dist), the parameter correlation index (E_corr), and the boundary matching degree penalty term (E_bound), with adaptive weight coefficients of 0.5, 0.3, and 0.2, respectively. The optimization process adopted a temperature self-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 iteration number N = 1000. In the implementation, it was divided into three optimization stages: in the global exploration stage (1-300 iterations), a large neighborhood search with a 15% amplitude was implemented to optimize the global parameter distribution; in the local optimization stage (301-700 iterations), a medium perturbation with an 8% amplitude was used to adjust the well coordinates and boundary parameters, and an elite reservation mechanism was introduced; in the fine tuning stage (701-1000 iterations), a small perturbation with a 3% amplitude was used to optimize the boundary sensitive area, especially for the river boundary area (model 35 columns 10-55 rows), which was subjected to intensive optimization. The optimization results are shown in Figure 4 .

[0147] After the optimized parameter set passed the strict physical consistency test, the Flopy library was used to realize automatic modeling, and finally the MODFLOW6 model input parameters were integrated. The model system contains 5400 active units (75 rows × 125 columns grid), with the basic parameters of top elevation 100m, bottom elevation 70m, and initial water head 40m, and forms a stable water flow gradient through the northwest side (100m) and southeast side (80m) constant head boundary, as shown in Figure 3 . The specific parameter configuration includes: the spatial coordinates and dynamic pumping capacity parameters of the 10 optimized well sites in the Well package, the river boundary water level parameters using a linear gradient function in the RIV module, the permeability coefficient 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 -9 m / s) based on the permeability coefficient zoning in the RCH / EVT module, etc.

[0148] The 75 × 125 grid data output by MODFLOW6 is converted into a spatiotemporal graph structure, with each node containing spatial coordinates (x, y, z), permeability coefficient K, water level h, pollutant concentration c, and state characteristics of the previous five historical time steps (prev_h, prev_c). Two types of edge connections are constructed: spatial adjacency edges use a 4-neighborhood + vertical connection method, and time series edges associate the same node at consecutive time steps. Boundary conditions are explicitly coded by a 1-dimensional mask vector, including constant head boundary identification, river boundary water level value, and pumping well flow rate parameters, ensuring that the physical boundary information is completely preserved.

[0149] As shown in Figure 2 The network adopts a two-stage prediction architecture: the headwater prediction module (HeadGNN) contains 3 layers of residual graph convolution, each with a hidden dimension of 128, using a ReLU activation function; the contaminant prediction module (ConcGNN) adopts 2 layers of gated graph convolution (GRU unit), with a hidden dimension of 64, and the input is spliced with the headwater prediction results and the original node features (total dimension 192). Both modules output the final prediction value with a linear layer, and the hydrodynamic-contaminant transport collaborative modeling is realized through cascading connection.

[0150] The training adopts a composite loss function: the headwater prediction uses a relative Huber loss (δ = 0.1, σ_h = 1.2m), and the contaminant prediction adopts a dimension-preserving MSE loss (the concentration is normalized to [0, 1]). The physical constraint terms such as flux conservation and boundary matching are introduced, and λ = 0.3 is used to balance the two prediction tasks, and β = 0.5 is used to control the strength of the physical constraint, so that the prediction results meet the data fitting and physical laws at the same time.

[0151] The dynamic updating mechanism of 4-head spatial attention and 5-step time attention window is adopted. The Adam optimizer (learning rate 1e-4) is used for training, the batch size is 16, and the optimal result is achieved after 200 rounds of training (early stopping strategy). The results achieved by the final model on the training set, the validation set and the test set are shown in Table 1, wherein MSE_Head is the headwater mean square error, and MSE_Conc is the concentration mean square error.

[0152] Table 1 Prediction accuracy of cascaded graph neural network

[0153]

[0154] Therefore, the above-mentioned groundwater pollution concentration prediction method adopting the spatio-temporal graph neural network fused with physical constraints realizes the collaborative prediction of headwater and contaminant concentration, the physical constraint guarantees the prediction rationality, only a small amount of data is needed to establish the prediction model, the calculation burden of constructing a complex mechanism model in the traditional method is avoided, and the prediction accuracy under the condition of data scarcity is significantly improved.

[0155] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application but not to limit them, although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present application can still be modified or equivalently replaced, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present application.

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, For adaptive threshold parameters, The standard deviation of water level observations. For the first The actual water level value of each node To predict water level values, 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, 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, For density, For Darcy velocity, For source and sink items, For the set of boundary nodes, The number of boundary points, Given the boundary values, It is the normal vector. This is the measured flux; The total loss function is: ； In the formula, the weighting coefficients are adopted. Adaptive adjustment, The standard deviation parameter, This 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, To update the door, To reset the door, Candidate state Output the final state; 、 、 For a trainable parameter matrix, For nodes i Spatial neighbor set, This is the weight matrix for spatial attention values. For the size of the time attention window, This is the weight matrix for time attention values. Spatial attention weights, For time-based attention weights.

Images