Groundwater level prediction method and system based on graph neural network and hydrogeological prior constraint

Abstract

The application provides a groundwater level prediction method and system based on a graph neural network and hydrogeological prior constraints, relates to the technical field of data processing, and comprises the following steps: obtaining groundwater level data and spatial coordinate data of each groundwater monitoring well in a target area; preprocessing the groundwater level data to obtain multi-well spatio-temporal sequence data; constructing a graph adjacency matrix between the groundwater monitoring wells according to the spatial coordinate data; inputting the constructed time sequence input sample into a spatio-temporal prediction network to obtain a water level prediction sequence; constructing a hydrogeological prior constraint residual according to the water level prediction sequence; performing constraint fine-tuning training on the spatio-temporal prediction network based on a joint loss function to obtain a trained spatio-temporal prediction network; and obtaining the groundwater level prediction results of each monitoring well according to the trained spatio-temporal prediction network, so that the multi-step prediction results considering data fitting and prediction sequence statistical consistency are formed under the condition of limited monitoring well data, and the prediction accuracy of the groundwater level is improved.

Description

Technical Field

[0001] This invention relates to data processing technology, and more particularly to a method and system for predicting groundwater levels based on graph neural networks and hydrogeological prior constraints. Background Technology

[0002] Groundwater is an important strategic resource for maintaining ecological balance, ensuring urban and rural water supply, and supporting economic and social development. Accurately predicting the dynamic changes in groundwater levels is of great significance for optimizing water resource allocation, preventing groundwater over-extraction, and protecting the ecological environment.

[0003] Currently, most existing groundwater level prediction methods employ data-driven models to learn statistical mapping relationships from historical water level sequences. However, they struggle to simultaneously consider both spatial correlations between wells and the temporal evolution of water levels. Spatially, existing methods typically model each monitoring well independently, failing to explicitly utilize the spatial hydraulic connections between wells and thus unable to perceive the impact of water level changes in one well on neighboring wells. Temporally, while existing methods can characterize the long-term dependence of water level evolution over time using recurrent neural networks, errors accumulate progressively with each prediction step, leading to a significant decrease in medium- to long-term prediction accuracy. Furthermore, when training samples are insufficient, monitoring wells are sparsely distributed, or external hydrological conditions change, the lack of explicit control over the relative relationships between wells and the smoothness of temporal evolution can result in predictions exhibiting asynchronous changes between adjacent wells, deviations from historical statistical relationships in relative water levels between wells, or short-term drastic fluctuations, affecting the stability of the prediction results.

[0004] Therefore, how to generate multi-step prediction results that take into account both data fitting and statistical consistency of prediction sequences under the condition of limited monitoring well data, and improve the accuracy of groundwater level prediction, has become an urgent problem to be solved. Summary of the Invention

[0005] This invention provides a groundwater level prediction method and system based on graph neural networks and hydrogeological prior constraints. It can generate multi-step prediction results that take into account both data fitting and statistical consistency of prediction sequences under the condition of limited monitoring well data, thereby improving the accuracy of groundwater level prediction.

[0006] A first aspect of the present invention provides a groundwater level prediction method based on graph neural networks and hydrogeological prior constraints, comprising:

[0007] The groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area are acquired, and the groundwater level data is preprocessed to obtain multi-well spatiotemporal sequence data.

[0008] A graph adjacency matrix is ​​constructed between each of the groundwater monitoring wells based on the spatial coordinate data, and a time-series input sample is constructed based on the multi-well spatiotemporal sequence data.

[0009] A spatiotemporal prediction network is constructed. The time-series input samples are input into the spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future. The spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer.

[0010] Hydrogeological prior constraint residuals are constructed based on the water level prediction sequence, and the hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term;

[0011] Based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, and the hydrogeological prior constraint residual, a joint loss function is constructed. The spatiotemporal prediction network is then fine-tuned and trained based on the joint loss function to obtain the trained spatiotemporal prediction network.

[0012] By inputting the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network, the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times are obtained.

[0013] Optionally, in one possible implementation of the first aspect, the preprocessing of the groundwater level data to obtain multi-well spatiotemporal sequence data includes:

[0014] Missing values ​​were removed and imputed from the groundwater level data of groundwater monitoring wells in various locations. The water level data of each monitoring well were converted into groundwater head values ​​or directly used as modeling data based on a unified elevation benchmark.

[0015] The water level data of each monitoring well are aggregated according to a preset time granularity to obtain the time series data of each monitoring well. The time series data of each monitoring well are then intersected and merged according to the monitoring time. The missing positions in the merged sequence are then interpolated to obtain a multi-well spatiotemporal sequence matrix.

[0016] Based on the water level data within the training time interval in the multi-well spatiotemporal sequence matrix, a normalization parameter is determined. The water level data in the multi-well spatiotemporal sequence matrix is ​​then mapped to a preset numerical range according to the normalization parameter to obtain multi-well spatiotemporal sequence data.

[0017] Optionally, in one possible implementation of the first aspect, constructing a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data includes:

[0018] The latitude and longitude coordinates of each groundwater monitoring well are projected onto a plane rectangular coordinate system, and the spatial distance between any two groundwater monitoring wells is calculated based on the projected plane rectangular coordinate system.

[0019] Based on the spatial distances described, the K-nearest neighbor method is used to determine the neighboring monitoring wells of each groundwater monitoring well, an initial adjacency matrix is ​​constructed, and the feature distance is determined based on the statistical characteristics of all spatial distances.

[0020] Based on the spatial distance and the feature distance, the adjacency weight between each groundwater monitoring well and its neighboring monitoring wells is determined. Each adjacency weight is filled into the corresponding adjacency position in the initial adjacency matrix. A self-loop is added to the initial adjacency matrix. The adjacency matrix after adding the self-loop is then symmetricized and normalized to obtain the graph adjacency matrix.

[0021] Optionally, in one possible implementation of the first aspect, constructing the time-series input sample based on the multi-well spatiotemporal sequence data includes:

[0022] From the multi-well spatiotemporal sequence data, extract the multi-well water level data corresponding to a first preset number of consecutive monitoring times as a sub-input sample, and determine the multi-well water level data corresponding to a second preset number of consecutive monitoring times after the first preset number of monitoring times as the multi-well water level prediction target corresponding to the sub-input sample.

[0023] According to the order of monitoring time, the original monitoring time intervals corresponding to the multi-well spatiotemporal sequence data are divided into training time intervals, verification time intervals, and testing time intervals.

[0024] Sub-input samples are constructed within the training time interval, verification time interval, and testing time interval respectively. Following the chronological order of the monitoring times of the multi-well spatiotemporal sequence data, one monitoring time is shifted backward each time. The operations of extracting the sub-input samples and determining the multi-well water level prediction targets are repeated to obtain multiple sets of sub-input samples and their corresponding multi-well water level prediction targets.

[0025] The time-series input sample is composed of multiple sets of corresponding sub-input samples and multiple well water level prediction targets.

[0026] Optionally, in one possible implementation of the first aspect, the construction of the spatiotemporal prediction network, which involves inputting the time-series input samples into the spatiotemporal prediction network to obtain a water level prediction sequence for each monitoring well at multiple consecutive monitoring times in the future, includes:

[0027] Based on the spatial feature extraction layer, the graph adjacency matrix, and the attention mechanism, the aggregation weight is determined, and the water level data of each monitoring well corresponding to each monitoring time are aggregated in the neighborhood to obtain multi-well spatial features.

[0028] Based on the time feature extraction layer, the spatial features of the multi-wells corresponding to each monitoring well are time-series modeled according to the order of monitoring time to obtain the time-series feature information of each monitoring well;

[0029] Based on the multi-step output layer and the temporal characteristic information of each monitoring well, the water level prediction sequence of each monitoring well is output for a second preset number of consecutive monitoring times in the future.

[0030] Optionally, in one possible implementation of the first aspect, the step of determining aggregation weights based on the spatial feature extraction layer, the graph adjacency matrix, and the attention mechanism, and performing neighborhood aggregation on the water level data of each monitoring well corresponding to each monitoring time to obtain multi-well spatial features, includes:

[0031] The water level data of each monitoring well at each monitoring time are sequentially linearly transformed, and the aggregation weight is determined based on the graph adjacency matrix and attention mechanism. Then, neighborhood feature aggregation, residual connection and nonlinear activation processing are performed to obtain the multi-well spatial features corresponding to each monitoring time.

[0032] The time feature extraction layer performs time-series modeling on the spatial features of the multi-wells corresponding to each monitoring well according to the chronological order of monitoring time, to obtain the time-series feature information of each monitoring well, including:

[0033] The spatial characteristics of each monitoring well at a first preset number of consecutive monitoring times are input into a loop network in chronological order of monitoring times. Based on the loop network, the temporal characteristic information of each monitoring well at the last historical monitoring time is output.

[0034] The step of outputting a water level prediction sequence for each monitoring well based on the multi-step output layer and the temporal characteristic information of each monitoring well for a second preset number of consecutive monitoring times includes:

[0035] Input the time-series characteristic information of each monitoring well at the last historical monitoring time into the linear mapping unit, and directly output the water level change of each monitoring well in the next second preset number of monitoring times;

[0036] Using the known water level corresponding to the last historical monitoring moment in the time-series input sample as a benchmark, the known water level is superimposed with the water level change to obtain the water level prediction sequence for each monitoring well in the next second preset number of monitoring moments.

[0037] Optionally, in one possible implementation of the first aspect, the construction of hydrogeological prior constraint residuals based on the water level prediction sequence, wherein the hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term, comprising:

[0038] The water level prediction sequence is concatenated with the water level data corresponding to the last historical monitoring moment in the time series input sample to obtain the continuous water level prediction sequence of each monitoring well, and the continuous water level prediction sequence is restored to the water level sequence with physical dimensions.

[0039] Based on the water level difference between adjacent monitoring times and the water level range of each monitoring well within the training time interval, the normalized water level change of each monitoring well is determined. Based on the deviation between the normalized water level change of each monitoring well at each monitoring time and the average normalized water level change of all monitoring wells at that monitoring time, and normalized by the overall change range of the normalized water level change of each monitoring well, the first constraint term is constructed.

[0040] Based on the difference in the average water level of each monitoring well within the training time interval, a benchmark for the relative difference between each pair of monitoring wells is determined. Based on the deviation between the water level difference of each pair of monitoring wells in the continuous water level prediction sequence and the benchmark for the relative difference between the pairs of monitoring wells, the second constraint term is constructed.

[0041] Based on the difference results of the values ​​of adjacent monitoring times of each monitoring well in the continuous water level prediction sequence, the third constraint term is constructed.

[0042] Based on the water level data within the training time interval, the reference scales corresponding to the first constraint term, the second constraint term, and the third constraint term are determined. The first constraint term, the second constraint term, and the third constraint term are divided by their respective reference scales. The first constraint term, the second constraint term, and the third constraint term after reference scale normalization are weighted and summed according to preset constraint weights to obtain the hydrogeological a priori constraint residuals.

[0043] Optionally, in one possible implementation of the first aspect, the step of constructing a joint loss function based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, and the hydrogeological prior constraint residuals, and performing constraint fine-tuning training on the spatiotemporal prediction network based on the joint loss function to obtain the trained spatiotemporal prediction network includes:

[0044] Based on the mean square error between the water level prediction sequence and the multi-well water level prediction target, the data prediction loss is obtained. The data prediction loss is then summed with the product of the preset prior constraint weights and the hydrogeological prior constraint residuals to obtain the joint loss function.

[0045] The spatiotemporal prediction network is iteratively trained based on the data prediction loss until the preset pre-training convergence condition is met, and the pre-trained network parameters are obtained.

[0046] The network parameters of the spatiotemporal prediction network are initialized to the pre-trained network parameters, and during the training iteration process, the prior constraint weights are gradually increased from the initial weight values ​​to the preset target weight values ​​according to the preset number of warm-up rounds.

[0047] The network parameters of the spatiotemporal prediction network are iteratively updated based on the joint loss function until the preset fine-tuning convergence condition is met, thus obtaining the trained spatiotemporal prediction network.

[0048] Optionally, in one possible implementation of the first aspect, after obtaining the predicted groundwater level for each monitoring well at multiple consecutive monitoring times, the method further includes:

[0049] The groundwater level prediction results of each monitoring well are converted into water level values ​​with physical dimensions, and the water level prediction values ​​of each monitoring well in each prediction time domain are determined according to multiple preset prediction time domains.

[0050] Based on the predicted water level and the actual water level of each monitoring well in each prediction time domain, the root mean square error and coefficient of determination of each monitoring well are calculated, and the root mean square error and coefficient of determination of all monitoring wells are taken as macro averages to obtain the macro average index corresponding to each prediction time domain.

[0051] The arithmetic mean of the macro-average indices corresponding to the third preset number of prediction time domains is used as the overall evaluation index for multi-step prediction, in order to evaluate the prediction accuracy of the spatiotemporal prediction network.

[0052] A second aspect of the present invention provides a groundwater level prediction system based on graph neural networks and hydrogeological prior constraints, comprising:

[0053] The processing module is used to acquire groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area, and to preprocess the groundwater level data to obtain multi-well spatiotemporal sequence data.

[0054] The construction module is used to construct a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data, and to construct a time-series input sample based on the multi-well spatiotemporal sequence data;

[0055] The prediction module is used to process the time-series input samples through a spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future; the spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer;

[0056] The constraint module is used to construct hydrogeological prior constraint residuals based on the water level prediction sequence. The hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term.

[0057] The training module is used to construct a joint loss function based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, as well as the hydrogeological prior constraint residual. First, the spatiotemporal prediction network is pre-trained based on the data prediction loss to obtain the pre-trained network parameters. Then, the spatiotemporal prediction network is fine-tuned based on the joint loss function using the pre-trained network parameters as initial parameters to obtain the trained spatiotemporal prediction network.

[0058] The output module is used to input the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network to obtain the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times.

[0059] A third aspect of the present invention provides a computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor, implements the groundwater level prediction method based on graph neural networks and hydrogeological prior constraints described in the first aspect.

[0060] The beneficial effects of this invention are as follows:

[0061] 1. This invention combines spatial modeling of graph structures, temporal modeling of cyclic networks, and statistical prior constraints determined by monitoring data during the training period. It adopts an incremental output method based on the last known water level. The spatial feature extraction layer adaptively determines the aggregation weight based on the characteristics of the monitoring wells through an attention mechanism within the adjacency range defined by the K-nearest neighbor graph. The temporal feature extraction layer performs temporal modeling of the spatial characteristics of each monitoring well. The multi-step output layer outputs the water level change of each monitoring well at multiple consecutive monitoring times in the future and superimposes it with the known water level at the last historical monitoring time to obtain the water level prediction sequence. This anchors the prediction results to the latest observation, thereby simultaneously utilizing the spatial correlation information between wells and the temporal evolution information of water level under sparse monitoring well conditions, improving the accuracy and stability of medium- and long-term water level prediction.

[0062] 2. This invention uses groundwater head as the modeling variable when a unified elevation benchmark is available. Based on monitoring data within the training time interval, it constructs a water level change consistency constraint, an inter-well head difference stability constraint, and a temporal smoothing constraint. These three constraints are normalized using a reference scale calculated within the training time interval. Statistical prior regularization is applied to the prediction sequence from three perspectives: consistency of changes across multiple wells, relative relationships between wells, and smoothness of temporal evolution. Thus, when training data is insufficient or hydrological conditions change, the constrained prediction results do not exhibit unreasonable phenomena such as drastically opposite water level changes among wells, significant deviations in inter-well water level differences from historical levels, or repeated oscillations in water level evolution trends. This improves the hydrogeological rationality and credibility of the prediction results.

[0063] 3. This invention adopts a two-stage training method of pre-training followed by constraint fine-tuning, and sets up a prior constraint weight warm-up mechanism. In the pre-training stage, the network is trained only based on the data prediction loss, so that the network can fully learn the water level change pattern in the data. In the fine-tuning stage, the pre-trained network parameters are used as the initial parameters, and the prior constraint weights are gradually increased from the initial weight values ​​to the preset values ​​according to the preset warm-up rounds. The optimal network parameters are determined based on the joint loss of the validation set, so that the network has both good prediction accuracy and hydrogeological consistency. Attached Figure Description

[0064] Figure 1 This is a flowchart illustrating a groundwater level prediction method based on graph neural networks and hydrogeological prior constraints provided in an embodiment of the present invention.

[0065] Figure 2 This is a schematic diagram comparing the prediction accuracy of each prediction model in different prediction time domains in the embodiments of the present invention;

[0066] Figure 3 This is a schematic diagram illustrating the changes in prediction accuracy of each prediction model under different training data conditions in this embodiment of the invention.

[0067] Figure 4 This is a schematic diagram of a groundwater level prediction system based on graph neural networks and hydrogeological prior constraints provided in an embodiment of the present invention. Detailed Implementation

[0068] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0069] See Figure 1 This is a flowchart illustrating a groundwater level prediction method based on graph neural networks and hydrogeological prior constraints provided in an embodiment of the present invention. Figure 1 The method shown can be implemented by software and / or hardware devices. It includes steps S1 to S6, as detailed below:

[0070] S1. Obtain groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area, and preprocess the groundwater level data to obtain multi-well spatiotemporal sequence data.

[0071] The target area refers to the geographical range within which groundwater level prediction is required. Several groundwater monitoring wells are deployed within this area. Groundwater monitoring wells are artificial or natural wells used to observe dynamic changes in groundwater levels. Each monitoring well has a unique number and a fixed geographical location. Groundwater level data refers to the groundwater head values ​​recorded at various groundwater monitoring wells at multiple consecutive monitoring times. Spatial coordinate data refers to the geographic spatial location information of various groundwater monitoring wells. Multi-well spatiotemporal sequence data refers to the multidimensional data structure obtained by arranging the groundwater level data according to the monitoring time and monitoring well number after preprocessing.

[0072] It is understandable that raw water level data obtained directly from the monitoring site may contain various quality issues, such as erroneous readings caused by monitoring equipment failure, data loss due to communication interruption, and numerical deviations caused by sensor drift. If these unprocessed data are used directly, the errors and noise in the data will interfere with the judgment of the actual water level change pattern. Therefore, the raw data obtained can be standardized first to ensure that the data from all monitoring wells meet the same standard in terms of quality and scale.

[0073] Specifically, groundwater level data and corresponding spatial coordinate data of each well can be obtained from historical monitoring records of groundwater monitoring wells in the target area at multiple consecutive monitoring times. The obtained groundwater level data is then preprocessed to eliminate quality defects and numerical scale differences between different monitoring wells, so that the processed data meets the input requirements of subsequent processing steps in terms of numerical quality and range. After preprocessing, a multi-well spatiotemporal sequence data with a uniform format and orderly arrangement can be obtained.

[0074] In some embodiments, step S1, "preprocessing the groundwater level data to obtain multi-well spatiotemporal sequence data," includes the following steps:

[0075] S11 performs missing value removal and interpolation processing on the groundwater level data of groundwater monitoring wells in various locations, and converts the water level data of each monitoring well into groundwater head values ​​or directly uses them as modeling data based on a unified elevation benchmark.

[0076] Among them, missing value removal and interpolation processing refers to the process of filling in the gaps in the data when water level data is not recorded at certain monitoring times due to monitoring equipment failure, communication interruption, etc., by using existing normal data before and after the monitoring time and estimating reasonable replacement values ​​through reasonable mathematical methods. For abnormal values ​​that deviate significantly from the reasonable range, they are identified and removed, and the removed positions are filled using neighboring data. Unified elevation benchmark refers to a unified reference elevation surface agreed upon for water level measurements of each monitoring well. The measurement results of each monitoring well are converted to this reference surface, making the water level data of different monitoring wells spatially comparable. Groundwater head value refers to the groundwater pressure head height measured with a unified elevation benchmark surface, determined by the difference between the elevation of the monitoring well measurement point and the water level burial depth, with the unit being meters. When a unified elevation benchmark cannot be obtained, the water level burial depth value is directly used as the modeling data, with the unit being meters. Modeling data refers to the data actually used to construct the multi-well spatiotemporal sequence matrix and subsequent model training and prediction after the above processing.

[0077] It is understandable that water level depth reflects the depth of the water surface relative to the surface benchmark. The larger the depth value, the lower the water level. Hydraulic head value, on the other hand, represents the height of the water level relative to a unified elevation benchmark. The two values ​​have different meanings. When a unified elevation benchmark for all monitoring wells cannot be obtained, the difference in depth values ​​between different monitoring wells cannot be interpreted as a difference in hydraulic head, hydraulic gradient, or steady-state water level difference. However, the depth sequence of each monitoring well can still accurately reflect the change of its water level over time. Therefore, water level depth can be directly used as a modeling and prediction variable. When a unified elevation benchmark is obtained, the hydraulic head value of each monitoring well can be determined according to the difference between the elevation of the measurement point and the water level depth for cross-well comparison.

[0078] Specifically, the original water level records of each monitoring well are traversed. For missing data locations, the missing values ​​are estimated based on the existing normal data before and after the location. For outliers that deviate significantly from the reasonable range, they are removed and filled in the same way. After the missing values ​​are removed and interpolated, it is determined whether a unified elevation datum can be obtained. If it can be obtained, the water level data of each monitoring well is uniformly converted into the head value measured at the elevation datum. If it cannot be obtained, the water level depth value is directly used as the modeling data.

[0079] S12, aggregate the water level data of each monitoring well according to the preset time granularity to obtain the time series data of each monitoring well, and then take the intersection and union of the time series data of each monitoring well according to the monitoring time. The missing positions in the merged sequence are interpolated to obtain the multi-well spatiotemporal sequence matrix.

[0080] Among them, the preset time granularity refers to a pre-set uniform time unit, such as a natural day; the time series data refers to the water level series data obtained by averaging the water level values ​​of each monitoring well on a daily basis; the monitoring time refers to the time corresponding to each time unit after aggregation, such as when the granularity is a natural day, each natural day is a monitoring time; and the multi-well spatiotemporal sequence matrix refers to the matrix structure formed by arranging the water level values ​​of each monitoring well at the same monitoring time according to the two dimensions of monitoring time and monitoring well number.

[0081] It is understandable that the original observation frequency and data recording period of different monitoring wells may differ. For example, some monitoring wells may record once a day, while others may record multiple times a day, and the start and end times of recording for each monitoring well may be different. If the original records of each monitoring well are directly spliced ​​together, the monitoring wells may not be aligned on the time axis, and a regular spatiotemporal matrix cannot be formed. Therefore, the water level data of each monitoring well can be aggregated to the same time granularity first, and then the intersection can be taken according to the monitoring time. Only the time period in which each monitoring well has a record can be retained, and the few missing positions that still exist after the intersection can be interpolated to obtain a time-aligned, complete and continuous multi-well spatiotemporal sequence matrix.

[0082] Specifically, the water level values ​​of each monitoring well can be averaged daily using natural days as the time granularity to obtain the daily water level sequence of that well. The daily water level sequences of each monitoring well are then intersected and merged according to date, retaining only the dates on which all monitoring wells have observation data. For individual dates that still have gaps in the merged sequence, normal data from adjacent dates are used for interpolation to obtain a complete multi-well spatiotemporal sequence matrix. Each row of this matrix corresponds to a monitoring time, and each column corresponds to a monitoring well.

[0083] For example, if there are 3 monitoring wells in the target area, and a total of 100 complete records of monitoring times are kept, then a multi-well spatiotemporal sequence matrix of 100 rows × 3 columns is constructed. Each row corresponds to a monitoring time, and each column corresponds to a monitoring well. Each element in the matrix stores the water level value of the monitoring well at that monitoring time. After the matrix is ​​constructed, it can be used for subsequent steps.

[0084] S13, determine the normalization parameter based on the water level data in the training time interval of the multi-well spatiotemporal sequence matrix, and map the water level data in the multi-well spatiotemporal sequence matrix to a preset value interval according to the normalization parameter to obtain multi-well spatiotemporal sequence data.

[0085] Among them, the training time interval refers to the time range defined from the multi-well spatiotemporal sequence matrix and used only for model training; the normalization parameter refers to the minimum and maximum values ​​used to map the original values ​​to the preset value interval proportionally; the preset value interval refers to the target interval range for normalization processing, such as [-1, 1]; and the normalization processing refers to scaling the water level data according to a uniform ratio so that the data falls within the preset value interval.

[0086] It is understandable that different monitoring wells are located in different geographical locations and have different hydrogeological conditions, so their water level values ​​may vary. If well data of different magnitudes are directly mixed and used, the monitoring well data with larger values ​​will dominate in subsequent processing, while the water level change information contained in the monitoring well data with smaller values ​​will be weakened due to the difference in magnitude. In addition, if data from all time periods is used to determine the normalization parameter, statistical information from future time periods will be leaked to the training process, affecting the objectivity of model evaluation. Therefore, the normalization parameter can be determined based only on the data within the training time interval, and this normalization parameter can be used to uniformly map the data from all time periods.

[0087] Specifically, data corresponding to the training time interval can be extracted from the multi-well spatiotemporal sequence matrix. The minimum and maximum values ​​of the water level in each monitoring well are statistically analyzed and used as normalization parameters for each monitoring well. Based on these normalization parameters, the water level values ​​of each monitoring well in the multi-well spatiotemporal sequence matrix are mapped proportionally to a preset numerical range. For example, assuming there are two monitoring wells in the target area, and the water level of well A fluctuates between 200 meters and 220 meters, then the minimum value of 200 meters and the maximum value of 220 meters for well A, and the minimum value of 10 meters and the maximum value of 15 meters for well B, are used as the normalization parameters for each monitoring well. When the preset numerical range is [-1, 1], each monitoring well is mapped according to its own minimum and maximum values. For well A, 200 meters and 220 meters correspond to -1 and 1, respectively, and for well B, 10 meters and 15 meters correspond to -1 and 1, respectively. After this processing, the water level data of different monitoring wells are mapped to the same numerical range, resulting in normalized multi-well spatiotemporal sequence data.

[0088] S2, construct a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data, and construct a time-series input sample based on the multi-well spatiotemporal sequence data.

[0089] The graph adjacency matrix is ​​a data structure stored in the form of a mathematical matrix to record whether there is a spatial correlation between each groundwater monitoring well and other groundwater monitoring wells and the strength of the correlation. Each element in the matrix represents the connection weight between the corresponding two groundwater monitoring wells. The larger the weight, the tighter the spatial correlation. The time series input sample refers to the data unit extracted from the multi-well spatiotemporal sequence data in chronological order. Each data unit contains water level data from multiple consecutive historical monitoring times, and the water level data from the next monitoring time as the corresponding prediction target.

[0090] It is understandable that groundwater is not static but flows slowly within aquifers. Therefore, monitoring wells in different locations are not independent of each other. The water level change of one monitoring well is often related to the water level changes of surrounding monitoring wells. Generally, the closer two monitoring wells are, the shorter the path of groundwater flow between them, and the stronger the mutual influence of water level changes. If only each monitoring well is analyzed individually, it is equivalent to ignoring the objective hydraulic connection between the monitoring wells, which may result in the loss of important spatial information. Therefore, we can first use the spatial coordinates of each monitoring well to calculate the distance between wells and construct a graph adjacency matrix that can quantify the spatial correlation strength between monitoring wells. This expresses the spatial hydraulic connection of groundwater in numerical form. At the same time, the spatiotemporal sequence data of multiple wells can be divided into several sets of input and prediction target paired samples to provide structured data input for subsequent steps.

[0091] Specifically, a graph adjacency matrix can be constructed between monitoring wells based on the spatial coordinate data of each monitoring well. This matrix records whether there is a spatial correlation between any two monitoring wells and the strength of the correlation. At the same time, several sets of data units can be extracted from the spatiotemporal sequence data of multiple wells in chronological order. Each set of data units contains water level data of multiple wells at a continuous historical monitoring time, as well as the target water level to be predicted after this historical data. These data units together constitute the time series input sample.

[0092] In some embodiments, step S2, "constructing a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data," includes the following steps:

[0093] S21, project the latitude and longitude coordinates of each groundwater monitoring well to a plane rectangular coordinate system, and calculate the spatial distance between any two groundwater monitoring wells based on the projected plane rectangular coordinate system.

[0094] Among them, latitude and longitude coordinates refer to the geographical location information of the monitoring well on the Earth's surface expressed in longitude and latitude, plane rectangular coordinates refer to the coordinates obtained after converting the latitude and longitude coordinates to a plane coordinate system with meters as the unit according to the equidistant projection method, and spatial distance refers to the straight-line distance between two monitoring wells on the plane.

[0095] It is understandable that latitude and longitude coordinates are angular quantities. Directly calculating the distance between wells using the difference between latitude and longitude will result in scale distortion due to different latitudes. The actual ground distance corresponding to the same difference in longitude at different latitudes is different. Therefore, we can first use the average value of the latitude and longitude of each monitoring well in the target area as the projection reference, project the latitude and longitude coordinates into plane rectangular coordinates, and then calculate the distance between wells in the plane rectangular coordinate system to obtain a distance value that accurately reflects the actual spatial distance between wells.

[0096] Specifically, the average latitude and longitude of all monitoring wells in the target area can be calculated. Using this average value as the projection reference, the latitude and longitude coordinates of each monitoring well can be converted into plane rectangular coordinates in meters. Based on the plane rectangular coordinates of each monitoring well, the Euclidean distance between any two monitoring wells can be calculated to obtain the spatial distance between the two groundwater monitoring wells.

[0097] S22, based on the spatial distances, the K-nearest neighbor method is used to determine the neighboring monitoring wells of each groundwater monitoring well, an initial adjacency matrix is ​​constructed, and the feature distance is determined based on the statistical characteristics of all spatial distances.

[0098] The K-nearest neighbor method refers to selecting the K nearest monitoring wells from the remaining monitoring wells as its neighboring monitoring wells for each monitoring well. K is a preset positive integer, such as K=2. The neighboring monitoring wells refer to the multiple monitoring wells that are spatially closest to the current monitoring well selected by the K-nearest neighbor method. The initial adjacency matrix is ​​a matrix structure initially established based on the determined adjacency relationships between each monitoring well. This matrix only retains monitoring well pairs that have mutual adjacency relationships. The feature distance refers to the feature value used to represent the spatial scale difference between monitoring wells in the target area. For example, it can be the median of all spatial distances.

[0099] Understandably, establishing connections between every monitoring well and all other monitoring wells in the area would be tantamount to assuming a direct hydraulic link between any two wells, which is unrealistic. Between two monitoring wells located far apart, groundwater requires a long flow path to influence each other, resulting in negligible correlation. Forcing connections would introduce irrelevant interference. Therefore, the K-nearest neighbor method can be used to filter connections, allowing each monitoring well to connect only with its closest few wells. An initial adjacency matrix is ​​then constructed based on these established adjacency relationships. This matrix retains only connections between wells with genuine spatial proximity, avoiding interference from irrelevant long-distance connections.

[0100] Specifically, based on the calculated spatial distance between any two monitoring wells, each monitoring well is sorted in ascending order of its distance from other monitoring wells. The K nearest monitoring wells are selected as its neighboring monitoring wells. For example, if there are 3 monitoring wells in the target area, and K=2, then for well 1, its corresponding neighboring monitoring wells are well 2 and well 3. Thus, the neighboring relationship of well 1 is determined as well 2 and well 3. This filtering operation is performed on all monitoring wells one by one to determine the neighboring monitoring wells of each monitoring well. Then, based on the determined neighboring relationships between the monitoring wells, an initial adjacency matrix is ​​constructed. In the matrix, only the positions corresponding to monitoring well pairs with neighboring relationships are connected, and no connections are established for the remaining positions.

[0101] In addition, the well distance can also be used to determine the K nearest neighbor range of each monitoring well, and then determine which positions in the adjacency matrix can participate in the subsequent neighborhood aggregation calculation. In the subsequent graph attention layer, the adjacency matrix limits the adjacent wells that can participate in aggregation through the positions of non-zero elements. The numerical distance weight is not directly added to the attention score. The actual aggregation weight within the neighborhood is determined by the monitoring well features through a learnable attention mechanism. When there are only 3 monitoring wells and K=2 in the embodiment, each well is connected to the other two wells, and the adjacency range is fully connected. Under this data condition, the well distance does not produce a further neighbor filtering effect.

[0102] Specifically, for each monitoring well, the K monitoring wells with the smallest distance between them are selected as its neighboring monitoring wells. The median of all spatial distances between the wells is calculated as the feature distance. Then, based on the determined proximity relationships between the monitoring wells, an initial adjacency matrix is ​​constructed. In the matrix, only the locations corresponding to the monitoring well pairs with proximity relationships are connected, and no connections are established for the remaining locations.

[0103] For example, if there are 3 monitoring wells in the target area, when K is set to 2, each monitoring well is connected to the other two monitoring wells. After the adjacency matrix constructed from the distance between the wells is symmetricized, self-loops are added and normalized, the positions of its non-zero elements are used to determine the range of aggregateable nodes of the graph attention layer. The attention layer then calculates and normalizes the aggregation weight between the wells based on the characteristics of each monitoring well at the current time.

[0104] S23, based on the spatial distance and the feature distance, determine the adjacency weight between each groundwater monitoring well and its neighboring monitoring wells, fill each adjacency weight into the corresponding adjacency position in the initial adjacency matrix, add a self-loop to the initial adjacency matrix, and perform symmetry and normalization processing on the adjacency matrix after adding the self-loop to obtain the graph adjacency matrix.

[0105] Among them, adjacency weight refers to the numerical value that quantifies the spatial correlation strength between two monitoring wells, adjacency position refers to the position of the matrix element corresponding to the monitoring well pair with a neighboring relationship in the initial adjacency matrix, self-loop refers to the connection between the monitoring well and itself, and normalization processing refers to scaling the rows or columns of the adjacency matrix to keep the feature scale of each monitoring well consistent during aggregation calculation.

[0106] It is understandable that the strength of hydraulic connection between wells usually weakens with increasing distance. Therefore, the adjacency weight can be set as a decreasing function that decreases with increasing spatial distance. In order to eliminate the influence of spatial scale differences in different regions, the ratio of distance to feature distance can be used as the independent variable for weight calculation. For example, the adjacency weight can be set as the reciprocal of this ratio plus one. Adding a self-loop can allow the monitoring well to retain its own water level information when the neighborhood is aggregated, avoiding the complete submersion of its own features during the aggregation process. Symmetry and normalization processing can ensure the symmetry and numerical stability of the adjacency matrix, which is convenient for subsequent matrix operations with the feature matrix.

[0107] Specifically, for two adjacent monitoring wells, the ratio of their spatial distance to their characteristic distance is calculated. The ratio is incremented by one and the reciprocal is taken to obtain the adjacency weight of the well pair. That is, the closer the two wells are, the greater the adjacency weight. The adjacency weight is filled into the corresponding position in the initial adjacency matrix, and the weight of each monitoring well and itself is set to 1 to complete the self-loop addition. The adjacency matrix after adding the self-loop is symmetric to make the two elements corresponding to any well pair have the same value. Based on the sum of the weights of each row of the adjacency matrix, the adjacency matrix is ​​symmetrically normalized to make the numerical scales in the row direction and column direction consistent, thus obtaining the graph adjacency matrix.

[0108] The adjacency weight can be calculated using the following formula:

[0109]

[0110] in, For monitoring wells With monitoring well Adjacency weights between them The spatial distance between the two monitoring wells. The feature distance is denoted as .

[0111] The adjacency matrix after adding self-loops is symmetrically normalized to obtain the graph adjacency matrix. :

[0112]

[0113] in, To add the adjacency weight matrix after adding self-loops, It is the identity matrix. It is a diagonal matrix and its diagonal elements are matrices. and The sum of the elements of each row of the sum, diagonal matrix The inverse square root matrix, This is the graph adjacency matrix after symmetric normalization.

[0114] In some embodiments, step S2, "constructing a time-series input sample based on the multi-well spatiotemporal sequence data", includes the following steps:

[0115] S24. Extract the multi-well water level data corresponding to a first preset number of consecutive monitoring times from the multi-well spatiotemporal sequence data as a sub-input sample. Determine the multi-well water level data corresponding to a second preset number of consecutive monitoring times after the first preset number of monitoring times as the multi-well water level prediction target corresponding to the sub-input sample.

[0116] The first preset quantity refers to the number of consecutive historical monitoring moments used for prediction, such as 14. Multi-well water level data refers to the set of water level values ​​of all monitoring wells in the target area at the same monitoring moment. Sub-input sample refers to a segment of multi-well water level data of consecutive historical monitoring moments extracted from multi-well spatiotemporal sequence data, which serves as the input basis for a prediction. The second preset quantity refers to the number of future consecutive monitoring moments to be predicted, such as 14. The multi-well water level prediction target refers to the actual water level data of multi-wells for the second preset quantity of consecutive monitoring moments immediately following the historical monitoring moments covered by the sub-input sample.

[0117] Understandably, predicting future water levels cannot rely solely on data from a single moment in time. A single-moment water level value lacks information about trends and cannot determine whether the level is rising or falling, or at what rate. A continuous historical data set is needed as input, allowing subsequent steps to analyze the trends and patterns of water level changes. Furthermore, groundwater level changes exhibit continuity and inertia; there is a correspondence between water levels at future monitoring moments and those at past monitoring moments. Therefore, a continuous historical water level data set can be paired with a subsequent continuous future water level data set, forming a complete input and prediction target. This allows the model to learn the mapping relationship from history to the future. Moreover, a direct multi-step prediction approach—using water level data from multiple consecutive historical monitoring moments as input and outputting predicted water levels for multiple consecutive future monitoring moments at once—rather than repeatedly predicting and feeding back input moment by moment, avoids the gradual accumulation and amplification of errors during moment-by-moment prediction, resulting in more stable predictions for multiple future monitoring moments.

[0118] Specifically, starting from the beginning of the multi-well spatiotemporal sequence data, the water level data of the multi-wells for a first preset number of monitoring times is extracted as the first sub-input sample. Taking the first preset number as an example, it can be 14, and the second preset number can also be 14. The water level data of the multi-wells from the 1st monitoring time to the 14th monitoring time is extracted as the sub-input sample. The actual water level data of the multi-wells from the 15th monitoring time to the 28th monitoring time is determined as the multi-well water level prediction target corresponding to the sub-input sample. The sub-input sample provides historical water level change information for 14 consecutive monitoring times, and the multi-well water level prediction target provides the actual water level for the next 14 consecutive monitoring times. The two are paired to form a complete training or evaluation data unit.

[0119] S25. According to the order of monitoring time, the original monitoring time interval corresponding to the multi-well spatiotemporal sequence data is divided into training time interval, verification time interval and test time interval.

[0120] The original monitoring time interval refers to the complete time range from the start time to the end time covered by the multi-well spatiotemporal sequence data. The training time interval refers to the time period defined from the original monitoring time interval for training the spatiotemporal prediction network parameters. The validation time interval refers to the time period defined from the original monitoring time interval for monitoring the performance of the spatiotemporal prediction network on unseen data and determining the optimal network parameters during the training process. The testing time interval refers to the time period defined from the original monitoring time interval for the final accuracy evaluation of the trained spatiotemporal prediction network.

[0121] Understandably, spatiotemporal prediction networks need to be trained repeatedly on large amounts of data to achieve good predictive capabilities. The prediction performance after training cannot be judged solely based on the performance on the training data, because the network may fit well on the training data but perform poorly on unseen new data. Therefore, the original data can be divided into three non-overlapping parts in time: the data in the training time interval is used to adjust the network parameters; the data in the validation time interval is used to test the network's adaptability to new data during training and select the optimal parameters accordingly; and the data in the test time interval is completely independent of the training and validation processes and is only used for final evaluation. The division is based on the chronological order of the monitoring time, rather than random division, in order to maintain the continuity of the time series data and make the training process simulate the way past data is used to predict future data in actual applications, avoiding the information leakage caused by future data being mixed into the training process in time.

[0122] Specifically, according to the chronological order of monitoring time from morning to night, the original monitoring time intervals corresponding to the multi-well spatiotemporal sequence data are divided into training time intervals, verification time intervals, and test time intervals. For example, the first 70% of the original monitoring time interval can be used as the training time interval, the next 15% as the verification time interval, and the last 15% as the test time interval. The three time intervals are adjacent to each other on the time axis and do not overlap.

[0123] S26. Sub-input samples are constructed in the training time interval, verification time interval and test time interval respectively. The monitoring time sequence of the multi-well spatiotemporal sequence data is shifted backward by one monitoring time. The sub-input sample extraction and multi-well water level prediction target determination operations are repeatedly performed to obtain multiple sets of sub-input samples and their corresponding multi-well water level prediction targets.

[0124] Understandably, a single set of sub-input samples and the corresponding multi-well water level prediction targets are insufficient to train a spatiotemporal prediction network. Multi-well spatiotemporal sequence data typically covers a long time range. If only one set of samples is extracted from each time interval, a large amount of available information from other time periods within that interval is wasted. By successively shifting one monitoring time forward along the chronological order of monitoring time and repeating the extraction operation, as many sample sets as possible can be extracted from a time interval, allowing the network to fully learn the water level change patterns. At the same time, the monitoring times covered by the multi-well water level prediction targets of the samples constructed in different time intervals do not overlap, thus ensuring that the validation samples and test samples correspond to strictly independent future time periods and avoiding the evaluation results being overly optimistic due to the overlap of prediction targets in time.

[0125] Specifically, for the training time interval, starting from its initial monitoring time, the first sub-input sample is extracted and the corresponding multi-well water level prediction target is determined in the same manner as in S24. Taking a first preset quantity of 14 and a second preset quantity of 14 as an example, starting from the first monitoring time of the training time interval, data from the 1st to the 14th monitoring time are extracted as the first sub-input sample, and data from the 15th to the 28th monitoring time are used as its multi-well water level prediction target. Then, the extraction start position is shifted backward by one monitoring time. Starting from the second monitoring time, data from the 2nd to the 15th monitoring time are extracted as the second sub-input sample, and data from the 16th to the 29th monitoring time are used as its multi-well water level prediction target. This process is repeated, shifting backward by one monitoring time and repeating the extraction and determination operations until the shifted sub-input sample and multi-well water level prediction target are all completely within the training time interval. For the verification time interval and the test time interval, the same method is used to construct multiple sets of sub-input samples and their corresponding multi-well water level prediction targets in their respective intervals.

[0126] S27, based on multiple sets of corresponding sub-input samples and multi-well water level prediction targets, constitutes time-series input samples.

[0127] Specifically, all sub-input samples obtained within the training time interval and their corresponding multi-well water level prediction targets are arranged in order of monitoring time from morning to evening within that interval to form training samples. All sub-input samples obtained within the validation time interval and their corresponding multi-well water level prediction targets are arranged in order of monitoring time to form validation samples. All sub-input samples obtained within the test time interval and their corresponding multi-well water level prediction targets are arranged in order of monitoring time to form test samples. The training samples, validation samples, and test samples together constitute the time-series input samples.

[0128] S3. Construct a spatiotemporal prediction network. Input the time-series input samples into the spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future. The spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer.

[0129] The spatiotemporal prediction network refers to a neural network model used to simultaneously extract spatial correlation features between monitoring wells and temporal evolution features of individual wells. It includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer. The spatial feature extraction layer is used to perform neighborhood aggregation on the water level data of each monitoring well based on the graph adjacency matrix to extract spatial correlation features between wells. The temporal feature extraction layer is used to perform temporal modeling on the time series of each monitoring well to extract the temporal features of water level evolution over time. The multi-step output layer is used to map the extracted temporal features into a water level prediction sequence for each monitoring well at multiple consecutive monitoring times in the future.

[0130] It is understandable that groundwater level changes are affected by both spatial and temporal factors: spatially, the water level change of one monitoring well will propagate through the aquifer and affect neighboring monitoring wells, and there is a hydraulic connection between the water levels of each well. Temporally, water level changes have continuity, trends, and a certain periodicity. The water level evolution process over many days in history contains information for predicting future water levels. Therefore, it is difficult to fully characterize the spatiotemporal evolution of groundwater levels by only modeling spatial correlation or only modeling temporal evolution. This embodiment links spatial feature extraction and temporal feature extraction. First, the spatial feature extraction layer aggregates information between wells in the spatial dimension. Then, the temporal feature extraction layer models the evolution law in the temporal dimension. Finally, the multi-step output layer outputs multi-step prediction results at once, taking into account both spatial correlation and temporal dynamics.

[0131] In some embodiments, step S3 can be implemented as follows:

[0132] S31. Based on the spatial feature extraction layer, the graph adjacency matrix, and the attention mechanism, the aggregation weight is determined, and the water level data of each monitoring well corresponding to each monitoring time are aggregated in the neighborhood to obtain multi-well spatial features.

[0133] Among them, the attention mechanism refers to a learnable weight allocation method that can adaptively calculate the correlation between each monitoring well and its neighboring monitoring wells based on the characteristics of each monitoring well at the current time. Neighbors with high correlation are given greater weights, and neighbors with low correlation are given smaller weights. The aggregation weight refers to the coefficient calculated by the attention mechanism to characterize the strength of the correlation between each monitoring well and its neighboring monitoring wells. The larger the aggregation weight, the greater the influence of the neighbor on the current monitoring well. Neighborhood aggregation refers to the operation of weighting and summing the features of each monitoring well within the range of the non-zero element positions in the graph adjacency matrix according to the corresponding aggregation weights, and then fusing the weighted sum with the monitoring well's own features. Multi-well spatial features refer to the feature expression after fusing the water level information of the monitoring well itself and the water level information of its neighboring monitoring wells. Each monitoring time corresponds to a set of multi-well spatial features.

[0134] Understandably, at a given monitoring moment, the water level of a monitoring well contains its own information and is also influenced by other monitoring wells within the graph's adjacency range. The graph adjacency matrix defines the neighborhood range of each monitoring well through the positions of its non-zero elements, i.e., which neighboring monitoring wells each monitoring well can interact with. However, the numerical weights in the adjacency matrix are not directly used as aggregation weights. The attention mechanism adaptively calculates the aggregation weights corresponding to each neighboring monitoring well within the neighborhood range defined by the graph adjacency matrix, based on the actual characteristics of each monitoring well at the current moment. If the characteristics of a neighboring monitoring well are highly correlated with the characteristics of the current monitoring well, a larger aggregation weight is assigned; otherwise, a smaller aggregation weight is assigned. This ensures that the spatial features aggregated by each monitoring well are a fusion result that is dynamically adjusted based on actual data.

[0135] Specifically, for each monitoring time, the water level data of each monitoring well at that time is used as the input feature. First, a linear transformation is performed to map the original water level values ​​to a higher-dimensional feature space. Then, based on the range of neighboring monitoring wells determined by the non-zero element positions in the graph adjacency matrix, an attention mechanism is used to calculate the attention coefficient between each monitoring well and its neighboring monitoring wells. This attention coefficient is the aggregation weight. The features of each neighboring monitoring well are weighted and summed according to the aggregation weight. The weighted sum is added to the feature of the monitoring well itself. Then, a nonlinear activation process is performed, and the activated result is randomly deactivated to obtain the feature of each monitoring well after fusing spatial information at that monitoring time, i.e., the multi-well spatial feature. The above process is performed for each monitoring time to obtain the multi-well spatial feature sequence corresponding to each monitoring time.

[0136] Taking a target area with 3 monitoring wells, where the water level of well 1 is 10.0 meters, well 2 is 5.7 meters, and well 3 is 3.3 meters at a certain monitoring time as an example, after linear transformation, each of the three wells obtains a feature vector. The positions of the non-zero elements in the graph adjacency matrix limit the neighborhood range of well 1 to well 2 and well 3, the neighborhood range of well 2 to well 1 and well 3, and the neighborhood range of well 3 to well 1 and well 2. The attention mechanism takes the feature vector of each well as input and calculates six aggregation weights for well 1 to well 2, well 1 to well 3, well 2 to well 1, well 2 to well 3, well 3 to well 1, and well 3 to well 2, respectively. Each weight reflects the degree of influence of the corresponding neighboring monitoring well on the current monitoring well. The features of the neighboring monitoring wells are weighted and summed according to the aggregation weights and then added to the features of the well itself to obtain the multi-well spatial features of each monitoring well at that monitoring time.

[0137] In some embodiments, step S31 may be implemented in the following ways:

[0138] S311, the water level data of each monitoring well at each monitoring time are sequentially linearly transformed, and the aggregation weight is determined based on the graph adjacency matrix and attention mechanism. Then, neighborhood feature aggregation, residual connection and nonlinear activation processing are performed to obtain the multi-well spatial features corresponding to each monitoring time.

[0139] Linear transformation refers to the process of multiplying the input water level data of each monitoring well by a learnable weight matrix and adding a bias vector, thereby mapping the original water level values ​​to a higher-dimensional feature space. Neighborhood feature aggregation refers to the operation of weighting and summing the features of each neighboring monitoring well after linear transformation according to the aggregation weight determined by the attention mechanism, based on the range of neighboring monitoring wells defined by the non-zero elements in the graph adjacency matrix, and then fusing the weighted sum with the features of the monitoring well itself. Residual connection refers to the operation of directly adding the features obtained after neighborhood feature aggregation to the input features after linear transformation. Its purpose is to retain the original input information when the network is deepened, avoid excessive dilution of the original information due to multiple transformations, and make the network training more stable. Nonlinear activation processing refers to the operation of applying a nonlinear activation function, such as the ReLU function, to the features after residual connection. Its purpose is to introduce nonlinear expressive power, enabling the network to learn the complex nonlinear spatial relationships in the water level data.

[0140] Understandably, the original water level data of each monitoring well is a low-dimensional, single value. This value alone is insufficient to fully express the spatial interrelationships between the wells. Linear transformation first maps the low-dimensional water level values ​​to a high-dimensional feature space, giving each monitoring well a richer feature vector. Based on this, neighborhood feature aggregation uses the neighborhood range defined by the graph adjacency matrix to sum the features of each neighboring monitoring well according to attention weights, so that the features of each monitoring well include information from neighboring wells. Residual connection adds the aggregated features to its own original features, ensuring that the network does not lose its own water level status while aggregating neighbor information. Nonlinear activation processing introduces nonlinearity into the entire feature extraction process, enabling the network to fit more complex inter-well correlation patterns than simple linear weighting. The above operations are chained together to form the complete processing flow of the spatial feature extraction layer for data at a monitoring time. Repeating this process for each monitoring time yields a sequence that changes with the monitoring time.

[0141] Specifically, for each monitoring time, the water level data of each monitoring well at that time is organized into an input vector. First, a linear transformation is performed to obtain the feature vector corresponding to each monitoring well. Then, based on the positions of non-zero elements in the graph adjacency matrix, the set of neighboring monitoring wells for each monitoring well is determined. An attention mechanism is used to calculate the aggregation weight between each monitoring well and its neighboring monitoring wells. The feature vectors of each neighboring monitoring well are weighted and summed according to the aggregation weight to obtain the neighborhood aggregation result. The neighborhood aggregation result is residually connected with the feature vector of the monitoring well itself, i.e., element-wise addition. A nonlinear activation function is applied to the summed result, and then random deactivation is performed to obtain the multi-well spatial features corresponding to each monitoring well at that monitoring time. The above operation is performed for each monitoring time to obtain the multi-well spatial feature sequence corresponding to each monitoring time.

[0142] S32, Based on the time feature extraction layer, the spatial features of the multi-well corresponding to each monitoring well are time-series modeled according to the order of monitoring time to obtain the time-series feature information of each monitoring well.

[0143] Temporal modeling refers to the process of sequentially inputting the spatial features of each monitoring well at multiple consecutive monitoring times into the temporal feature extraction layer in the order of monitoring time from early to late. The temporal feature extraction layer processes and accumulates historical information moment by moment. Temporal feature information refers to the fixed-length feature expression output by the temporal feature extraction layer after processing the features of a monitoring time and combining it with the previously accumulated historical information. This feature expression is the network's internal summary of the current and previous water level changes, which includes the trend information of water level evolution over time.

[0144] It is understandable that the spatial feature extraction layer obtains the spatial features of each monitoring time. The spatial features of each monitoring time only reflect the spatial relationship between each monitoring well at that time. The spatial features of a single time do not contain the trend of water level change over time. To determine the evolution law in the time dimension, it is necessary to process the spatial features of multiple consecutive monitoring times in chronological order. The time feature extraction layer can adopt a recurrent network structure, such as a long short-term memory network, which has a gated memory mechanism. When processing features time by time, it can selectively retain the trend of water level change and forget short-term random fluctuations. After processing the features of each monitoring time, the network updates its internal memory. When the features of the last historical monitoring time are processed, the network memory accumulates a complete understanding of the entire historical change process.

[0145] Specifically, the multi-well spatial features corresponding to 14 consecutive monitoring times for each monitoring well are sent to the time feature extraction layer in chronological order from morning to night. Taking 3 monitoring wells and 14 days of historical data as an example, each well has one multi-well spatial feature for each day. The time feature extraction layer processes each well separately. For well 1, the multi-well spatial features from day 1 to day 14 are input into the time feature extraction layer sequentially. When processing the features of day 1, the network outputs the internal state of day 1 based on the initial memory and updates the memory. When processing the features of day 2, it outputs the internal state of day 2 based on the updated memory and updates the memory again. This process is repeated day by day until the features of day 14 are processed. At this point, the network memory has accumulated all the information on the water level evolution of well 1 over the past 14 days, including the direction, magnitude, and trend of change. The time-series feature information corresponding to well 1 is output. The same processing is performed on wells 2 and 3 to obtain their respective time-series feature information. The time-series feature information of each monitoring well constitutes the input of the subsequent multi-step output layer.

[0146] In some embodiments, step S32 may be implemented as follows:

[0147] S321, input the multi-well spatial characteristics corresponding to each monitoring well at a first preset number of consecutive monitoring times into a loop network in chronological order of monitoring times, and output the temporal characteristic information of each monitoring well at the last historical monitoring time based on the loop network.

[0148] Among them, recurrent networks refer to a type of neural network with a recurrent connection structure. It contains hidden states that are passed over time and can update the hidden states and retain historical information at each time step when processing sequential data. Examples include long short-term memory networks or gated recurrent units. Historical monitoring time refers to a first preset number of consecutive monitoring times contained in the sub-input sample, such as 14 consecutive monitoring times. The water level data of these monitoring times are known and used to predict future water levels.

[0149] Understandably, the spatial feature extraction layer has generated corresponding multi-well spatial features for each monitoring moment. However, these features are independent of each other in time. Looking at the spatial features of a single monitoring moment alone cannot reveal the direction of water level evolution over time. The recurrent network transmits the hidden state moment by moment in the time dimension, continuously accumulating the information from previous monitoring moments to the current moment. After processing the last historical monitoring moment, the hidden state of the recurrent network has integrated the trend information of the entire historical sequence. Therefore, this hidden state can be used as the temporal feature information corresponding to the monitoring well for subsequent multi-step output layers to predict future water levels. The reason for using a recurrent network instead of simply averaging the features of all moments is that water level changes have a sequential relationship, and the recurrent network can capture this sequential dependency.

[0150] Specifically, taking a first preset number of 14 as an example, for each monitoring well, the multi-well spatial features corresponding to the 1st to 14th historical monitoring times are sequentially input into the recurrent network in the order of monitoring times from early to late. When processing the features of the 1st monitoring time, the recurrent network calculates and updates the hidden state based on the initial hidden state. When processing the features of the 2nd monitoring time, it calculates and updates the hidden state again by combining the hidden state updated at the previous time with the current features. By performing time-by-time processing, until the features of the 14th monitoring time are processed, the hidden state output by the recurrent network is the temporal feature information of the monitoring well at the last historical monitoring time. The above processing is performed on each monitoring well to obtain the temporal feature information corresponding to each monitoring well.

[0151] S33, based on the multi-step output layer and the temporal characteristic information of each monitoring well, output the water level prediction sequence of each monitoring well for a second preset number of consecutive monitoring times in the future.

[0152] Specifically, the temporal characteristic information of each monitoring well at the last historical monitoring moment is input into the multi-step output layer. The linear mapping unit inside the multi-step output layer maps it into a second preset number of values, which correspond to the water level changes at the next second preset number of monitoring moments. Taking the second preset number as 14 as an example, the multi-step output layer outputs 14 values, which represent the changes of the monitoring well relative to the last known water level from day 1 to day 14 in the future. The 14 water level changes are added to the known water level of the monitoring well at the last historical monitoring moment to obtain the water level prediction sequence of the monitoring well for the next 14 consecutive monitoring moments.

[0153] Taking a target area with 3 monitoring wells and using a continuous 14-day forecast for the next 14 days as an example, the spatial feature extraction layer first performs neighborhood aggregation on the water level data of the 3 monitoring wells for each day to obtain the multi-well spatial features that integrate the spatial information between wells for each day. The temporal feature extraction layer performs temporal modeling on each monitoring well along the 14-day time sequence to obtain the temporal feature information of each monitoring well on the 14th day, which is the last historical monitoring moment. The multi-step output layer maps the temporal feature information of well 1 into 14 water level changes, and then adds them to the known water level of well 1 on the 14th day item by item to obtain the daily water level prediction sequence of well 1 for the next 14 days. Wells 2 and 3 are processed in the same way to finally obtain the daily water level prediction sequence of the 3 monitoring wells for the next 14 days.

[0154] In some embodiments, step S33 may be implemented in the following ways:

[0155] S331 inputs the time-series characteristic information of each monitoring well at the last historical monitoring time into the linear mapping unit, and directly outputs the water level change of each monitoring well at the next second preset number of monitoring times.

[0156] The linear mapping unit refers to a computational unit composed of a learnable weight matrix and a bias vector. Its function is to directly map the input time-series feature information into an output vector of a specified length through a matrix multiplication plus a bias operation. The water level change refers to the change in the known water level of each monitoring well at each future monitoring time relative to the last historical monitoring time. Positive values ​​indicate a rise in water level, and negative values ​​indicate a fall in water level. The unit of its value is consistent with the original water level data.

[0157] It is understandable that the time-series feature information is an abstract summary of the historical water level change process by the recurrent network. Its dimension is usually inconsistent with the number of monitoring moments that need to be predicted in the future. It needs to be converted into an output vector corresponding to the number of prediction steps. The linear mapping unit takes the time-series feature information as input and directly outputs a vector of length two preset numbers through a linear calculation. Each element in this vector corresponds to the water level change at a future monitoring moment. By using a direct mapping instead of step-by-step iterative prediction, the error introduced by using the prediction result of the previous step as the input of the next step can be avoided, making the prediction results for multiple future monitoring moments more stable.

[0158] Specifically, taking a second preset quantity of 14 as an example, for each monitoring well, its temporal feature information at the last historical monitoring moment is used as an input vector and sent to the linear mapping unit. The linear mapping unit stores a pre-trained weight matrix and bias vector. Through matrix multiplication and bias operation, the input temporal feature information is mapped into an output vector of length 14. The 14 elements in the output vector are the water level changes of the monitoring well in the next 14 consecutive monitoring moments. The above operation is performed on each monitoring well to obtain the water level change sequence corresponding to each monitoring well.

[0159] S332, using the known water level corresponding to the last historical monitoring time in the time-series input sample as a benchmark, the known water level is superimposed with the water level change to obtain the water level prediction sequence for each monitoring well in the next second preset number of monitoring times.

[0160] Specifically, taking a second preset number of 14 as an example, for each monitoring well, the known water level at the last historical monitoring time (i.e., the 14th monitoring time) in the sub-input sample is obtained. The water level change of the monitoring well in the next 14 consecutive monitoring times is added to the known water level in chronological order to obtain the absolute water level prediction value of the monitoring well in the next 14 consecutive monitoring times. The 14 absolute water level prediction values ​​are arranged in chronological order to form the water level prediction sequence of the monitoring well. The above operation is performed on each monitoring well to obtain the water level prediction sequence of each monitoring well in the next second preset number of monitoring times.

[0161] S4. Construct hydrogeological prior constraint residuals based on the water level prediction sequence. The hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term.

[0162] Among them, the hydrogeological prior constraint residual refers to the numerical index that measures the degree to which the water level prediction sequence deviates from the hydrogeological prior conditions. The smaller the value, the better the hydrogeological consistency of the prediction results. The hydrogeological prior constraint residual includes a first constraint term to constrain the water level change of each monitoring well to remain consistent at each monitoring time during the prediction period, a second constraint term to constrain the relative difference of water level between each monitoring well to remain stable, and a third constraint term to constrain the water level evolution trend over time to remain smooth.

[0163] It is understandable that the water level prediction using a spatiotemporal prediction network in step S3 essentially infers future water levels based on statistical patterns in historical data. When training data is insufficient or hydrological conditions differ from those during training, the network may output statistically reasonable but hydrogeologically unreasonable predictions. For example, the water level change rates of each monitoring well may fluctuate drastically and in opposite directions, or the relative differences in water levels between monitoring wells may deviate significantly from their historical statistical relationships. Therefore, three types of hydrogeological priors—regional consistency, inter-well relative differences, and temporal smoothing—can be mathematically incorporated into the prediction process. Constraints are imposed on the prediction results from three perspectives: water level change consistency, inter-well relative differences, and water level evolution smoothness. This ensures that the prediction results meet the data fitting accuracy requirements while exhibiting good hydrogeological consistency. It should be noted that the above three constraints are prior regularization terms constructed based on the statistical characteristics of monitoring data and are not equivalent to rigorous physical equations such as the groundwater flow control equation, Darcy's law, or water balance equation. This embodiment does not claim to rigorously simulate the physical processes of groundwater movement based on these constraints.

[0164] In some embodiments, step S4 can be implemented as follows:

[0165] S41, the water level prediction sequence is spliced ​​with the water level data corresponding to the last historical monitoring moment in the time series input sample to obtain the continuous water level prediction sequence of each monitoring well, and the continuous water level prediction sequence is restored to the water level sequence with physical dimensions.

[0166] Among them, the continuous water level prediction sequence refers to the sequence formed by splicing the known water level at the last historical monitoring moment with the predicted water levels at multiple consecutive future monitoring moments in chronological order, while the water level sequence with physical dimensions refers to the sequence obtained after restoring the normalized values ​​to water level values ​​in meters.

[0167] It is understandable that the spatiotemporal prediction network outputs predicted values ​​within a normalized interval, while the calculation of prior constraints requires water level values ​​with physical dimensions. At the same time, prior constraints need to examine the water level change process within the prediction period, and the known water level at the beginning of the prediction period is an important benchmark for calculation. Therefore, the known water level at the last historical monitoring moment can be concatenated with the predicted water level sequence to form a continuous sequence, and then the normalization parameter can be used to restore the sequence to a water level sequence with physical dimensions, so as to perform constraint calculations based on physical dimensions.

[0168] Specifically, the predicted water level sequence for each monitoring well over the next 14 consecutive monitoring times can be concatenated with the known water level of that well at the last historical monitoring time to obtain a continuous water level prediction sequence of length 15. Using the normalization parameter determined in step S13, the normalized values ​​of each monitoring well in this continuous water level prediction sequence are restored to water level values ​​in meters, resulting in a water level sequence with physical dimensions. The inverse normalization calculation formula is as follows:

[0169]

[0170] in, The first in the continuous water level prediction sequence The monitoring well in the first Normalized values ​​at each monitoring time point For the first The maximum water level value of the monitoring well during the training time interval. For the first The minimum water level value of the monitoring well during the training time interval. For the restored first The monitoring well in the first The formula for the physical dimension water level value at each monitoring time is as follows: First, normalize the value from the interval [-1,1] back to the interval [0,1], then multiply it by the water level fluctuation range of the monitoring well during the training period, i.e., the maximum value minus the minimum value, and finally add the lowest water level value of the monitoring well to obtain a value with the same dimension as the original water level data.

[0171] S42, based on the water level difference between adjacent monitoring times and the water level range of each monitoring well within the training time interval, determine the normalized water level change of each monitoring well. Based on the deviation between the normalized water level change of each monitoring well at each monitoring time and the average normalized water level change of all monitoring wells at that monitoring time, normalize the change by the overall change range of the normalized water level of each monitoring well, and construct the first constraint term.

[0172] Here, the water level difference refers to the difference between the predicted water level at the next monitoring time and the actual water level at the current monitoring time. The normalized water level change is a dimensionless value obtained by dividing the water level difference between adjacent monitoring times by the water level range within the training time interval of the corresponding monitoring well. The water level range refers to the difference between the maximum and minimum water level values ​​within the training time interval of the monitoring well. Adjacent monitoring times refer to two monitoring times that are adjacent according to the monitoring frequency, such as two consecutive days. To prevent the calculation from being unstable due to an excessively small range, a lower limit ε is set for the water level range. s , ε s =1×10 -3 m, that is, when the range of values ​​is less than ε s Press ε s calculate.

[0173] It is understandable that monitoring wells within the same region are situated within the same aquifer system and are affected by shared rainfall and groundwater extraction. Their water level changes are typically synchronous; for example, after a rainfall event, the water levels in several monitoring wells within the region will rise, rather than one well rising while another falls. If the water level change of a well output by the spatiotemporal prediction network differs significantly from the average change of other wells in the same region, it indicates that the prediction result violates the hydrogeological principle of regional consistency. However, the water level fluctuation range varies among different monitoring wells; for example, some wells may experience fluctuations of up to 10 meters. Some wells have a water level difference of only 1 meter. Even with the same change of 0.1 meters, the change is significant for a well with a fluctuation range of 1 meter, but very small for a well with a fluctuation range of 10 meters. Therefore, it is necessary to first divide the water level difference by the water level range of the monitoring well itself, converting the change into a proportion relative to its own fluctuation range, eliminating the influence of the fluctuation range differences between different monitoring wells, and ensuring that the constraint is applied fairly to each monitoring well. Then, the average value of the normalized water level change of all monitoring wells at the same monitoring time is used as the benchmark to measure the degree to which the change of each monitoring well deviates from the average value. The greater the deviation, the worse the regional consistency.

[0174] Specifically, the water level difference between adjacent monitoring times for each monitoring well can be calculated. This difference is then divided by the water level range within the corresponding training time interval of the monitoring well, with the lower limit of this range being εs. This yields the normalized water level change between adjacent monitoring times for each monitoring well. For each predicted monitoring time, the average of the normalized water level changes of all monitoring wells at that time is calculated as the average normalized water level change for that monitoring time. The square of the difference between the normalized water level change of each monitoring well at that monitoring time and its corresponding average normalized water level change is calculated and averaged over all monitoring wells and all monitoring times to obtain the degree of deviation. The square mean of the normalized water level changes of all monitoring wells is calculated as the overall change amplitude. To prevent the overall change amplitude from being too small and resulting in division by zero, a lower limit is set for it. , =1×10 -8 Furthermore, this overall change magnitude is used only as a normalization constant and does not participate in the gradient calculation during backpropagation. The deviation is divided by the overall change magnitude and then by the reference scale determined within the training time interval to obtain the first constraint term. The relevant calculation formula is as follows:

[0175]

[0176]

[0177]

[0178] in, For the first The monitoring well in the first Normalized water level change at each monitoring time For the first The monitoring well in the first Water level values ​​at each monitoring time. For the first The monitoring well in the first Water level values ​​at each monitoring time. For the first The range of water level values ​​for the monitoring well during the training period. For the first The maximum water level value of the monitoring well during the training time interval. For the first The minimum water level value of the monitoring well during the training time interval. This represents the lower limit of the water level range. The total number of monitoring wells, To predict the total number of monitoring times, For the first The monitoring well in the first Normalized water level change at each monitoring time For the first The average normalized water level change of all monitoring wells at each monitoring time point. This represents the lower limit of the overall change range. This is the first constraint term.

[0179] In the above formula, This measures the degree to which the changes in each monitoring well deviate from the average changes in the same area; the greater the deviation, the worse the regional consistency. It measures the overall magnitude of change and is used to convert the numerator into a relative value. For example, if all monitoring wells change synchronously, the change of each monitoring well is close to the average value, the numerator is close to zero, and the first constraint term is close to zero, indicating that the regional consistency is very good. If the change of a certain well deviates greatly from the average value, the numerator increases, the first constraint term increases, indicating that the regional consistency is poor.

[0180] S43. Based on the difference in the average water level of each monitoring well within the training time interval, determine the relative difference benchmark between each monitoring well pair. Based on the deviation between the water level difference of each monitoring well pair in the continuous water level prediction sequence and the relative difference benchmark between the monitoring well pair, construct the second constraint term.

[0181] Among them, the difference in average water level refers to the difference between the average water level of one monitoring well and the average water level of another monitoring well during the training time interval. The inter-well head difference benchmark refers to the historical statistical benchmark of the inter-well head difference determined by the difference in the average head of two monitoring wells under a unified elevation benchmark during the training time interval. Since the modeling variables of each monitoring well have been uniformly converted to a unified elevation benchmark, this benchmark reflects the average head difference between the two monitoring wells and has a physical meaning under a unified elevation benchmark.

[0182] It is understandable that when two monitoring wells are located in the same area, their water levels are affected by common recharge, discharge, and aquifer conditions. The relative difference in water levels between the two wells usually fluctuates around a certain historical average level. For example, if the average water level difference between well A and well B during the training period is 3 meters, then during the prediction period, the water level difference between well A and well B should also be maintained at approximately 3 meters. Even if the overall groundwater level rises or falls, the water level difference between the two wells should remain relatively stable because they are two fixed observation points in the same aquifer. The water level difference is mainly determined by the spatial location and elevation of the two points and will not change drastically over time. If the predicted water level difference between two monitoring wells deviates significantly from the historical average level, it indicates that the prediction result may have disrupted the statistical consistency of the water level relationship between the wells. Therefore, constraint terms can be constructed based on the deviation between the water level difference of each monitoring well pair and the benchmark of the relative difference between the wells to maintain the relative difference relationship between the wells in the prediction result.

[0183] Specifically, the mean water level value of each monitoring well within the training time interval is calculated. The difference between the mean values ​​of two monitoring wells is used to obtain the inter-well head difference benchmark for that well pair. The mean of the water level range of all monitoring wells within the training time interval is calculated as the inter-well difference normalization scale. For any two monitoring wells in the continuous water level prediction sequence, the water level difference between them at each predicted monitoring time is calculated. This water level difference is subtracted from the corresponding inter-well head difference benchmark, and then divided by the inter-well difference normalization scale to obtain the deviation value. The average of the squared deviation values ​​of all well pairs and all predicted monitoring times is taken, and then divided by the reference scale determined within the training time interval to obtain the second constraint term. The relevant calculation formulas are as follows:

[0184]

[0185]

[0186] in, For the first Monitoring well and the first The benchmark for the head difference between monitoring wells is determined by the difference in the average water level of the two monitoring wells during the training time interval, reflecting the average head difference under a unified elevation benchmark. For the first The average water level of the monitoring well during the training time period. For the first The average water level of the monitoring well during the training time period. For the first The monitoring well in the first Water level values ​​at each monitoring time. For the first The monitoring well in the first Water level values ​​at each monitoring time. The average value of the water level range across all monitoring wells during the training time interval. To monitor the total number of well pairs, To predict the total number of monitoring times, summation index. This indicates that the summation is performed only over different pairs of monitoring wells. This is the second constraint term.

[0187] The above formula means that for each pair of monitoring wells, at each prediction time, the water level difference between them is calculated, and then the historical average water level difference during the training period is subtracted. The difference reflects the degree to which the current predicted water level relationship between wells deviates from the historical level. The second constraint term is obtained by averaging the squares of the deviations of all well pairs and all times. For example, if the average water level difference between well A and well B during the training period is 3 meters, and the water level difference every day during the prediction period is also around 3 meters, then the deviation is close to zero, and the second constraint term is close to zero, indicating that the water level relationship between wells has maintained historical consistency. If the water level difference on a certain day during the prediction period becomes 5 meters, then the deviation is 2 meters, and the second constraint term increases, indicating that the water level relationship between wells has deviated abnormally.

[0188] S44. Based on the difference results of the values ​​of adjacent monitoring times of each monitoring well in the continuous water level prediction sequence, the third constraint term is constructed.

[0189] Among them, the difference result refers to the difference between the numerical changes between two adjacent monitoring times, that is, the second difference, which is used to measure the degree of fluctuation of the water level change trend itself.

[0190] It is understandable that the change process of groundwater level is not only stable in itself, but its trend is also usually gradual. For example, the water level may drop slowly for several days in a roughly stable rate, rather than fluctuating repeatedly like dropping 0.1 meters today, rising 0.1 meters tomorrow, and dropping 0.1 meters the day after. The second difference measures the degree of fluctuation in the water level change trend itself, that is, the difference between the water level changes in two adjacent time periods. If the water level change trend is smooth, the changes in adjacent time periods are roughly equal, and the second difference is close to zero. If the water level change trend fluctuates sharply, the changes in adjacent time periods are significantly different, and the amplitude of the second difference increases. Therefore, constraints can be imposed on the second difference of the continuous water level prediction sequence of each monitoring well to keep the water level evolution trend smooth.

[0191] Specifically, for each monitoring well, the difference between the values ​​at adjacent monitoring times in the continuous water level prediction sequence is calculated, and then the difference between adjacent differences is calculated again to obtain the second-order difference result for that monitoring well. Taking the values ​​at times t, t+1, and t+2 in the continuous water level prediction sequence as an example, first calculate the difference between time t+1 and time t, then calculate the difference between time t+2 and time t+1, and then subtract these two differences to obtain the second-order difference result corresponding to time t. The average of the squares of the second-order difference results of all monitoring wells is then divided by the reference scale determined within the training time interval to obtain the third constraint term. The relevant calculation formula is as follows:

[0192]

[0193] in, The first in the continuous water level prediction sequence The monitoring well in the first Normalized values ​​at each monitoring time point For the first The monitoring well in the first Normalized values ​​at each monitoring time point For the first The monitoring well in the first Normalized values ​​at each monitoring time point This represents the total number of all monitored wells and the predicted time. The expression within the parentheses is the third constraint term. This is the formula for calculating the second-order difference.

[0194] S45, determine the reference scales corresponding to the first constraint term, the second constraint term, and the third constraint term based on the water level data within the training time interval, divide the first constraint term, the second constraint term, and the third constraint term by their respective reference scales, and perform a weighted summation of the first constraint term, the second constraint term, and the third constraint term after reference scale normalization according to preset constraint weights to obtain the hydrogeological a priori constraint residuals.

[0195] Among them, the preset constraint weights refer to the coefficients of the proportion of each constraint in the hydrogeological prior constraint residuals. For example, the weights of the three constraints are set to 0.40, 0.30 and 0.30 respectively. The reference scale refers to the magnitude benchmark determined according to the calculation method of the corresponding constraint item based on the water level data in the training time interval, which is used to make the magnitude of each constraint consistent.

[0196] Understandably, the three constraints characterize the hydrogeological consistency of the prediction results from three perspectives: water level change consistency, inter-well relative difference, and time smoothness. Their impact on the prediction results varies. By setting the weights of each constraint, their relative importance can be adjusted. For example, if water level change consistency is considered the most important, the weight of the first constraint can be set to the maximum. Furthermore, although all three constraints are dimensionless statistics, their calculation methods and numerical ranges differ, resulting in significant differences in magnitude. For instance, the first constraint, after normalization of the overall change magnitude, may be between 0 and 1; the second constraint, after normalization of the inter-well difference scale, may be between 0 and a few values; and the third constraint, based on normalized numerical calculations, may fall within a very small range. If weighted directly without processing, the contribution of these small-magnitude constraints to the sum would be negligible, failing to exert a binding effect. Therefore, each constraint is first divided by a reference scale obtained using the same calculation method within the training time interval to ensure consistency in magnitude before weighted summation, thus allowing the three constraints to function together in a balanced manner.

[0197] Specifically, the first constraint term after reference scale normalization is multiplied by 0.40, the second constraint term by 0.30, and the third constraint term by 0.30. The three weighted results are then summed to obtain the hydrogeological a priori constraint residual. The reference scale normalization calculation formula is as follows:

[0198]

[0199] in, As a reference to the k-th constraint after scale normalization, For the k-th constraint calculated according to the aforementioned formula, The reference scale is calculated based on the water level data within the training time interval using the corresponding formula. k takes values ​​of 1, 2, and 3, which correspond to the first constraint term, the second constraint term, and the third constraint term, respectively.

[0200] S5. Based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, and the hydrogeological prior constraint residual, a joint loss function is constructed. Based on the joint loss function, the spatiotemporal prediction network is subjected to constraint fine-tuning training to obtain the trained spatiotemporal prediction network.

[0201] Among them, data prediction loss refers to the magnitude of the deviation between the water level prediction sequence and the multi-well water level prediction target, and the joint loss function refers to the function that integrates the data prediction loss and the hydrogeological prior constraint residual into a unified optimization objective.

[0202] It is understandable that an untrained spatiotemporal prediction network has poor predictive ability due to its initial parameters. It needs to be improved through a large number of iterations. The training process is essentially a cyclical iterative process of prediction, checking, adjustment, and re-prediction. However, if data prediction loss and hydrogeological prior constraint residuals are introduced at the beginning of training, the prior constraints may interfere too strongly or too early with the prediction results, making it difficult for the network to simultaneously meet the two objectives of data fitting and prior constraints. Therefore, this embodiment adopts a two-stage training method of pre-training followed by constraint fine-tuning. First, the network is allowed to fully learn the water level change patterns in the data. Then, prior constraints are gradually introduced on the basis of data fitting, so that the network can eventually have both prediction accuracy and hydrogeological consistency.

[0203] In some embodiments, step S5 can be implemented as follows:

[0204] S51. Based on the mean square error between the water level prediction sequence and the multi-well water level prediction target, the data prediction loss is obtained. The data prediction loss is then summed with the product of the preset prior constraint weights and the hydrogeological prior constraint residuals to obtain the joint loss function.

[0205] Among them, mean square error refers to the error index obtained by first calculating the square of the difference between each predicted value and the corresponding target value, and then averaging all the squared values. Prior constraint weight refers to the pre-set coefficient used to control the proportion of hydrogeological prior constraint residuals in the joint loss function.

[0206] Understandably, the data prediction loss reflects the numerical accuracy of the prediction, while the hydrogeological prior constraint residual reflects the hydrogeological consistency of the prediction results. Both evaluate the quality of the prediction results from different perspectives. A smaller data prediction loss indicates that the predicted value is closer to the actual water level; a smaller hydrogeological prior constraint residual indicates that the prediction results better conform to hydrogeological priors such as regional consistency, stability of inter-well relative differences, and temporal smoothness. Weighting and summing the two results yields a unified joint loss function. A smaller joint loss function value indicates better overall performance of the spatiotemporal prediction network. The role of the prior constraint weights is to adjust the importance of hydrogeological consistency in the overall optimization objective. Larger weights make the network more inclined to satisfy prior constraints during training; smaller weights make the network more focused on fitting historical data.

[0207] Specifically, the square of the difference between the predicted value of each monitoring well at each predicted monitoring time and the corresponding multi-well water level prediction target is calculated in the water level prediction sequence. Then, the average of all squared values ​​is calculated to obtain the data prediction loss. The preset prior constraint weights are multiplied by the hydrogeological prior constraint residuals, and then added to the data prediction loss to obtain the joint loss function. The relevant calculation formulas are as follows:

[0208]

[0209]

[0210] in, To predict losses based on data, For the first The actual water level corresponding to each predicted value. For the first One predicted value, The total number of predicted values. For prior constraint weights, For hydrogeological a priori constraint residuals, This is the joint loss function.

[0211] S52, the spatiotemporal prediction network is iteratively trained based on the data prediction loss until the preset pre-training convergence condition is met, and the pre-trained network parameters are obtained.

[0212] Among them, the pre-training convergence condition refers to the criterion for determining whether the pre-training phase can end, such as the loss on the validation samples no longer decreasing for several consecutive rounds, or reaching the preset maximum number of iterations. The pre-trained network parameters refer to the network parameters that are saved at the end of the pre-training phase to minimize the validation loss.

[0213] Understandably, the goal of the pre-training phase is to allow the spatiotemporal prediction network to fully learn the changing patterns in historical water level data. This phase only optimizes the data prediction loss and does not introduce prior constraints. Since there are no prior constraint terms, this phase can use a relatively large learning rate and a long number of training epochs to allow the network parameters to converge fully. During the training process, validation samples can be used to monitor the network's performance on unseen data and save the network parameters with the minimum validation loss. The pre-training phase ends when the validation loss no longer improves for several consecutive epochs or reaches the maximum number of iterations.

[0214] Specifically, training samples are input into the spatiotemporal prediction network in batches according to a preset batch size. The data prediction loss for each batch is calculated, the gradient of the network parameters is calculated based on the data prediction loss, and the network parameters are updated along the loss descent direction to complete one iteration. After each iteration, the validation loss is calculated using validation samples. If the validation loss is improved compared to the historical best value, the current network parameters are saved. If the validation loss does not improve for several consecutive iterations and the number of iterations reaches the preset minimum number of iterations, training is stopped early. When the validation loss no longer improves or reaches the preset maximum number of iterations, the pre-training phase ends, and the saved optimal network parameters are used as the pre-trained network parameters.

[0215] S53, initialize the network parameters of the spatiotemporal prediction network to the pre-trained network parameters, and during the training iteration process, gradually increase the prior constraint weights from the initial weight values ​​to the preset target weight values ​​according to the preset number of warm-up rounds.

[0216] Among them, the preset warm-up rounds refer to the number of iteration rounds required for the prior constraint weights to increase from the initial weight value to the target weight value, the initial weight value refers to the value of the prior constraint weights at the beginning of the fine-tuning phase, and the target weight value refers to the final value of the prior constraint weights after the fine-tuning phase ends.

[0217] Understandably, the goal of the fine-tuning phase is to gradually introduce hydrogeological prior constraints based on the pre-trained network's existing predictive capabilities, thereby improving hydrogeological consistency while maintaining prediction accuracy. If the prior constraint weights are set to the target weight values ​​at the beginning of the fine-tuning phase, the prior constraints will suddenly have a significant impact on the network, potentially causing drastic adjustments to network parameters and disrupting the data patterns learned during the pre-training phase. Therefore, a weight warm-up strategy is adopted. In the first few iterations of the fine-tuning phase, the prior constraint weights are gradually increased from a small initial value to the target weight value, allowing the influence of the prior constraints on the network to be applied gradually from weak to strong. This enables the network to gradually adapt to the prior constraints while maintaining its existing predictive capabilities.

[0218] Specifically, the network parameters of the spatiotemporal prediction network are set as those of the pre-trained network. The target weight values ​​of the prior constraint weights can be determined through parameter tuning on the validation set grid. For example, in the candidate set {0.05, 0.1, 0.2, 0.3, 0.5, 0.7}, a target weight value of 0.7 is selected based on the root mean square error of the 14-day mean of the validation set. A warm-up round is set, for example, 10 rounds. In the first 10 iterations of the fine-tuning phase, the prior constraint weights are linearly increased from the initial weight value of 0.07 to the target weight value of 0.7 according to the number of iterations. After the 10th round, the weights remain at 0.7.

[0219] S54. The network parameters of the spatiotemporal prediction network are iteratively updated based on the joint loss function until the preset fine-tuning convergence condition is met, and the trained spatiotemporal prediction network is obtained.

[0220] Among them, the fine-tuning convergence condition refers to the criterion for determining whether the fine-tuning phase can end.

[0221] Understandably, the fine-tuning phase uses the joint loss function as the optimization objective, while considering both data prediction accuracy and hydrogeological consistency. In each iteration, the network adjusts its parameters based on the gradient of the joint loss function, gradually decreasing the value of the joint loss function. This improves hydrogeological consistency while maintaining prediction accuracy. During training, validation samples are used to monitor the network's performance on unseen data, and the network parameters with the minimum validation joint loss are saved to ensure that the final network has good generalization ability. To prevent excessively large gradient magnitudes from causing training instability, the gradients of the network parameters can be pruned. An adaptive moment estimation optimizer can be used to adaptively adjust the step size based on the statistical information of the historical gradients of each parameter, which helps to accelerate convergence and improve training stability.

[0222] Specifically, in each iteration, training samples are input into the spatiotemporal prediction network in batches. The data prediction loss and hydrogeological prior constraint residuals are calculated, a joint loss function is constructed, and the partial derivatives of the joint loss function with respect to each network parameter are calculated. The update amount of each parameter is determined in combination with the preset adjustment step size, and the network parameters are updated in the direction that makes the joint loss function decrease. At the same time, the gradients of the network parameters are clipped to prevent the gradient magnitude from being too large and causing training instability. After each iteration, the validation joint loss, which is composed of the weighted sum of the data prediction loss and the prior constraint residuals, is calculated using validation samples, and the network parameter with the minimum validation joint loss is saved. When the validation joint loss does not improve for several consecutive iterations and reaches the preset minimum number of iterations or the preset maximum number of iterations, the fine-tuning stage ends, and the trained spatiotemporal prediction network is obtained. In this embodiment, an adaptive moment estimation optimizer can be used to iteratively update the network parameters. This optimizer can adaptively adjust the step size according to the statistical information of the historical gradients of each parameter, which helps to accelerate convergence and improve training stability. In each iteration, all network parameters can be updated synchronously based on the gradient obtained by backpropagation.

[0223] S6. Input the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network to obtain the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times.

[0224] Among them, the groundwater level prediction result refers to the predicted water level value of each monitoring well at multiple consecutive monitoring times in the future. This result can be further converted into a water level value with physical dimensions and output.

[0225] Understandably, in practical applications, the known water level data of the target area at recent continuous monitoring times are processed into an input format acceptable to the network according to steps S1 and S2 and input into the trained spatiotemporal prediction network. Then, the predicted water level values ​​of each monitoring well for multiple consecutive monitoring times in the future can be output according to step S3. The accuracy of the output results and the three statistical prior indicators can be evaluated according to the test method described in this embodiment.

[0226] Specifically, known water level data of groundwater monitoring wells in the target area for 14 consecutive monitoring times can be obtained. The data is preprocessed and sample constructed according to steps S1 and S2. The processed data is then input into the trained spatiotemporal prediction network. The network outputs the water level prediction results of each monitoring well for 14 consecutive monitoring times immediately following this known data. The prediction results are then restored to water level values ​​in meters using normalization parameters.

[0227] For example, when it is necessary to understand the water level change trend of each monitoring well in the target area in the next two weeks, the known water level data of each monitoring well for the most recent 14 consecutive days can be obtained. After preprocessing, normalization and sample construction are completed in the manner of steps S1 and S2, the data is input into the trained spatiotemporal prediction network, and the network can output the daily water level prediction value of each monitoring well for the next 14 consecutive days.

[0228] After obtaining the predicted groundwater levels for each monitoring well at multiple consecutive monitoring times in the future, this scheme also includes the following embodiments:

[0229] A1 converts the groundwater level prediction results of each monitoring well into physical dimension water level values, and determines the water level prediction value of each monitoring well in each prediction time domain according to multiple preset prediction time domains.

[0230] The prediction time domain refers to different time spans used to evaluate the prediction accuracy in stages, calculated from the prediction starting point. For example, day 1 is a prediction time domain, days 1 to 3 are a prediction time domain, days 1 to 7 are a prediction time domain, and days 1 to 14 are a prediction time domain. The predicted water level value refers to the predicted water level value of each monitoring well at each monitoring time within a certain prediction time domain.

[0231] It is understandable that spatiotemporal prediction networks output a sequence of water level predictions for multiple consecutive monitoring times in a single operation. However, the difficulty of prediction varies across different time spans. Taking a 14-day prediction as an example, the most available information and the least accumulated error make the prediction less difficult when predicting the water level on day 1. Conversely, the available recent information is relatively limited when predicting the water level on day 14, making the prediction more difficult. If only one overall indicator is used to measure the prediction effect for all 14 days, it is impossible to distinguish the differences in the network's performance in short-term and long-term predictions, nor can it be determined which time period the error mainly originates from. Therefore, by determining the water level prediction values ​​for each monitoring well in each prediction time domain according to multiple preset prediction time domains, the prediction accuracy of the network can be evaluated in stages, providing a more comprehensive reflection of the network's performance across different prediction spans.

[0232] Specifically, using the normalization parameters determined in step S13, the groundwater level prediction results of each monitoring well are restored from normalized values ​​to physical dimension water level values ​​in meters. According to multiple preset prediction time domains, the water level prediction values ​​of each monitoring well in the corresponding time domain are extracted from the restored water level prediction results. For example, if four prediction time domains are preset, namely day 1, day 1 to 3, day 1 to 7, and day 1 to 14, the water level prediction values ​​of each monitoring well in these four time domains are extracted respectively. The prediction values ​​in each time domain constitute a set of prediction values ​​to be evaluated.

[0233] A2. Based on the predicted water level and the actual water level of each monitoring well in each prediction time domain, calculate the root mean square error and coefficient of determination of each monitoring well, and take the macro average of the root mean square error and coefficient of determination of all monitoring wells to obtain the macro average index corresponding to each prediction time domain.

[0234] Among them, the true water level value refers to the actual water level value obtained by actual measurement of each monitoring well at the corresponding monitoring time; the root mean square error is the error index obtained by first calculating the square of the difference between each predicted value and the corresponding true value, averaging all the square values ​​and then taking the square root; it is used to measure the degree of deviation of the predicted value from the true value as a whole. The smaller the value, the more accurate the prediction; the coefficient of determination is an index used to measure the explanatory power of the predicted value on the trend of the change of the true value. The closer the value is to 1, the more consistent the trend of the predicted value with the change of the true value; the macro average index is obtained by first calculating the root mean square error and the coefficient of determination of each monitoring well separately, then averaging the root mean square errors of all monitoring wells to obtain the macro average root mean square error, and averaging the coefficients of determination of all monitoring wells to obtain the macro average coefficient of determination.

[0235] Understandably, if all predicted values ​​from all monitoring wells are combined to calculate a uniform root mean square error, wells with large water level fluctuations (e.g., those fluctuating within a few meters) will dominate the calculation results. Conversely, wells with small water level fluctuations (e.g., those fluctuating only within tens of centimeters), even if accurately predicted, will have their contributions to the overall error overshadowed by the large fluctuations. This results in an evaluation that does not fairly reflect the network's predictive capabilities for different monitoring wells. The macro-averaging approach calculates the index for each monitoring well separately and then averages it across all wells. This ensures that each monitoring well, regardless of the magnitude of water level fluctuation, has the same weight in the final index, allowing for a fairer evaluation of the network's predictive capabilities for different monitoring wells. The root mean square error measures prediction accuracy from the perspective of numerical deviation, while the coefficient of determination measures prediction accuracy from the perspective of trend fit. Combining the two provides a more comprehensive evaluation of the network's predictive performance.

[0236] Specifically, for each prediction time domain, the root mean square error (RMSE) and coefficient of determination (COP) between the predicted water level and the corresponding actual water level at each monitoring time for each monitoring well are calculated. The RMSE is calculated as follows: for each monitoring well, the difference between the predicted and actual water level at each monitoring time is calculated, the squares of the differences are averaged, and then the square root is taken to obtain the RMSE of the well in that time domain. The COP is calculated as follows: for each monitoring well, the average of the actual water level at that time domain is calculated, then the sum of squares of the differences between the predicted and actual water levels, and the sum of squares of the differences between the actual water level and the average of the actual water levels are calculated. The former is subtracted from the latter by 1 to obtain the COP of the well in that time domain. The arithmetic mean of the RMSEs of all monitoring wells is taken to obtain the macro-average RMSE corresponding to that prediction time domain. The arithmetic mean of the COPs of all monitoring wells is taken to obtain the macro-average COP corresponding to that prediction time domain. The macro-average RMSE and the macro-average COP together constitute the macro-average index corresponding to that prediction time domain.

[0237] A3 uses the arithmetic mean of the macro-average indices corresponding to the third preset number of prediction time domains as the overall evaluation index for multi-step prediction, in order to evaluate the prediction accuracy of the spatiotemporal prediction network.

[0238] The third preset quantity refers to the number of preset prediction time domains. For example, if four prediction time domains are preset: day 1, day 1 to 3, day 1 to 7, and day 1 to 14, the third preset quantity is 4. The overall evaluation index is the value obtained by taking the arithmetic mean of the macro average index corresponding to each prediction time domain. It is used to comprehensively evaluate the overall prediction accuracy of the network in the multi-step prediction task. Prediction accuracy refers to the degree of agreement between the water level prediction value output by the network and the actual water level.

[0239] Understandably, the macro-average indices for each prediction time domain reflect the network's performance across different prediction spans. For example, the macro-average root mean square error (RMSE) on day 1 might be smaller, while the macro-average RMSE on day 14 might be larger. This aligns with the rule that the more prediction steps, the greater the error. However, to provide a comprehensive conclusion about the network's overall prediction capability, these indices need to be further aggregated into a single overall evaluation index. Taking the arithmetic mean ensures that each prediction time domain has equal weight in the overall evaluation, preventing any single time domain's index from dominating the overall evaluation result. This allows for a comprehensive reflection of the network's overall prediction accuracy in the short, medium, and long term. If the network predicts very accurately on day 1 but has a large deviation on day 14, the overall evaluation index can simultaneously reflect both short-term advantages and long-term shortcomings. If the network performs stably across all time domains, the overall evaluation index can also provide a correspondingly high score.

[0240] Specifically, the arithmetic mean of the macro average root mean square error corresponding to the third preset number of prediction time domains is taken as the macro average root mean square error in the overall evaluation index, and the arithmetic mean of the macro average coefficient of determination corresponding to the third preset number of prediction time domains is taken as the macro average coefficient of determination in the overall evaluation index. The prediction accuracy of the spatiotemporal prediction network is evaluated by using the macro average root mean square error and the macro average coefficient of determination in the overall evaluation index together. The smaller the macro average root mean square error and the closer the macro average coefficient of determination is to 1, the higher the prediction accuracy of the network.

[0241] To verify the effectiveness of the method of the present invention, experimental verification was conducted using a real-world dataset, as detailed below.

[0242] I. Experimental Data and Setup

[0243] This experiment used daily water level observation records from three groundwater monitoring wells in the target area for approximately three years (July 2023 to July 2026) as experimental data. The data came from publicly available groundwater monitoring stations (USGS), and the original records were water level depths in feet. The elevations of the measurement points (ground surface) relative to the unified elevation datum (North American Vertical Datum 1988, NAVD88) obtained from the three monitoring wells were 111.1 ft, 114.9 ft, and 153.7 ft, respectively. After conversion to meters using a conversion factor of 0.3048, these are approximately 33.86 m, 35.02 m, and 46.85 m, respectively. In this embodiment, the water level data was converted into groundwater head values ​​based on the unified elevation datum according to the formula "water head = measurement point elevation − water level depth", and this value was used as a modeling and prediction variable. The water head values ​​of the three monitoring wells ranged from approximately 29.70 to 31.31 m, 32.51 to 33.53 m, and 43.28 to 44.61 m, respectively. The differences in water head between the different monitoring wells have physical meaning under the unified elevation datum. After preprocessing the raw water level data, including missing value removal and imputation, unit conversion, daily aggregation, and multi-well merging, a complete multi-well daily head sequence of 1096 days was obtained. According to the chronological order of monitoring time, the original monitoring time interval was divided into a training time interval (767 days), a validation time interval (164 days), and a test time interval (165 days) in a ratio of 7:1.5:1.5. Time series samples were constructed in each time interval, resulting in 740 training samples, 151 validation samples, and 152 test samples. The monitoring times covered by the prediction targets of the samples in different time intervals do not overlap. The normalization parameters, inter-well head difference benchmark, prior reference scale, and upper limit of change rate are determined only based on the data in the training time interval, and the test samples are fixed for evaluation.

[0244] The spatiotemporal prediction network adopts the structure provided in this embodiment. The first preset quantity is 14, and the second preset quantity is also 14. This means that it inputs water level data from multiple wells for 14 consecutive days and outputs a water level prediction sequence for each monitoring well for the next 14 consecutive days. The multi-step output layer uses the known water level at the last historical monitoring moment of each monitoring well as a reference, outputting the water level change for each monitoring well for the next 14 consecutive days. The known water level and the water level change are then superimposed to obtain the water level prediction sequence. In the spatial feature extraction layer, the K-nearest neighbor parameter K=2. With 3 monitoring wells, each well is connected to the other two. The non-zero elements of the graph adjacency matrix are used to limit the graph attention aggregation range. The aggregation weight within the neighborhood is determined by the monitoring well features through an attention mechanism. The system is adaptively determined, and the numerical distance weights in the adjacency matrix are not directly added to the attention score. The hidden feature dimension is 48, the random inactivation rate is 0.15, and the internal weights of the three statistical prior constraints are 0.40, 0.30, and 0.30, respectively. They are then weighted after being divided by a reference scale determined within the training time interval. The target value of the prior constraint weights is selected as 0.7 from the candidate set {0.05, 0.1, 0.2, 0.3, 0.5, 0.7} based on the root mean square error of the 14-day mean of the validation set. The number of warm-up rounds is 10. The maximum number of iteration rounds in the pre-training stage is 180, and the maximum number of iteration rounds in the fine-tuning stage is 80. Both stages adopt an early stopping mechanism based on validation loss improvement.

[0245] During training, in the pre-training phase, training samples are input into the network in batches of 32, with an initial learning rate of 1e-3, a weight decay coefficient of 1e-4, a maximum number of iterations of 180, an early stopping patience period of 35, and a minimum number of iterations of 55. In the constraint fine-tuning phase, samples are also input in batches of 32, with a learning rate of 3e-4, a maximum number of iterations of 80, an early stopping patience period of 25, and a minimum number of iterations of 15. The optimal network parameters are determined based on the joint loss of the data prediction loss on the validation samples and the prior constraint residuals. In both phases, the gradients of the network parameters are clipped, with a gradient clipping threshold of 5.0 to prevent excessively large gradient magnitudes from causing training instability.

[0246] To objectively evaluate different network structures and training schemes, seven methods were compared under the same data partition: the complete model (PINN-RGNN-ST, i.e., graph attention aggregation, recurrent network, incremental multi-step output, three statistical prior constraints, and fine-tuning after pre-training), RGNN-ST (same prediction network, single-stage training with data prediction loss), MV-LSTM (multivariate recurrent network with direct multi-step output); PINN-MV-LSTM (multivariate recurrent network with three statistical prior constraints and two-stage training); GCN-Direct (graph convolutional features and direct multi-step output head, without recurrent network), PINN-RGNN-ST-nopre (no pre-training parameters, direct training with joint loss), and the persistent baseline (using the water level at the last known monitoring time as the predicted value for the next 14 consecutive days). Except for the persistent baseline, all methods used the same input length, data partition, and evaluation metrics, trained independently with five random seeds, and reported the mean ± standard deviation of the root mean square error and the coefficient of determination.

[0247] The root mean square error (RMSE) is the square root of the mean square of the difference between the predicted and actual values, measured in meters (m). A smaller RMSE indicates higher prediction accuracy. The coefficient of determination (COD) measures the degree to which the prediction results explain the changes in the actual values; a COD closer to 1 indicates better prediction performance. The RMSE and COD for each prediction time domain are calculated for each monitoring well, and then a macro-average is taken for all monitoring wells. The 14-day average index is the arithmetic mean of the macro-average indices corresponding to the 1st to 14th prediction time domains. In addition, individual indices for each monitoring well and the overall RMSE calculated by combining all monitoring wells, all test samples, and all prediction time domains are supplemented. The prediction time domains are divided into three groups: 1-3 days, 4-7 days, and 8-14 days. The arithmetic mean of the macro-average indices for each group is reported to examine the effect of prior constraints on error accumulation in the medium to long time domains.

[0248]

[0249]

[0250] in, The root mean square error, For the i-th water level, Let be the predicted water level for the i-th sample, and n be the total number of samples. R is the average of the true water level values. 2 The coefficient of determination.

[0251] II. Overall Prediction Accuracy Comparison

[0252] Under 767 days of training data, the prediction accuracy (mean ± standard deviation of 5 random seeds) of each method on the 14-day mean index is compared as follows: The root mean square error (RMSE) of the complete model PINN-RGNN-ST is 0.0804 ± 0.0015m, and the coefficient of determination (COP) is 0.7803 ± 0.0080; RGNN-ST is 0.0879 ± 0.0020m, and the COP is 0.7365 ± 0.0117; MV-LSTM is 0.0749 ± 0.0044m, and the COP is 0.7365 ± 0.0117. The root mean square error (RMSE) of the complete model is 0.8144±0.0226; for GCN-Direct it is 0.0831±0.0003m, and the RMSE is 0.7653±0.0018; for PINN-RGNN-ST-nopre it is 0.0888±0.0036m, and the RMSE is 0.7428±0.0129; for PINN-MV-LSTM it is 0.0767±0.0069m, and the RMSE is 0.8046±0.0369; and for the persistent baseline it is 0.0839m, with a RMSE of 0.7605. The complete model reduces the RMSE by approximately 4.3% compared to the persistent baseline and by approximately 9.5% compared to PINN-RGNN-ST-nopre. Compared to the single-stage trained RGNN-ST, the complete model has lower RMSEs in all comparisons with the same seed across 5 random seeds, with an average difference of approximately 8.6%. The above comparison reflects the overall differences between the complete two-stage constrained training scheme and each control scheme. Since RGNN-ST did not undergo the same additional fine-tuning stage as the complete model, the differences cannot be entirely attributed to prior constraints based solely on this set of comparisons. MV-LSTM achieved the best 14-day average accuracy in this embodiment, indicating that under the condition of only 3 monitoring wells and relatively stable sequences, multivariate recurrent networks can directly utilize multi-well information. This embodiment does not advocate that the graph module brings accuracy gains under all data conditions.

[0253] III. Comparison of Prediction Time Domain Accuracy

[0254] Under 767 days of training data, the prediction accuracy (mean of 5 random seeds) of each method in the prediction time domain at 1 day, 3 days, 7 days, and 14 days is as follows: Figure 2As shown. The root mean square errors (RMSEs) of the complete model in the four prediction time domains are 0.0161 m, 0.0388 m, 0.0787 m, and 0.1357 m, respectively, with coefficients of determination of 0.9914, 0.9568, 0.8272, and 0.4814, respectively. These RMSEs are lower than those of the persistent baseline (0.0163 m, 0.0399 m, 0.0816 m, 0.1434 m) and the single-stage trained RGNN-ST (0.0166 m, 0.0411 m, 0.0854 m, 0.1502 m). In the 14-day prediction time domain, the RMSE of MV-LSTM is 0.1201 m, lower than that of the complete model. The above results are used to illustrate the actual performance of the seven methods in this embodiment on a fixed test set; the difference between the complete model and RGNN-ST includes the effects of both two-stage training and prior constraints, and is not considered as causal proof of the independent contribution of prior constraints.

[0255] IV. Two-stage scheme and statistical prior indicator analysis

[0256] Under 767 days of training data, the 14-day root mean square error (RMSE) of the complete model PINN-RGNN-ST is 0.0804 m, lower than 0.0879 m for the single-stage RGNN-ST and 0.0888 m for PINN-RGNN-ST-nopre without pre-trained parameters. The above results only illustrate the overall performance of the complete two-stage constrained training scheme in this embodiment and cannot separate the accuracy difference into the independent contributions of prior constraints and additional fine-tuning. Statistical prior terms on the test set show that the average values ​​of the five random seeds for the complete model and RGNN-ST are: water level change consistency (0.167692 vs. 0.121636), inter-well head difference stability (1.334490 vs. 1.397537), and time series smoothing (0.002182 vs. 0.004597). The complete model has lower values ​​for the latter two terms and higher values ​​for the first term. Therefore, this embodiment only objectively reports the differences in three indicators and does not advocate for improvements in all three statistical prior indicators.

[0257] V. Analysis of Uncertainty in Training Data Volume

[0258] To examine the changes in results under different training history lengths, seven methods were trained using historical data from approximately 50, 100, 200, 300, 500, and 767 days prior to the validation period, with the validation and testing periods remaining constant. The root mean square error of each method on the fourteen-day mean index was compared. Figure 3As shown, with approximately 50, 100, and 200 days of training data, the results of each deep learning method fluctuated significantly, with the persistence baseline (0.0839 m) being competitive. With approximately 300 days of training data, the root mean square error of the complete model was 0.0776 m, lower than the persistence baseline (0.0839 m) and RGNN-ST (0.0828 m). With approximately 500 days of training data, the complete model had an error of 0.0806 m, lower than RGNN-ST (0.0829 m). With approximately 767 days of training data, the complete model had an error of 0.0780 m, lower than the persistence baseline (0.0839 m) and RGNN-ST (0.0863 m). The above results only indicate that, under the fixed data partitioning and random training conditions of this embodiment, the complete model achieved lower errors in several settings with approximately 300 days or more of training history, and should not be used to infer the minimum number of training days required for other regions or other datasets.

[0259] VI. Training Process Observation

[0260] Under 767 days of training data, the records show that the constraint fine-tuning phase of the full model achieved the lowest joint validation loss on average in about one iteration, while RGNN-ST achieved the lowest validation error on average in about 35 iterations during single-stage training. Since the fine-tuning phase of the full model starts from pre-trained parameters, and the model selection metrics, prior weight warm-up, and early stopping settings differ from those in single-stage training, this difference in iterations is only recorded as part of the training process and does not prove that constraint fine-tuning has an independent fast or stable convergence effect.

[0261] VII. Summary of Example Results. Based on the above example results, the method of this invention uses groundwater head under a unified elevation benchmark as the modeling variable. Within the adjacency range defined by the K-nearest neighbor graph, it determines the aggregation weights based on the characteristics of the monitoring wells through an attention mechanism, extracts temporal features using a recurrent network, and outputs the future water level change based on the last known water level. Simultaneously, it constructs three statistical prior constraints—water level change consistency, inter-well head difference stability, and temporal smoothing—from the monitoring data during the training period, employing a two-stage training approach of pre-training followed by joint fine-tuning. The example results show that the complete model's 14-day mean error on a fixed test set is lower than that of the persistent baseline, single-stage RGNN-ST, and PINN-RGNN-ST-nopre without pre-training parameters, but higher than MV-LSTM and PINN-MV-LSTM. Among the three statistical prior indices, the complete model's inter-well head difference stability and temporal smoothing terms are lower than RGNN-ST, while its water level change consistency term is higher than RGNN-ST. The above results support the idea that this method is an feasible joint training scheme for spatiotemporal prediction and statistical priors, but we do not interpret the overall differences in different training stages as independent accuracy contributions of prior constraints, nor do we claim that this method is superior to all comparative methods under all data conditions.

[0262] See Figure 4 This is a schematic diagram of a groundwater level prediction system based on graph neural networks and hydrogeological prior constraints provided in an embodiment of the present invention. The system includes:

[0263] The processing module is used to acquire groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area, and to preprocess the groundwater level data to obtain multi-well spatiotemporal sequence data.

[0264] The construction module is used to construct a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data, and to construct a time-series input sample based on the multi-well spatiotemporal sequence data;

[0265] The prediction module is used to construct a spatiotemporal prediction network. The time-series input samples are input into the spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future. The spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer.

[0266] The constraint module is used to construct hydrogeological prior constraint residuals based on the water level prediction sequence. The hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term.

[0267] The training module is used to construct a joint loss function based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, as well as the hydrogeological prior constraint residuals. Based on the joint loss function, the spatiotemporal prediction network is fine-tuned and trained to obtain the trained spatiotemporal prediction network.

[0268] The output module is used to input the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network to obtain the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times.

[0269] The system can be used to execute accordingly. Figure 1 The steps in the method embodiments shown are implemented in a similar manner and have similar technical effects, and will not be repeated here.

[0270] In a third aspect, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the above-described groundwater level prediction method based on graph neural networks and hydrogeological prior constraints. The implementation principle and technical effects are similar to those of the aforementioned method embodiments, and will not be repeated here.

[0271] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A groundwater level prediction method based on graph neural networks and hydrogeological prior constraints, characterized in that, include: The groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area are acquired, and the groundwater level data is preprocessed to obtain multi-well spatiotemporal sequence data. A graph adjacency matrix is ​​constructed between each of the groundwater monitoring wells based on the spatial coordinate data, and a time-series input sample is constructed based on the multi-well spatiotemporal sequence data. A spatiotemporal prediction network is constructed. The time-series input samples are input into the spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future. The spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer. Hydrogeological prior constraint residuals are constructed based on the water level prediction sequence, and the hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term; Based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, and the hydrogeological prior constraint residual, a joint loss function is constructed. The spatiotemporal prediction network is then fine-tuned and trained based on the joint loss function to obtain the trained spatiotemporal prediction network. By inputting the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network, the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times are obtained.

2. The method according to claim 1, characterized in that, The preprocessing of the groundwater level data to obtain multi-well spatiotemporal sequence data includes: Missing values ​​were removed and imputed from the groundwater level data of groundwater monitoring wells in various locations. The water level data of each monitoring well were converted into groundwater head values ​​or directly used as modeling data based on a unified elevation benchmark. The water level data of each monitoring well are aggregated according to a preset time granularity to obtain the time series data of each monitoring well. The time series data of each monitoring well are then intersected and merged according to the monitoring time. The missing positions in the merged sequence are then interpolated to obtain a multi-well spatiotemporal sequence matrix. Based on the water level data within the training time interval in the multi-well spatiotemporal sequence matrix, a normalization parameter is determined. The water level data in the multi-well spatiotemporal sequence matrix is ​​then mapped to a preset numerical range according to the normalization parameter to obtain multi-well spatiotemporal sequence data.

3. The method according to claim 1, characterized in that, The step of constructing a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data includes: The latitude and longitude coordinates of each groundwater monitoring well are projected onto a plane rectangular coordinate system, and the spatial distance between any two groundwater monitoring wells is calculated based on the projected plane rectangular coordinate system. Based on the spatial distances described, the K-nearest neighbor method is used to determine the neighboring monitoring wells of each groundwater monitoring well, an initial adjacency matrix is ​​constructed, and the feature distance is determined based on the statistical characteristics of all spatial distances. Based on the spatial distance and the feature distance, the adjacency weight between each groundwater monitoring well and its neighboring monitoring wells is determined. Each adjacency weight is filled into the corresponding adjacency position in the initial adjacency matrix. A self-loop is added to the initial adjacency matrix. The adjacency matrix after adding the self-loop is then symmetricized and normalized to obtain the graph adjacency matrix.

4. The method according to claim 1, characterized in that, The construction of time-series input samples based on the multi-well spatiotemporal sequence data includes: From the multi-well spatiotemporal sequence data, extract the multi-well water level data corresponding to a first preset number of consecutive monitoring times as a sub-input sample, and determine the multi-well water level data corresponding to a second preset number of consecutive monitoring times after the first preset number of monitoring times as the multi-well water level prediction target corresponding to the sub-input sample. According to the order of monitoring time, the original monitoring time intervals corresponding to the multi-well spatiotemporal sequence data are divided into training time intervals, verification time intervals, and testing time intervals. Sub-input samples are constructed within the training time interval, verification time interval, and testing time interval respectively. Following the chronological order of the monitoring times of the multi-well spatiotemporal sequence data, one monitoring time is shifted backward each time. The operations of extracting the sub-input samples and determining the multi-well water level prediction targets are repeated to obtain multiple sets of sub-input samples and their corresponding multi-well water level prediction targets. The time-series input sample is composed of multiple sets of corresponding sub-input samples and multiple well water level prediction targets.

5. The method according to claim 1, characterized in that, The construction of the spatiotemporal prediction network involves inputting the time-series input samples into the spatiotemporal prediction network to obtain a water level prediction sequence for each monitoring well at multiple consecutive monitoring times in the future, including: Based on the spatial feature extraction layer, the graph adjacency matrix, and the attention mechanism, the aggregation weight is determined, and the water level data of each monitoring well corresponding to each monitoring time are aggregated in the neighborhood to obtain multi-well spatial features. Based on the time feature extraction layer, the spatial features of the multi-wells corresponding to each monitoring well are time-series modeled according to the order of monitoring time to obtain the time-series feature information of each monitoring well; Based on the multi-step output layer and the temporal characteristic information of each monitoring well, the water level prediction sequence of each monitoring well is output for a second preset number of consecutive monitoring times in the future.

6. The method according to claim 5, characterized in that, The aggregation weights are determined based on the spatial feature extraction layer, the graph adjacency matrix, and the attention mechanism. Neighborhood aggregation is then performed on the water level data of each monitoring well corresponding to each monitoring time to obtain multi-well spatial features, including: The water level data of each monitoring well at each monitoring time are sequentially linearly transformed, and the aggregation weight is determined based on the graph adjacency matrix and attention mechanism. Then, neighborhood feature aggregation, residual connection and nonlinear activation processing are performed to obtain the multi-well spatial features corresponding to each monitoring time. The time feature extraction layer performs time-series modeling on the spatial features of the multi-wells corresponding to each monitoring well according to the chronological order of monitoring time, to obtain the time-series feature information of each monitoring well, including: The spatial characteristics of each monitoring well at a first preset number of consecutive monitoring times are input into a loop network in chronological order of monitoring times. Based on the loop network, the temporal characteristic information of each monitoring well at the last historical monitoring time is output. The step of outputting a water level prediction sequence for each monitoring well based on the multi-step output layer and the temporal characteristic information of each monitoring well for a second preset number of consecutive monitoring times includes: Input the time-series characteristic information of each monitoring well at the last historical monitoring time into the linear mapping unit, and directly output the water level change of each monitoring well in the next second preset number of monitoring times; Using the known water level corresponding to the last historical monitoring moment in the time-series input sample as a benchmark, the known water level is superimposed with the water level change to obtain the water level prediction sequence for each monitoring well in the next second preset number of monitoring moments.

7. The method according to claim 1, characterized in that, The hydrogeological prior constraint residuals are constructed based on the water level prediction sequence. These residuals include a first constraint term, a second constraint term, and a third constraint term, comprising: The water level prediction sequence is concatenated with the water level data corresponding to the last historical monitoring moment in the time series input sample to obtain the continuous water level prediction sequence of each monitoring well, and the continuous water level prediction sequence is restored to the water level sequence with physical dimensions. Based on the water level difference between adjacent monitoring times and the water level range of each monitoring well within the training time interval, the normalized water level change of each monitoring well is determined. Based on the deviation between the normalized water level change of each monitoring well at each monitoring time and the average normalized water level change of all monitoring wells at that monitoring time, and normalized by the overall change range of the normalized water level change of each monitoring well, the first constraint term is constructed. Based on the difference in the average water level of each monitoring well within the training time interval, a benchmark for the relative difference between each pair of monitoring wells is determined. Based on the deviation between the water level difference of each pair of monitoring wells in the continuous water level prediction sequence and the benchmark for the relative difference between the pairs of monitoring wells, the second constraint term is constructed. Based on the difference results of the values ​​of adjacent monitoring times of each monitoring well in the continuous water level prediction sequence, the third constraint term is constructed. Based on the water level data within the training time interval, the reference scales corresponding to the first constraint term, the second constraint term, and the third constraint term are determined. The first constraint term, the second constraint term, and the third constraint term are divided by their respective reference scales. The first constraint term, the second constraint term, and the third constraint term after reference scale normalization are weighted and summed according to preset constraint weights to obtain the hydrogeological a priori constraint residuals.

8. The method according to claim 4, characterized in that, The process involves constructing a joint loss function based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, and the hydrogeological prior constraint residuals. The spatiotemporal prediction network is then subjected to constraint fine-tuning training based on this joint loss function to obtain the trained spatiotemporal prediction network, including: Based on the mean square error between the water level prediction sequence and the multi-well water level prediction target, the data prediction loss is obtained. The data prediction loss is then summed with the product of the preset prior constraint weights and the hydrogeological prior constraint residuals to obtain the joint loss function. The spatiotemporal prediction network is iteratively trained based on the data prediction loss until the preset pre-training convergence condition is met, and the pre-trained network parameters are obtained. The network parameters of the spatiotemporal prediction network are initialized to the pre-trained network parameters, and during the training iteration process, the prior constraint weights are gradually increased from the initial weight values ​​to the preset target weight values ​​according to the preset number of warm-up rounds. The network parameters of the spatiotemporal prediction network are iteratively updated based on the joint loss function until the preset fine-tuning convergence condition is met, thus obtaining the trained spatiotemporal prediction network.

9. The method according to claim 1, characterized in that, After obtaining the predicted groundwater levels for each monitoring well at multiple consecutive monitoring times, the method further includes: The groundwater level prediction results of each monitoring well are converted into water level values ​​with physical dimensions, and the water level prediction values ​​of each monitoring well in each prediction time domain are determined according to multiple preset prediction time domains. Based on the predicted water level and the actual water level of each monitoring well in each prediction time domain, the root mean square error and coefficient of determination of each monitoring well are calculated, and the root mean square error and coefficient of determination of all monitoring wells are taken as macro averages to obtain the macro average index corresponding to each prediction time domain. The arithmetic mean of the macro-average indices corresponding to the third preset number of prediction time domains is used as the overall evaluation index for multi-step prediction, in order to evaluate the prediction accuracy of the spatiotemporal prediction network.

10. A groundwater level prediction system based on graph neural networks and hydrogeological prior constraints, characterized in that, include: The processing module is used to acquire groundwater level data and spatial coordinate data of groundwater monitoring wells in the target area, and to preprocess the groundwater level data to obtain multi-well spatiotemporal sequence data. The construction module is used to construct a graph adjacency matrix between each of the groundwater monitoring wells based on the spatial coordinate data, and to construct a time-series input sample based on the multi-well spatiotemporal sequence data; The prediction module is used to construct a spatiotemporal prediction network. The time-series input samples are input into the spatiotemporal prediction network to obtain the water level prediction sequence of each monitoring well at multiple consecutive monitoring times in the future. The spatiotemporal prediction network includes a spatial feature extraction layer, a temporal feature extraction layer, and a multi-step output layer. The constraint module is used to construct hydrogeological prior constraint residuals based on the water level prediction sequence. The hydrogeological prior constraint residuals include a first constraint term, a second constraint term, and a third constraint term. The training module is used to construct a joint loss function based on the data prediction loss between the water level prediction sequence and the multi-well water level prediction target, as well as the hydrogeological prior constraint residuals. Based on the joint loss function, the spatiotemporal prediction network is fine-tuned and trained to obtain the trained spatiotemporal prediction network. The output module is used to input the groundwater level data of each monitoring well in the target area at multiple consecutive historical monitoring times into the trained spatiotemporal prediction network to obtain the groundwater level prediction results of each monitoring well at multiple consecutive future monitoring times.

Images