Groundwater dynamic evolution prediction method

Abstract

The invention provides an underground water dynamic evolution prediction method, and belongs to the technical field of underground water prediction based on deep learning. Hydrological parameters and geographic position data of monitoring points are collected, and a regional hydrogeological input tensor is constructed; then establishing an initial finite element model to carry out space extrapolation and flow field completion on underground water level distribution; based on an initial finite element model structure, designing a global optimization algorithm based on target decomposition, carrying out real-time structure optimization on the finite element model, then simulating underground water evolution of a non-monitoring area through a finite element numerical solution, and generating multi-node physical consistency hydrological time series data after resolution is improved; and inputting the generated hydrological time series data into a designed perceptual water level prediction model, and outputting underground water level prediction results at multiple moments and multiple positions in the future and corresponding underground water level thermodynamic diagrams. According to the method, more accurate, continuous and interpretable high-precision prediction of the underground water system is realized, and reliable support is provided for scientific management of underground water resources.

Description

Technical Field

[0001] The present invention belongs to the technical field of groundwater prediction based on deep learning, and particularly relates to a method for predicting the dynamic evolution of groundwater. Background Art

[0002] As an important underground fresh water resource, the dynamic changes of groundwater are directly related to various social and economic activities such as agricultural irrigation, water resource scheduling, ecological protection, and urban water supply. Affected by multiple factors such as rainfall, evaporation, soil properties, geological structure, and human activities, the groundwater system usually has significant spatio-temporal heterogeneity and highly non-linear evolution characteristics. Currently, the commonly used groundwater prediction methods mostly rely on the time series analysis of monitoring points, lacking the modeling of the physical mechanism of groundwater and the full exploration of spatial flow field information, making it difficult to achieve reliable inference for non-monitored areas. At the same time, due to the limited density of sensor deployment and the incompleteness of data collection, the existing methods still have great limitations in data completion, boundary response, and evolution trend prediction, and it is difficult to meet the requirements of refined management and regional regulation.

[0003] Currently, the main methods for predicting the dynamic evolution of groundwater mainly include the following: Prediction methods based on statistical regression and time series: These methods are based on traditional mathematical statistics and time series modeling, including multiple linear regression (MLR), autoregressive moving average (ARMA / ARIMA), grey prediction model, etc. They mainly fit the time change trend through historical water level data to achieve short-term prediction of the groundwater level. Such methods have certain practicality in early groundwater research, with simple modeling and high computational efficiency. However, they generally rely on the stationarity of data and single-time dimension modeling, unable to effectively describe the non-linear evolution characteristics of groundwater in complex hydrogeological environments, and have limited response capabilities to sudden events and boundary condition changes, resulting in poor prediction accuracy and adaptability; Numerical simulation methods based on physical modeling: This method establishes the control equations of groundwater movement (such as Darcy's law, groundwater continuity equation, Richards equation, etc.), and combines regional boundary conditions, initial water level fields, and geological parameters, and uses numerical methods such as finite difference (FDM), finite element (FEM), or finite volume (FVM) to solve, so as to realize the spatial distribution simulation of the groundwater flow and evolution process. Its advantage lies in having an explicit physical mechanism interpretation ability and being suitable for studying groundwater systems with multiple boundaries and multi-layer structures. However, such methods are usually sensitive to the setting of the model structure, rely on a large number of measured parameters, and the cost of obtaining geological information is relatively high, making it difficult to be quickly deployed in practical projects. At the same time, its extrapolation ability in data missing areas is limited, and the simulation results are easily affected by initial and boundary errors; Data-driven machine learning prediction method: With the improvement of data acquisition and computing capabilities, machine learning and deep learning technologies have been gradually introduced into groundwater level prediction, such as support vector regression (SVR), random forest (RF), long short-term memory network (LSTM), gated recurrent unit (GRU), etc. These methods are good at dealing with non-linear sequence modeling problems and can directly learn the evolution laws from historical data without explicit physical equations. However, they generally have the following problems: insufficient ability to model spatial information, only able to predict single points or a few points; lack of physical constraints in the extrapolation area, and the interpretability of the results is poor; in addition, they have poor adaptability to boundary changes and extreme working conditions, the model stability is easily affected by the quality of training data, and it is difficult to be extended to highly heterogeneous formation areas.

[0004] In summary, the existing groundwater dynamic evolution prediction methods have obvious deficiencies in physical mechanism description, spatial completion ability and model stability, and it is difficult to simultaneously take into account the spatio-temporal evolution characteristics and engineering application requirements under complex geological structures. Summary of the Invention

[0005] In view of the above problems, the present invention aims to design a new groundwater prediction method that integrates physical constraint modeling and data-driven prediction, and has good spatial generalization ability and dynamic evolution expression ability, so as to achieve more accurate, continuous and interpretable high-precision prediction of the groundwater system, and thus provide reliable support for the scientific management and regulation of groundwater resources.

[0006] The present invention proposes a groundwater dynamic evolution prediction method, which includes the following steps: S1, collect hydrological parameter data, including groundwater level monitoring data, hydrogeological parameters and initial water level information, and collect the geographical location data of the monitoring points, and construct them into a regional hydrogeological input tensor; S2, based on the input data constructed in S1, establish an initial finite element model FEM to perform spatial extrapolation and flow field completion of the groundwater level distribution; S3, based on the initial finite element model FEM structure designed in S2, design a global optimization algorithm based on objective decomposition to perform real-time structural optimization on the regional division method, grid density distribution and boundary condition setting of the finite element model; S4, simulate the groundwater evolution in the non-monitoring area through the finite element numerical solution method to generate multi-node and high-resolution physically consistent hydrological time series data; S5, input the hydrological time series data generated in S4 into the constructed and trained flow field-aware water level prediction model, and output the groundwater level prediction results and their corresponding groundwater level heat maps at multiple future moments and multiple positions; the flow field-aware water level prediction model includes a spatio-temporal feature fusion encoder, a physically constrained groundwater level prediction module and a heat map generation module.

[0007] Preferably, the hydrological parameter data includes the groundwater level monitoring value collected by a pressure-type water level gauge , the soil moisture content collected by a soil moisture sensor , the pore water pressure collected by a pore water pressure sensor , the conductivity monitored by a conductivity sensor , the surface precipitation collected by deploying a rain gauge , the surface evaporation estimated by an evaporation pan combined with a temperature and humidity sensor , the temperature obtained by a meteorological sensor and humidity ; the above parameters together constitute the hydrological parameter data ; The geographical location data of the collection and monitoring points is selected in the area where the dynamic evolution of groundwater is to be carried out monitoring points, record the longitude and latitude position coordinates of the monitoring points, and obtain the geographical location sequence of the monitoring points .

[0008] Preferably, the specific content of S2 includes: S21, regional modeling and boundary initial setting: According to the obtained multi-monitoring point position coordinate data , establish a two-dimensional spatial network coordinate system including all monitoring points in the area , and obtain the initial hydrological data of all monitoring points at the initial time point based on the hydrological parameter data at the complete sampling time points, denoted as the initial groundwater level field , and use this initial groundwater level field as the starting state of the finite element simulation; in addition, set the boundary conditions as a mixed boundary including a constant boundary, a no-flow boundary, and a connected boundary; S22, grid structure design and parameter presetting: Initialization of the regional division method and optimization range: The parameter of the regional division method is defined as , and its value range includes: regular network division, Delaunay triangulation, quad-tree partitioning, and geological unit-guided partitioning, that is ; Initialization of the grid density parameter and optimization range: A total of sub-regions are obtained under the selected regional division method, and the grid density is set for different sub-regions , represents the grid density of the th region, and ; and set the highest grid density to , the lowest grid density is ; Therefore, the parameter to be optimized for the grid density is , and ; S23, Joint finite element modeling driven by multiple parameters: Construct a finite element model of multi-physical field coupling, based on the regional grid and the initial field as the input, and construct a system of multi-physical field coupling control equations according to the existing physical constraint equations, including: groundwater flow control equation, Richards equation, head-pressure conversion formula, conductivity migration equation, surface evaporation estimation equation, and heat conduction equation.

[0009] Preferably, the S3 establishes a multi-index optimization problem for the structure of the multi-physical field coupling finite element model based on the established decision variables and the established objective function; Among them, the decision variables are defined as: Based on the determined parameters to be optimized , , Determine the decision variables, the boundary parameter , where , 、 is the boundary flux adjustment parameter; the regional division method parameter , and ; The parameter to be optimized for the grid density is , and ; The complete set of decision variables is defined as: ; Among them, the objective function is designed as a multi-index weighted objective function:

[0010] Among them , , are hyperparameter weights; Completion accuracy loss , used to measure the predicted hydrological parameters of the monitored position points during the FEM completion process and the measured monitoring values ; Physical consistency error , used to measure the residual measure that violates the control equation in the solution obtained by FEM inference; Grid complexity penalty , used to limit the computational burden caused by excessive grid encryption.

[0011] Preferably, the specific process of the S4 includes: S41, Initial value assignment and boundary condition loading: According to the initial groundwater level field , and set the initial boundary parameter combination as , the initial regional division method is Regular network division, and the initial grid density are all preset as ; S42, Adopt an implicit finite element time-step advancing strategy, and use the coupled control equations to perform step-by-step evolutionary solution for each time step , and obtain the estimated values of all hydrological parameters in the entire regional coordinate system at the corresponding moment; S43, Based on the complete estimated results of the hydrological parameters in the obtained region, select the parameter supplementation degree as , and perform sampling in the regional coordinate system to obtain a total of groups of supplementary hydrological data , as well as the corresponding positions of the supplementary data ; Concatenate the supplemented data with the original collected data to form groups of complete hydrological data and ; And , , and form physically consistent water level time series data .

[0012] Preferably, the spatio-temporal feature fusion encoder in the flow field perception water level prediction model is used to extract enhanced spatio-temporal features, including a graph convolutional network layer, a dilated convolutional layer, a gated unit layer, and a spatio-temporal cross-attention mechanism layer; For the obtained physically consistent water level time series data , where is the time step, is the number of water level spatial position nodes, is the feature dimension, and the specific processing steps are as follows: S51, Input into the graph convolutional network layer to gradually extract spatial topological features. The graph convolutional network layer contains two graph convolutional layers, and each graph convolutional layer contains a graph convolutional kernel with 32 channels, a batch normalization layer, and a GeLU activation function; S52, Input into the dilated convolutional layer, capture long-period patterns through dilated convolution, learn the seasonal change trend of the water level, and obtain long-period time series features; S53, Input into the gated unit layer to learn the short-term fluctuation pattern of the data, learn the extreme change pattern of the water level, and obtain short-period time series features; S54. Input the spliced time series features obtained by channel splicing of the spatial topological features, long-term time series features, and short-term time series features into the spatio-temporal cross-attention mechanism to obtain spatio-temporal feature weights; S55. Multiply the spatio-temporal feature weights by the spliced time series feature points and then add them to the spatial topological features to obtain enhanced spatio-temporal features.

[0013] Preferably, the physical constraint groundwater level prediction module in the flow field perception water level prediction model is used to obtain the groundwater level prediction results at multiple future times; this module is based on the gated unit and the physical constraint of Darcy's law. The cell state update formula of the physical constraint groundwater level prediction module is specifically as follows:

[0014] Where is the cell state at is the output of the forget gate at is the output of the input gate at is the estimated cell state at is the Darcy's law constraint intensity coefficient, which controls the constraint intensity of Darcy's law on the gated unit's prediction of the groundwater level; is the discretized residual of Darcy's law, which reflects the degree to which the groundwater level predicted by the gated unit violates the physical law. The specific calculation formula is as follows:

[0015] Where, represents the gradient symbol, is the groundwater permeability, which is obtained by pumping tests; is the hydraulic gradient calculated based on the spatial position, which is calculated from the groundwater level monitoring values of adjacent spatial nodes in , is the source-sink term, which represents the net increase or decrease of water volume per unit area per unit time, and is the difference between the surface precipitation and surface evaporation in .

[0016] Preferably, the heat map generation module in the flow field perception water level prediction model outputs a groundwater level heat map; this module includes a grid conversion layer and a hierarchical upsampling decoder, specifically as follows: Input the groundwater level prediction results into the grid conversion layer, and use the inverse distance weighted algorithm to convert the groundwater level prediction results into a 64×64 low-resolution grid; The upsampling decoder for the low-resolution raster input level obtains the spatial distribution heat map of the groundwater level prediction result; the level upsampling decoder consists of four sequentially connected upsampling networks, which extract the features of the groundwater level prediction result layer by layer and improve the resolution to obtain the spatial distribution heat map of the 1024×1024 high-resolution water level prediction result; each upsampling network consists of a transposed convolutional layer, a sub-pixel convolutional layer, and a bilinear interpolation layer connected in sequence; the transposed convolutional layer is used to capture the features of the non-uniform groundwater distribution; the sub-pixel convolutional layer improves the resolution by channel reorganization of the features; the bilinear interpolation layer is used to ensure the smoothness of the groundwater level at different spatial positions and eliminate the stepped artifacts in the geological unit transition zone.

[0017] Preferably, the flow field-aware water level prediction model uses a model training strategy optimized by Bayesian optimization, specifically including: Construct a joint training parameter set for the spatio-temporal feature fusion encoder and the physically constrained groundwater level prediction module , where is the network structure parameter set of the spatio-temporal feature fusion encoder, is the network structure parameter set of the physically constrained groundwater level prediction module, is the learning rate, is the Darcy's law constraint strength coefficient; Use the Latin hypercube sampling method to generate M joint training parameter samples ; Construct a parameter design space optimization problem based on the KL divergence, specifically as follows:

[0018] where is the prior distribution of the groundwater level, , where is the mean of the groundwater level in is the standard deviation of the groundwater level in represents the Gaussian distribution; is the variational posterior distribution of the groundwater level, , where is the mean of the groundwater level prediction result obtained by predicting multiple future moments after configuring the spatio-temporal feature fusion encoder and the physically constrained groundwater level prediction module using , is the standard deviation of the groundwater level prediction result; represents the Gaussian distribution; For each sample in , solve the parameter design space optimization problem based on the KL divergence respectively to obtain the parameter set ; Taking the root mean square error between the predicted groundwater level result output by the physical constraint groundwater level prediction module and the true groundwater level result as the loss function, calculate the root mean square error obtained by separately configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module, and denote the parameter set that minimizes the root mean square error as the optimal model parameter set The root mean square error obtained by separately configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module is calculated, and the parameter set that minimizes the root mean square error is denoted as the optimal model parameter set 。

[0019] Compared with the prior art, the present invention has the following beneficial effects: (1) Flow field completion mechanism based on multi-physical field coupling: Integrate the multi-physical field control equations into the finite element simulation process, construct a groundwater system completion mechanism with strong physical consistency and tight spatio-temporal coupling, and significantly improve the integrity of groundwater level data and the credibility of modeling; (2) Structure-adaptive finite element modeling optimization strategy: A multi-objective joint optimization method integrating boundary adjustment, grid density control and regional division selection is proposed to realize the automatic adaptation and optimal configuration of the finite element modeling structure, and significantly improve the generalization ability and evolution simulation accuracy of the model under heterogeneous strata and complex boundary conditions; (3) Flow field-aware water level prediction model: It can fully capture the spatio-temporal evolution characteristics of hydrological data and generate a thermal map of the future multi-moment groundwater level spatial distribution that satisfies physical constraints. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 It is the overall technical route flow chart of the present invention.

[0021] Figure 2 It is the overall framework diagram of the physical modeling optimization module.

[0022] Figure 3 It is the structural schematic diagram of the flow field-aware water level prediction model.

[0023] Figure 4 It is the thermal map of groundwater dynamic prediction at different time steps in the embodiment.

[0024] Figure 5 It is the comparison chart of the groundwater level prediction accuracy of three models under the RMSE index in the embodiment.

[0025] Figure 6 It is the comparison chart of the groundwater level prediction accuracy of three models under the MAE index in the embodiment. DETAILED DESCRIPTION OF THE INVENTION

[0026] The present invention proposes a method for predicting the dynamic evolution of groundwater. The overall technical path includes the following core steps: First, a regional hydrogeological data set containing multi-source sensor data such as groundwater level, soil moisture, pore water pressure, conductivity, rainfall, evaporation, temperature and humidity is constructed, and based on the geographical location and time series structure, an input tensor with spatio-temporal constraints is formed; Subsequently, a finite element model containing multi-physics field control equations is established to perform spatial extrapolation and dynamic completion of the groundwater level distribution in the non-monitored area, enhancing the spatial integrity and physical consistency of the observation data; In addition, the present invention further designs an adjustable modeling structure, and through jointly optimizing the regional division method, grid density layout and boundary condition parameters, realizes the adaptive alignment of the finite element model structure with the actual hydrogeological evolution process; Finally, a deep learning model with flow direction perception ability is constructed, taking the completed high-quality time series data as input, jointly modeling the evolution trend of groundwater, realizing accurate prediction of multiple future moments and multiple spatial positions, and realizing a closed-loop modeling process from physical driving to intelligent prediction.

[0027] The overall process of this embodiment is as Figure 1 shown: Construction of hydrogeological data set: Collect hydrogeological parameter data, including groundwater level monitoring data, hydrogeological parameters and initial water level information, and collect geographical location data of monitoring points to construct a regional hydrogeological input tensor; Provide multi-dimensional input data for finite element simulation and structural optimization. At the same time, taking the future groundwater level observation results as labels, construct a prediction target set for supervised training; Design a finite element-driven flow field completion module: Based on the constructed input data, establish an initial finite element model (FEM) to perform spatial extrapolation and flow field completion of the groundwater level distribution; Provide more comprehensive data after completion for the prediction task; Physically modeled optimization module with adjustable structure: Based on the designed initial finite element simulation structure, design a global optimization algorithm based on objective decomposition to perform real-time structural optimization on the regional division method, grid density distribution and boundary condition settings of the finite element model; Make the finite element simulation results more consistent with the measured water level data, and improve the model accuracy and spatial adaptability; Generate multi-node, high-resolution physically consistent hydrogeological time series data: Simulate the evolution of groundwater in the non-monitored area through finite element numerical solutions to generate multi-node, high-resolution physically consistent hydrogeological time series data; Flow field perception water level prediction model: Input the generated hydrogeological time series data into the constructed and trained flow field perception water level prediction model, and output the predicted results of the groundwater level at multiple future moments and multiple positions and their corresponding groundwater level heat maps.

[0028] The following specifically describes the specific implementation process of the present invention in combination with specific embodiments.

[0029] I. Construction of Hydrogeological Dataset The present invention aims to realize the prediction of the dynamic evolution of groundwater based on the groundwater monitoring data within a region. Therefore, a hydrogeological dataset is first constructed. The specific content includes: Selection of hydrogeological parameters: According to the requirements of regional groundwater system modeling, key physical parameters affecting groundwater flow and evolution are selected, including: collecting groundwater level monitoring values using a pressure water level gauge , collecting soil moisture content using a soil moisture sensor , collecting pore water pressure using a pore water pressure sensor , monitoring electrical conductivity using an electrical conductivity sensor , deploying a rain gauge to collect surface precipitation , estimating surface evaporation using an evaporation pan combined with temperature and humidity sensors , obtaining temperature and humidity using a meteorological sensor; These parameters together constitute a hydrogeological parameter combination [gl, sc, wp, cd, sp, ei, tp,hm].

[0030] Multi-point monitoring data collection: Since the groundwater corresponding to different spatial positions is different and there is a direct influence between them, monitoring points are selected in the area where groundwater dynamic evolution is to be carried out, and multi-point monitoring data is collected. Specifically: (1) Record the longitude and latitude position coordinates of monitoring points to obtain a monitoring point geographical location sequence , that is , where represents the longitude and latitude position coordinates of the th monitoring point, and (2) At time point collect the hydrogeological parameter data of monitoring points , where represents the hydrogeological parameters collected at the th time point and the th monitoring point; (3) For monitoring points, collect the hydrogeological data at complete sampling time points , and ; Therefore, the overall input data includes the geographical location sequence of all monitoring points within the region and the hydrogeological parameters at complete sampling time points .

[0031] Groundwater level label: To achieve supervised modeling of the dynamic evolution process of groundwater, the present invention constructs a groundwater level label corresponding to the monitoring input as the target output for training the prediction model. Specifically: (1) Set the prediction time step to , and collect the groundwater level sequences of each monitoring point at the time point , where represents the groundwater level of the th monitoring point at the time point , and , ; (2) Collect the groundwater levels of the monitoring points over the complete prediction time step , and use them as the output data of the prediction model.

[0032] Hydrogeological dataset construction: The obtained input data and , as well as the obtained output data , are used as a set of hydrogeological datasets. Based on this method, a total of sets of datasets are collected, and finally the complete dataset collected and constructed is used as the hydrogeological dataset.

[0033] II. Design of a flow field completion model driven by the finite element method Due to the spatio-temporal heterogeneity of the groundwater system, it is difficult to comprehensively reflect the regional groundwater dynamic process only relying on the hydrogeological data collected at limited monitoring points. At the same time, the complexity of densely deploying sensors is high and the cost is expensive. Therefore, to solve the problems of incomplete data collection information and high cost caused by too many monitoring points, the present invention designs a flow field completion model driven by the finite element method. This model takes the multi-monitoring point data and constructed in stages as input, establishes a regional groundwater finite element model (FEM) with physical constraints, and conducts spatial deduction and dynamic completion of the water level state in the unobserved area under the control of boundary conditions and geological parameters, so as to obtain more complete hydrogeological information.

[0034] 1. Regional modeling and initial boundary setting: According to the multi-monitoring point position coordinate data obtained, establish a two-dimensional spatial network coordinate system that includes all monitoring points in the region, and obtain the initial hydrogeological data of all monitoring points at the initial time point according to the data, denoted as the initial groundwater level field , and use this initial groundwater level field as the starting state of the finite element simulation; In addition, to ensure the physical rationality of the simulation, the boundary conditions are set as a mixed boundary including a constant boundary, a no-flow boundary, and a communication boundary:

[0035] Among them, represents the hydrological data corresponding to the monitoring point with the coordinate point at the time point ; represents the transmission coefficient, is the outer normal vector of the boundary, and are the constant reference value and the far-field reference value of each hydrological parameter respectively; , 、 are the boundary flux adjustment parameters, and these three parameters together constitute the boundary parameters to be optimized , that is .

[0036] 2. Mesh structure design and parameter preset: After setting the regional boundary conditions, to support the stability of subsequent finite element calculations, spatial mesh division of the region is carried out and the structural parameters are initially set; specifically: (1) Initialization of the regional division method and optimization range: The parameter of the regional division method is defined as , and its value range includes: (regular network division), (Delaunay triangulation), (quad-tree partitioning), and (geological unit-guided partitioning), that is ; (2) Initialization of the mesh density parameter and optimization range: Under the selected regional division method, a total of sub-regions are obtained, and the mesh density is set for different sub-regions, represents the mesh density of the th region, and ; and the highest mesh density is set as , and the lowest mesh density is ; Therefore, the parameter of the mesh density to be optimized is , and .

[0037] 3. Joint finite element modeling driven by multiple parameters: Construct a finite element model of multi-physical field coupling, based on the regional mesh and the initial field Taking [input] as the input, a multi-physics coupling control equation system is constructed based on the existing physical constraint equations, including: groundwater flow control equation, Richards equation, head-pressure conversion formula, conductivity migration equation, surface evaporation estimation equation, and heat conduction equation; finally expressed as:

[0038] Among them, is the state vector, which contains all the hydrological parameters in ; represents the time partial derivative of all variables, represents the spatial gradient of all variables; is the two-dimensional spatial coordinate, is the time variable; is the total coupling control equation system containing the above physical constraint equations.

[0039] III. Construction of the physically-based modeling optimization module with adjustable structure Since the simulation accuracy of the finite element completion process highly depends on its modeling structure (including regional division, mesh density, and boundary condition settings), and it is difficult to simultaneously consider multi-region adaptability and the fitting effect of hydrological evolution characteristics by manual setting, the present invention designs a physically-based modeling optimization module with adjustable structure. This module automatically adjusts the finite element structure parameters by introducing a joint optimization strategy , , and , making the completion result closer to the measured data; the physically-based modeling optimization module is as shown in Figure 2 The specific process is as follows: 1. Definition of decision variables: Based on the determined parameters to be optimized , , the decision variables are determined. Specifically, the boundary parameter , where , 、 is the boundary flux adjustment parameter; the regional division method parameter , and ; the mesh density parameter to be optimized is , and ; Therefore, the complete decision variable set is defined as: .

[0040] 2. Design of the objective function: To measure the consistency between the finite element completion result and the measured data under the current structure setting, a multi-index weighted objective function is constructed, where , , is the hyperparameter weight, and the objective function specifically includes: (1) Completion accuracy loss , which is used to measure the predicted hydrological parameters at the monitored position points during the FEM completion process and the measured monitoring values The error between them is:

[0041] where and respectively represent the predicted and measured hydrological data at the time point and the th monitoring point; (2) Physical consistency error , which is used to measure the residual measure that violates the control equation in the solution obtained by FEM inference, that is:

[0042] where is the total coupled control equation system, is the regional spatial domain, is all time points, the set of predicted hydrological data at the monitoring points; (3) Mesh complexity penalty , which is used to limit the computational burden caused by excessive mesh refinement and is defined as:

[0043] Therefore, the optimization objective of the physically adjustable physical modeling optimization module is ; Based on the objective function established with the established decision variables, a multi-index optimization problem of the multi-physics field coupled finite element model structure is established.

[0044] 3. Design a global optimization algorithm based on objective decomposition to solve the multi-index optimization problem of the established physical field coupled finite element model structure, and obtain the optimal decision variable set , and the specific steps are as follows: 1) First, use the Latin hypercube sampling method to sample , , , and within their value ranges times to obtain Group decision variables; secondly, each group of decision variables is used as the parameters of the physical field coupling finite element model in turn, and the modified volume and fracture permeability corresponding to each group of decision variables are obtained by performing simulation calculations on the physical field coupling finite element model; thirdly, each group of decision variables and their corresponding completion accuracy loss, physical consistency error, and mesh complexity penalty are used to calculate the objective function value; finally, the decision variables and the objective values are combined to obtain the population , where the th solution , where is the objective value of the th solution; 2) Ascendingly sort the completion accuracy loss, physical consistency error, and mesh complexity penalty obtained for all solutions in . For the solutions ranked in the top in terms of completion accuracy loss, assign the skill factor to form the population . For the solutions ranked in the top in terms of physical consistency error, assign the skill factor to form the population . For the solutions ranked in the top in terms of mesh complexity penalty, assign the skill factor to form the population . The solutions not assigned the skill factor are formed into a diversity population ; ; 3) Initialize the number of evaluations, set the maximum number of evaluations, and set the cross-task evolution value ; 4) Generate a random number between 0 and 1 ; 5) When , perform crossover and mutation operations on , and respectively to obtain the offspring populations , and ; This way of evolving within the population enables the solutions of the three skill factors to explore the decision space with the three types of indicators of completion accuracy loss, physical consistency error, and mesh complexity penalty as the optimization objectives respectively, and obtain solutions that can achieve smaller values in terms of completion accuracy loss, physical consistency error, and mesh complexity penalty; When , randomly select the solutions of two skill factors to perform crossover and mutation operations to obtain the offspring populations , and ; This way of migrating and evolving between populations can promote the knowledge transfer between different skill factor populations and guide the population to converge to the global optimal objective value in the target space; ; 6) Update the evaluation count by adding the evaluation count to ; 7) Perform indicator-based natural selection on the merged populations of and with the merged population , and with the merged population as well as and with the merged population ; the smaller the indicator value, the higher the ranking of the solution. Among them, the complement precision loss indicator is used as the selection indicator, select the top population individuals to form the next generation population and update ; the physical consistency error indicator is used as the selection indicator, select the top population individuals to form the next generation population and update ; the grid complexity penalty indicator is used as the selection indicator, select the top population individuals to form the next generation population and update ; 8) Perform crossover and mutation operations on to obtain the offspring population ; 9) Update the evaluation count by adding the evaluation count to the size of ; 10) Perform natural selection based on the objective function value on the merged populations of and . The smaller the objective function value, the higher the ranking of the solution. Select the top solutions as the next generation population and update ; 11) Perform natural selection based on the objective function value on the merged populations of , , and . Select the top solutions as the next generation population and update ; 12) Judge whether the evaluation count is greater than the maximum evaluation count. If it is greater, output , otherwise repeat steps 2)-11); The output in step 12) is the optimal solution set of the multi-index optimization problem of the physical field coupling finite element model structure. Select the solution with the smallest objective function value in Decision variables As the optimal physical field coupling finite element model structure parameters; input the optimal physical field coupling finite element model structure parameters into the finite element model to obtain the optimal physical field coupling finite element model.

[0045] 4. Use the optimal physical field coupling finite element model to simulate the groundwater evolution in the non-monitored area and generate multi-node, high-resolution physically consistent hydrological time series data. Combine the time series simulation and the replenishment process: (1) Initial value assignment and boundary condition loading: According to the initial groundwater level field And set the initial boundary parameter combination as The initial regional division method Is (Regular network division), the initial grid density Are all preset as ; (2) Adopt the implicit finite element time step advancement strategy and use the coupled control equations For each time step Perform step-by-step evolution and solution to obtain the estimated values of all hydrological parameters in the entire regional coordinate system At the corresponding moment; (3) Complementary data output: Based on the complete estimated results of the hydrological parameters in the obtained area, select the parameter supplementation degree as And sample under the regional coordinate system To obtain a total of Groups of supplementary hydrological data And the corresponding positions of the supplementary data ; Jointly splice the supplemented data with the original collected data to form Groups of complete hydrological data And ; And ; ; Therefore, the supplemented data finally output by the finite element-driven flow field completion module is And And form physically consistent water level time series data .

[0046] IV. Design of the flow field sensing water level prediction model The flow field sensing water level prediction model designed by the present invention can intelligently generate high-precision groundwater level prediction results and heat maps at multiple future moments based on the FEM supplemented data after structure optimization; the model mainly consists of three parts: a spatio-temporal feature fusion encoder, a physically constrained groundwater level prediction module, and a heat map generation module, as specifically Figure 3 Shown.

[0047] 1. Spatiotemporal Feature Fusion Encoder The spatiotemporal feature fusion encoder extracts enhanced spatiotemporal features. The spatiotemporal feature fusion encoder includes a graph convolutional network layer, a dilated convolutional layer, a gated unit layer, and a spatiotemporal cross-attention mechanism layer; for the physically consistent water level time series data obtained from S4 , where is the time step, is the number of water level spatial position nodes, is the feature dimension. The specific steps are as follows: 1) Input into the graph convolutional network layer to gradually extract spatial topological features and explicitly model the spatial heterogeneity of the aquifer; the graph convolutional network layer contains two layers of graph convolutional layers, and each graph convolutional layer contains a graph convolutional kernel with 32 channels, a batch normalization layer, and a GeLU activation function; 2) Input into the dilated convolutional layer to capture long-term periodic patterns through dilated convolution, learn the seasonal change trend of the water level, and obtain long-term time series features; 3) Input into the gated unit layer to learn the short-term fluctuation pattern of the data, learn the extreme change pattern of the water level, and obtain short-term time series features; 4) Input the concatenated time series features obtained by channel concatenating the spatial topological features, long-term time series features, and short-term time series features into the spatiotemporal cross-attention mechanism to obtain spatiotemporal feature weights. The specific calculation formula is as follows:

[0048] where is the attention weight of the -th dimension of the spatial topological feature and the -th dimension of the concatenated time series feature, is the -th dimension of the spatial topological feature, is the -th dimension of the concatenated time series feature; 5) Multiply the spatiotemporal feature weights by the concatenated time series features and then add them to the spatial topological features to obtain enhanced spatiotemporal features.

[0049] 2. Physically Constrained Groundwater Level Prediction Module The physically constrained groundwater level prediction module is used to obtain the groundwater level prediction results at multiple future moments; this module is based on the gated unit and the physical constraint of Darcy's law. Different from the traditional gated unit, the cell state update formula of the physically constrained groundwater level prediction module is specifically as follows:

[0050] where is the -th moment cell state, is the output of the forget gate at time is the output of the input gate at time is the estimated cell state at time is the Darcy's law constraint strength coefficient, which controls the constraint strength of Darcy's law on the prediction of the groundwater level by the gated unit; is the discretization residual of Darcy's law, which reflects the degree to which the groundwater level predicted by the gated unit violates the physical law. The specific calculation formula is as follows:

[0051] where represents the gradient sign, is the groundwater permeability, which is obtained by pumping test; is the hydraulic gradient calculated based on the spatial position, which is calculated from the groundwater level monitoring values of adjacent spatial nodes in ; is the source-sink term, which represents the net increase or decrease of water volume per unit area per unit time, and is the difference between the surface precipitation and surface evaporation in .

[0052] 3. Heatmap Generation Module The heatmap generation module outputs the groundwater level heatmap; this module includes a raster conversion layer and a hierarchical upsampling decoder, as follows: 1) Input the groundwater level prediction result into the raster conversion layer, and use the inverse distance weighted algorithm to convert the groundwater level prediction result into a 64×64 low-resolution raster; 2) Input the low-resolution raster into the hierarchical upsampling decoder to obtain the spatial distribution heatmap of the groundwater level prediction result; the hierarchical upsampling decoder consists of four sequentially connected upsampling networks, which extract the features of the groundwater level prediction result layer by layer and increase the resolution to obtain a 1024×1024 high-resolution spatial distribution heatmap of the groundwater level prediction result; each upsampling network consists of a transposed convolutional layer (kernel size: 5×5, stride: 2, dilation rate: 2), a sub-pixel convolutional layer, and a bilinear interpolation layer connected in sequence; the transposed convolutional layer is used to capture the features of the non-uniform groundwater distribution and enhance the feature differences in the spatial position distribution of the groundwater level; the sub-pixel convolutional layer increases the resolution by re-organizing the channels of the features; the bilinear interpolation layer is used to ensure the smoothness of the groundwater level at different spatial positions and eliminate the stepped artifacts in the geological unit transition area.

[0053] 4. Model Training Strategy Based on Bayesian Optimization 1) Construct the joint training parameter set of the spatio-temporal feature fusion encoder and the physically constrained groundwater level prediction module , where is the set of network structure parameters of the spatio-temporal feature fusion encoder, is the set of network structure parameters of the physical constraint groundwater level prediction module, is the learning rate, is the Darcy's law constraint strength coefficient; 2) Use the Latin hypercube sampling method to generate M joint training parameter samples ; 3) Construct the parameter design space optimization problem based on KL divergence, specifically as follows:

[0054] where is the prior distribution of groundwater level, , where is the mean of the groundwater level in is the standard deviation of the groundwater level in represents the Gaussian distribution; is the variational posterior distribution of groundwater level, , where is the mean of the groundwater level prediction results obtained by predicting for multiple future moments after configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module using , is the standard deviation of the groundwater level prediction results; represents the Gaussian distribution; 4) For each sample in , solve the parameter design space optimization problem based on KL divergence respectively to obtain the parameter set ; 5) Use the root mean square error between the groundwater level prediction result output by the physical constraint groundwater level prediction module and the true groundwater level result as the loss function, calculate the root mean square error obtained by configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module using respectively, and denote the parameter set that minimizes the root mean square error as the optimal model parameter set .

[0055] V. Experimental Result Description To verify the effectiveness and advancement of the multi-physical-field-driven groundwater dynamic evolution prediction method (hereinafter referred to as "this method") proposed by the present invention, two types of comparative experiments are designed in this embodiment, and the method performance is evaluated from two aspects of the spatial visualization of the dynamic evolution result and the prediction accuracy. The comparative models include: (1) LSTM model: using the traditional long short-term memory network to perform time series prediction on the groundwater level; (2) GRU model: a time modeling method based on the gating mechanism is used as another comparative model.

[0056] 1. Analysis of thermal maps for groundwater dynamic distribution prediction This example selects 25 representative monitoring points within the region and, based on a high-dimensional hydrogeological dataset, constructs a multi-physics field coupled finite element prediction model. This model then performs spatiotemporal completion and evolution simulation of the groundwater flow field. This model then integrates a water level prediction module that incorporates flow direction sensing to output groundwater level distributions at multiple future moments. To ensure data security and model consistency, both the monitoring point locations and water level values are standardized. Figure 4 The heat map results of groundwater level prediction of this method at the 10th, 20th and 30th time steps are shown.

[0057] Results demonstrate that the proposed method demonstrates good continuity and physical rationality in spatial completion and time-series prediction, clearly restoring the primary groundwater flow direction and accurately depicting trends in areas of significant groundwater fluctuation (e.g., the northwest region). Compared to traditional methods, the proposed method maintains high distribution stability even in heterogeneous strata and complex boundary conditions, validating the effectiveness of its finite element structure construction and multi-physics field-driven modeling.

[0058] 2. Comparative experiment on groundwater level prediction accuracy To further quantify the prediction performance of the present invention, this example selects 10 monitoring points as target prediction points. The present method, LSTM, and GRU are used to predict groundwater levels at several future moments and compared with the measured values. The evaluation indicators used are: (1) Root Mean Square Error (RMSE), which reflects the error magnitude; and (2) Mean Absolute Error (MAE), which reflects the overall deviation. Figure 5 and Figure 6 The comparison results of RMSE and MAE indicators of the three models at 10 monitoring points are shown respectively. The horizontal axis is the monitoring point number and the vertical axis is the corresponding prediction error value.

[0059] The results show that the proposed method achieves the lowest prediction error across all monitoring points, with a stable RMSE of approximately 0.15 and a MAE between 0.12 and 0.13, significantly outperforming both the LSTM model (RMSE approximately 0.18, MAE approximately 0.15) and the GRU model (RMSE approximately 0.22, MAE approximately 0.19). Furthermore, the proposed method exhibits minimal error fluctuation across different monitoring points, demonstrating greater robustness and consistency across time and space.

[0060] Comprehensive analysis shows that the multi-physical-field coupling prediction framework proposed by the present invention not only maintains the model accuracy, but also has structural physical consistency and cross-regional generalization ability, and is especially suitable for groundwater evolution modeling tasks under heterogeneous strata and complex boundaries. Compared with traditional neural network methods, it has obvious advantages in error control, stability and practical deployment feasibility, providing effective technical support for high-precision groundwater prediction and regional water resources management.

[0061] The above are only the preferred embodiments of the present application and are not used to limit the present application. For those skilled in the art, various changes and modifications can be made to the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included within the protection scope of the present application.

[0062] Although the specific implementation manners of the present invention have been described above, they are not limitations on the protection scope of the present invention. Those skilled in the art should understand that various modifications or deformations that can be made without creative labor on the basis of the technical solution of the present invention are still within the protection scope of the present invention.

Claims

1. A method for predicting the dynamic evolution of groundwater, characterized in that, It includes the following steps: S1. Collect hydrological parameter data, including groundwater level monitoring data, hydrogeological parameters and initial water level information, and collect geographical location data of monitoring points to construct a regional hydrogeological input tensor; S2. Based on the input data constructed in S1, establish an initial finite element model FEM to perform spatial extrapolation and flow field completion of the groundwater level distribution; S3. Based on the initial finite element model FEM structure designed in S2, design a global optimization algorithm based on objective decomposition to perform real-time structural optimization on the regional division method, grid density distribution and boundary condition setting of the finite element model; S4. Simulate the evolution of groundwater in the non-monitoring area through the finite element numerical solution method to generate physically consistent hydrological time series data with multiple nodes and improved resolution; S5. Input the hydrological time series data generated in S4 into the constructed and trained flow field-aware water level prediction model to output the groundwater level prediction results and their corresponding groundwater level heat maps at multiple future times and multiple locations; the flow field-aware water level prediction model includes a spatio-temporal feature fusion encoder, a physically constrained groundwater level prediction module and a heat map generation module.

2. The groundwater dynamic evolution prediction method according to claim 1, characterized in that: The hydrological parameter data includes the groundwater level monitoring values collected by a pressure-type water level gauge , the soil moisture content collected by a soil moisture sensor , the pore water pressure collected by a pore water pressure sensor , the conductivity monitored by a conductivity sensor , the surface precipitation collected by deploying a rain gauge , the surface evaporation estimated by an evaporation pan combined with temperature and humidity sensors , the temperature obtained by a meteorological sensor and humidity ; The above parameters together constitute the hydrological parameter data ; The geographical location data of the collection monitoring points are selected in the area where groundwater dynamic evolution is to be carried out monitoring points, and record the longitude and latitude position coordinates of each monitoring point to obtain the geographical location sequence of the monitoring points .

3. The groundwater dynamic evolution prediction method according to claim 1, characterized in that: The specific content of S2 includes: S21, Regional Modeling and Initial Boundary Setting: Based on the acquired position coordinate data of multiple monitoring points , establish a two-dimensional spatial network coordinate system that includes all monitoring points within the region , and obtain the initial hydrological data of all monitoring points at the initial time point according to the hydrological parameter data at the complete sampling time points , which is denoted as the initial groundwater level field , and use this initial groundwater level field as the starting state for finite element simulation; in addition, set the boundary conditions as a mixed boundary including constant boundaries, no-flow boundaries, and connected boundaries; S22. Grid structure design and parameter presetting: Initialization of the regional division method and optimization scope: The parameter of the regional division method is defined as , and its value range includes: Regular network division, Delaunay triangulation, Quadtree partitioning, and Geological unit-guided partitioning, that is ; Grid density parameter initialization and optimization range: Under the selected area division method, a total of sub-regions are obtained, and grid densities are set for different sub-regions , denotes the grid density of the th region, and ; and the highest grid density is set to , and the lowest grid density is ; Therefore, the grid density parameter to be optimized is , and ; S23, Multi-parameter-driven combined finite element modeling: Construct a finite element model for multi-physical field coupling, with the regional grid as the basis and the initial field as the input, and construct a system of multi-physical field coupling control equations based on the existing physical constraint equations, including: groundwater flow control equation, Richards equation, head-pressure conversion formula, conductivity migration equation, surface evaporation estimation equation, and heat conduction equation.

4. The groundwater dynamic evolution prediction method according to claim 3, wherein: S3 establishes a multi-physical field coupling finite element model structure multi-index optimization problem based on the established decision variables and the established objective function; Among them, the decision variables are defined as: Based on the determined parameters to be optimized , , Determine the decision variables and boundary parameters , where , 、 is the boundary flux adjustment parameter; the regional division method parameter , and ; the parameter to be optimized for the grid density is , and ; the complete decision variable set is defined as: ; Among them, the objective function is designed as a multi-index weighted objective function: ； Among them , , are hyperparameter weights; Compensate for precision loss for measuring the predicted hydrological parameters at the monitored location points during the FEM compensation process and the measured monitoring values the error between them; Physical consistency error , which is used to measure the residual measure that violates the governing equation in the solution obtained by FEM inference; Mesh complexity penalty , which is used to limit the computational burden brought by excessive mesh refinement.

5. The groundwater dynamic evolution prediction method according to claim 3, characterized in that: The specific process of S4 includes: S41, Initial value assignment and boundary condition loading: According to the initial groundwater level field , and set the initial boundary parameter combination as , the initial regional division method is regular network division, and the initial grid density are all preset as ; S42. When adopting the implicit finite element time-step advancement strategy, use the coupled control equations For each time step Perform step-by-step evolution and solution to obtain the estimated values of all hydrological parameters in the entire regional coordinate system at the corresponding moment ; S43, based on the complete estimation results of the hydrological parameters in the obtained area, select the parameter supplementation degree as , and conduct sampling in the regional coordinate system to obtain a total of groups of supplementary hydrological data , as well as the corresponding positions of the supplementary data ; splice the supplemented data and the original collected data together to form groups of complete hydrological data and ; and , , and form physically consistent water level time series data .

6. The groundwater dynamic evolution prediction method according to claim 1, characterized in that: The spatio-temporal feature fusion encoder in the flow field perception water level prediction model is used to extract enhanced spatio-temporal features, including a graph convolutional network layer, a dilated convolutional layer, a gated unit layer, and a spatio-temporal cross-attention mechanism layer; for the obtained physically consistent water level time series data , where is the time step, is the number of water level spatial position nodes, is the feature dimension, and the specific processing steps are as follows: S51, input into the graph convolutional network layer to gradually extract spatial topological features. The graph convolutional network layer includes two graph convolutional layers, and each graph convolutional layer includes a graph convolutional kernel with 32 channels, a batch normalization layer, and a GeLU activation function; S52, input into the dilated convolutional layer, capture long-term patterns through dilated convolution, learn the seasonal variation trend of water level, and obtain long-term time series features; S53, input the short-term fluctuation pattern of the learning data of the input gating unit layer, learn the extreme change pattern of the water level, and obtain short-period time series features; S54. Input the concatenated time series features obtained by channel concatenation of spatial topological features, long-period time series features and short-period time series features into the spatio-temporal cross-attention mechanism to obtain spatio-temporal feature weights; S55. Multiply the spatio-temporal feature weights by the concatenated time series features and then add them to the spatial topological features to obtain enhanced spatio-temporal features.

7. The groundwater dynamic evolution prediction method according to claim 1, wherein: The physically constrained groundwater level prediction module in the flow field-aware water level prediction model is used to obtain the groundwater level prediction results at multiple future times; this module is implemented based on the gated unit and the physical constraint of Darcy's law. The cell state update formula of the physically constrained groundwater level prediction module is specifically as follows: ； where is the cell state at time is the output of the forget gate at time is the output of the input gate at time is the estimated cell state at time; is the Darcy's law constraint strength coefficient, which controls the constraint strength of Darcy's law on the gated unit for predicting the groundwater level; is the discretization residual of Darcy's law, which reflects the degree to which the groundwater level predicted by the gated unit violates the physical law. The specific calculation formula is as follows: ； Among them, represents the gradient symbol, is the groundwater permeability, obtained by pumping tests; is the hydraulic gradient calculated based on spatial positions, calculated from the groundwater level monitoring values of adjacent spatial nodes in is the source-sink term, representing the net increase or decrease in water volume per unit area per unit time, which is the difference between the surface precipitation and surface evaporation in 8. The groundwater dynamic evolution prediction method according to claim 1, characterized in that: The heat map generation module in the flow field-aware water level prediction model outputs the groundwater level heat map; this module includes a raster conversion layer and a hierarchical upsampling decoder, specifically as follows: Input the groundwater level prediction results into the raster conversion layer, and use the inverse distance weighted algorithm to convert the groundwater level prediction results into a 64×64 low-resolution raster; The upsampling decoder for the low-resolution raster input level obtains the spatial distribution heat map of the groundwater level prediction result; the upsampling decoder for the level consists of four successively connected upsampling networks, which extract the features of the groundwater level prediction result layer by layer and improve the resolution to obtain the spatial distribution heat map of the groundwater level prediction result with a high resolution of 1024×1024; each upsampling network is composed of a transposed convolutional layer, a sub-pixel convolutional layer, and a bilinear interpolation layer connected in sequence; the transposed convolutional layer is used to capture the features of the non-uniform groundwater distribution; the sub-pixel convolutional layer improves the resolution by channel recombination of the features; the bilinear interpolation layer is used to ensure the smoothness of the groundwater level at different spatial positions and eliminate the stepped artifacts in the geological unit transition zone.

9. A method for predicting the dynamic evolution of groundwater as described in claim 1, characterized in that: The flow field-aware water level prediction model uses a model training strategy optimized by Bayesian optimization, specifically including: Construct a combined training parameter set for a spatio-temporal feature fusion encoder and a physically constrained groundwater level prediction module , where is the network structure parameter set of the spatio-temporal feature fusion encoder, is the network structure parameter set of the physically constrained groundwater level prediction module, is the learning rate, is the Darcy's law constraint intensity coefficient; Generate M joint training parameter samples using the Latin hypercube sampling method ; Construct a parameter design space optimization problem based on the KL divergence, specifically as follows: ； Among them The prior distribution of the groundwater level , where is The mean value of the groundwater level in is The standard deviation of the groundwater level in Indicates a Gaussian distribution The variational posterior distribution of the groundwater level , where is the mean value of the groundwater level prediction results obtained by predicting multiple future moments after configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module using ; is the standard deviation of the groundwater level prediction results Indicates a Gaussian distribution For each sample in , solve the parameter design space optimization problem based on KL divergence respectively to obtain a set of parameters that satisfy Darcy's law ; Taking the root mean square error between the predicted groundwater level result output by the physical constraint groundwater level prediction module and the true groundwater level result as the loss function, calculate the root mean square error obtained by separately configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module, and record the parameter set that minimizes the root mean square error as the optimal parameter set of the model The root mean square error obtained by separately configuring the spatio-temporal feature fusion encoder and the physical constraint groundwater level prediction module. Denote the parameter set that minimizes the root mean square error as the optimal parameter set of the model .

Images