A groundwater level simulation method based on lag response characteristics of multi-aquifer system
By combining cross-wavelet transform and multi-head self-attention mechanism network with gated cyclic units, the problems of identifying the main control factors and quantifying the hysteresis response relationship in multi-aquifer systems are solved, achieving efficient and accurate groundwater level simulation and improving the physical interpretability and simulation accuracy of the model.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- CAPITAL NORMAL UNIVERSITY
- Filing Date
- 2025-10-28
- Publication Date
- 2026-04-17
AI Technical Summary
Existing GRU models and improved attention mechanism models suffer from insufficient ability to identify controlling factors, lack of quantification of interlayer hysteresis response relationships, and low simulation efficiency in multi-aquifer systems, resulting in limited simulation accuracy and poor physical interpretability.
The hysteresis response features are extracted using the cross-wavelet transform method. Combined with multi-head self-attention mechanism and gated recurrent unit, the input features are linearly mapped and weighted through the multi-head self-attention mechanism network. Combined with GRU for time series simulation, a groundwater level simulation model is constructed.
It enables efficient and accurate simulation of multi-aquifer systems, improves the physical interpretability and simulation accuracy of the model, simplifies the model construction process, and reduces complexity.
Smart Images
Figure CN121389767B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of groundwater level simulation technology, and in particular to a groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system. Background Technology
[0002] Groundwater serves as a core freshwater source for human daily life, agricultural irrigation, and industrial production. The temporal and spatial variations in groundwater levels directly impact regional water resource security and ecological stability. Accurately simulating groundwater level dynamics is a crucial prerequisite for implementing water use regulation, groundwater over-extraction control, land subsidence prevention, and emergency drought response. It can provide core technological support for the sustainable and efficient utilization of groundwater resources.
[0003] Current groundwater level simulation methods are mainly divided into two categories: numerical models based on physical mechanisms and data-driven machine learning models. Traditional numerical simulation methods rely on clear physical mechanisms of groundwater flow, using numerical techniques such as the finite element method and finite difference method to discretize the generalized geological structure, and then solve the groundwater flow control equations to simulate the groundwater flow process. This approach can achieve high-precision simulations in areas with detailed geological data. Machine learning models, on the other hand, are based on a data-driven approach, starting from the time-series data itself to uncover the implicit nonlinear relationships between groundwater level and influencing factors (such as precipitation and extraction volume). They have advantages such as fast computation speed and strong nonlinear mapping capabilities, and can, to some extent, compensate for the high dependence of numerical models on geological data.
[0004] Gated recurrent unit (GRU) models, in particular, possess unique advantages in processing long-term time-series data and have shown significant application value in areas with scarce hydrogeological data. However, existing GRU models and conventional attention mechanism models based on them (such as Attention-GRU) still have certain shortcomings in simulating groundwater levels in multi-aquifer systems:
[0005] Insufficient ability to identify key controlling factors: The GRU model treats all input variables (such as precipitation, extraction volume, and overlying aquifer water level) with equal weights, which cannot effectively distinguish the degree of influence of different variables on groundwater level. This results in the model being insensitive to key influencing factors and limiting simulation accuracy.
[0006] The interlayer hysteresis response relationship is not quantified: In a multi-aquifer system, unconfined aquifers are mainly recharged by atmospheric precipitation, while confined aquifers rely on the overflow recharge of overlying aquifers. There are obvious water level hysteresis response characteristics between different aquifers. However, existing models do not quantitatively calculate and incorporate this hysteresis characteristic, which cannot truly reflect the hydrogeological processes of multi-aquifer systems and results in poor physical interpretability of the models.
[0007] Low efficiency in multi-aquifer simulation: Existing models often require the construction of independent simulation models for unconfined aquifers and each confined aquifer, which not only increases the complexity of model construction, but may also lead to the loss of correlation information between different aquifers, making it difficult to achieve efficient and accurate simulation of water levels in multiple aquifers simultaneously. Summary of the Invention
[0008] The purpose of this invention is to solve the problems mentioned in the background art by proposing a groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system.
[0009] To achieve the above objectives, the present invention adopts the following technical solution:
[0010] A groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system includes the following steps:
[0011] Step 1: Extraction of hysteresis response features: Based on the cross wavelet transform method, calculate the hysteresis time between the unconfined aquifer level and precipitation, as well as the hysteresis time between the confined aquifer level and the overlying aquifer level.
[0012] Step 2: Model Input Data Preparation: Prepare an input dataset containing hysteresis information for the target simulated aquifer.
[0013] Step 3, Multi-head Self-Attention Mechanism Feature Recognition: Input the input dataset prepared in step S2 into the multi-head self-attention mechanism network. The multi-head self-attention mechanism network performs linear mapping on the input features to obtain the query matrix, key matrix and value matrix, calculates the self-attention of each attention head, and concatenates and linearly transforms the outputs of all attention heads. Then, it enhances the feature expression through residual connection and nonlinear activation function to obtain the feature-enhanced output state.
[0014] Step 4: Timing simulation of the gated loop unit: The output state of the feature enhancement obtained in step S3 is used as input to the gated loop unit for timing simulation. The gated loop unit controls the retention and updating of historical information through update gate and reset gate, and outputs the hidden state.
[0015] Step 5, Output Results: The hidden state output by the gated loop unit is converted into a time series format by the decoder to obtain the groundwater level simulation results of the target aquifer.
[0016] Preferably, when preparing the input dataset containing lag information in step two, the specific steps include:
[0017] For unconfined aquifers, the input dataset includes: precipitation, groundwater extraction from the aquifer, groundwater level of the aquifer in the previous month, and precipitation data adjusted based on the lag time obtained in step S1.
[0018] For confined aquifers, the input dataset includes: precipitation, groundwater extraction from the aquifer, the groundwater level of the aquifer in the previous month, and groundwater level data of the overlying aquifer adjusted based on the lag time obtained in step S1.
[0019] Furthermore, in step one, Morlet wavelet is used as the wavelet function of the cross wavelet transform. The lag time is determined by calculating the product of the wavelet transform coefficients of the two time series and their complex conjugates, and based on the obtained phase difference.
[0020] Furthermore, in step three, the operation of the multi-head self-attention mechanism network includes:
[0021] Perform linear mapping on the input features E to generate the query matrix. Key matrix Sum matrix The specific formula is as follows:
[0022]
[0023] in To focus on the number of heads, , and These are the corresponding linear mapping matrices;
[0024] For each attention head, calculate the self-attention output using the following formula:
[0025]
[0026] in The dimension of a single attention head;
[0027] The outputs of all attention heads are concatenated, and a linear transformation is applied to obtain the multi-head self-attention output, as shown in the formula:
[0028]
[0029] in For the output of each attention head, It is a linear transformation matrix;
[0030] The multi-head self-attention output is added to the original input feature E through a residual network and then processed by a nonlinear activation function to obtain the final feature-enhanced output state, as shown in the formula:
[0031]
[0032] in It is a non-linear activation function.
[0033] Preferably, in step four, the operation of the gated loop unit includes:
[0034] At each time step t, receive the current input. The hidden state of the previous time step ; Calculate the reset gate and the update gate ;
[0035]
[0036]
[0037] in It's the sigmoid function, which transforms data to a value in the range of 0-1, serving as a gating signal. and To reset the gate weight, and To update the gate weights;
[0038] Calculate candidate hidden state The formula is:
[0039]
[0040] in For activation function, , As weight;
[0041] Calculate the hidden state at the current time step. The formula is:
[0042] .
[0043] Furthermore, the precipitation data in step one is monthly data with a spatial resolution of 1 kilometer, the groundwater level data is monthly monitoring data of unconfined aquifer, second confined aquifer, and third confined aquifer, and the groundwater extraction data is monthly data of the first to third aquifers in different districts and counties.
[0044] Furthermore, in step four, when constructing the coupled model, the dataset is divided into a training set and a validation set. 70% of the data is used as the training set, 30% of the data is used as the validation set, and the data from the last 12 months in the time series is used as the validation set, with the remainder used as the training set.
[0045] Compared with existing technologies, this invention provides a groundwater level simulation method based on the hysteresis response characteristics of multi-aquifer systems, which has the following beneficial effects:
[0046] 1. This invention uses Morlet wavelet transform with cross-wavelet transform as the basis function to quantitatively calculate the lag time between unconfined aquifers and precipitation, and between confined aquifers and overlying aquifers. This lag characteristic is then incorporated into the model input and training process, enabling the model to accurately reflect the replenishment lag relationship of "precipitation-unconfined water-confined water" in multi-aquifer systems. This fills the gap in existing machine learning models that lack physical mechanism support and effectively improves the physical interpretability of the model.
[0047] 2. This invention introduces a multi-head self-attention mechanism and sets multiple attention heads to linearly map and weight the input features, which can automatically identify the influence weight of different input variables on the groundwater level, thus avoiding the defect of "equalization of variable weights" in the traditional GRU model.
[0048] 3. This invention utilizes MSAM to assign differentiated weights to the input features of each aquifer, combined with GRU's efficient processing capabilities for long-term time-series data, to simultaneously simulate the water levels of unconfined aquifers, second confined aquifers, and third confined aquifers. Compared to the traditional "one aquifer, one model" approach, this significantly simplifies the model construction process, reduces repetitive work in data preprocessing and model training, and while ensuring simulation accuracy, significantly improves the efficiency of groundwater level simulation in multi-aquifer systems, while reducing the complexity of model application. Attached Figure Description
[0049] Figure 1 This is a flowchart of the multi-head attention mechanism in a groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system proposed in this invention.
[0050] Figure 2 This is a flowchart of the gated loop unit structure in a groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system proposed in this invention.
[0051] Figure 3 This is a simulation flowchart of a groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system proposed in this invention. Detailed Implementation
[0052] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.
[0053] Example 1:
[0054] Reference Figures 1-3Taking the multi-aquifer system of the Huairou emergency water source area as the research object, this paper details the invention "A groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system". This method can accurately simulate the groundwater level changes of the unconfined aquifer (first aquifer), the second confined aquifer, and the third confined aquifer in the Huairou emergency water source area. The specific implementation steps are as follows:
[0055] The basic data required for this embodiment are all based on actual monitoring and publicly available data from the Huairou emergency water source area, specifically including:
[0056] Precipitation data: Monthly precipitation data, sourced from the National Tibetan Plateau Scientific Data Center;
[0057] Groundwater level data: Monthly groundwater level monitoring well data of Huairou emergency water source area, covering unconfined aquifer, second confined aquifer and third confined aquifer, with a monitoring frequency of once a month, and the time span is consistent with precipitation data;
[0058] Groundwater extraction data: Monthly extraction data of the first to third aquifers in the Huairou emergency water source area.
[0059] Cross-wavelet transform for calculating lag time:
[0060] Using the cross-wavelet transform method with Morlet wavelet as the wavelet function, the two types of lag times are calculated respectively:
[0061] Lag time between unconfined aquifer (first aquifer) water level and precipitation: The time series of the monthly average water level of the unconfined aquifer is denoted as... The time series of average monthly precipitation during the same period is denoted as , Calculate the cross-wavelet transform coefficients, where for wavelet transform coefficients, for The complex conjugate of the wavelet transform coefficients is used to determine the lag time between the water level of the unconfined aquifer and precipitation based on the obtained phase difference. The calculation results show that the lag time is 1-2 months.
[0062] The lag time between the water level of the confined aquifer and the water level of its overlying aquifer: The monthly average water level time series of the second confined aquifer is denoted as... The monthly average water level time series of the first aquifer (overlying aquifer) is denoted as Calculate the cross-wavelet transform coefficients using the same formula as above. Combining the phase difference, the lag time between the second confined aquifer and the first aquifer is found to be 2-3 months; similarly, the monthly average water level time series of the third confined aquifer is denoted as... The time series of the monthly average water level of the second aquifer is denoted as... The calculation showed that the lag time between the third confined aquifer and the second aquifer was 3-4 months.
[0063] Based on the lag time obtained in step one, input datasets containing lag information are constructed for the unconfined aquifer, the second confined aquifer, and the third confined aquifer. The specific composition of the input dataset for each aquifer is as follows:
[0064] Input dataset for unconfined aquifers: monthly precipitation ( The data includes the average monthly extraction volume of the unconfined aquifer, the groundwater level of the unconfined aquifer in the previous month, and the precipitation data adjusted based on a "1-2 month lag time" (i.e., the precipitation input at time t is the actual precipitation at time t-1 or t-2).
[0065] Input dataset for the second confined aquifer: monthly precipitation ( The data includes the average monthly extraction volume of the second confined aquifer, the groundwater level of the second confined aquifer in the previous month, and the water level data of the first aquifer adjusted based on a "2-3 month lag time" (i.e., the water level of the first aquifer input at time t is the actual water level at time t-2 or t-3).
[0066] Input dataset for the third confined aquifer: monthly precipitation ( The data includes the average monthly extraction volume of the third confined aquifer, the groundwater level of the third confined aquifer in the previous month, and the water level data of the second aquifer adjusted based on a "3-4 month lag time" (i.e., the water level of the second aquifer input at time t is the actual water level at time t-3 or t-4).
[0067] The above input dataset is preprocessed, including missing value imputation (using linear interpolation) and outlier removal (using 3D modeling). The criteria and data standardization (mapping data to the [0,1] interval) are used to obtain a standardized dataset that can be directly input into the model.
[0068] The standardized datasets of each aquifer after preprocessing in step two are used as input features E and input into a multi-head self-attention mechanism network. The specific process is as follows:
[0069] Linear mapping generates query, key, and value matrices: With the number of attention heads h=8 and the dimension of each attention head d=64, linear mapping is performed on the input feature E to generate query matrices corresponding to the 8 attention heads. Key matrix Sum matrix The mapping process satisfies the formula:
[0070]
[0071] in , and These are the linear mapping matrices for the h-th attention head;
[0072] For each attention head, calculate the self-attention output according to the formula:
[0073]
[0074] in As a scaling factor, avoid The calculation result is too large, leading to Gradient vanishing; The function is used to normalize the attention weights to the [0,1] interval, highlighting the characteristic variables that have a significant impact on the groundwater level;
[0075] The self-attention output of 8 attention heads ( - The components are then concatenated and then subjected to a linear transformation matrix. By fusing the results, we obtain multi-head self-attention outputs. It satisfies the formula:
[0076]
[0077] The multi-head self-attention output is added to the original input feature E through a residual network and then processed by a nonlinear activation function to obtain the final feature-enhanced output state. It satisfies the formula:
[0078]
[0079] The final feature-enhanced output state will be obtained. As input The input is fed into a gated recurrent unit (GRU) for timing simulation. The specific process is as follows:
[0080] Gating signals and candidate hidden states calculation:
[0081] At each time step t, the GRU receives the current input. The hidden state of the previous time step First, the reset gate is calculated using the sigmoid activation function. and the update gate :
[0082]
[0083]
[0084] in It is the sigmoid function, which maps the calculation result to the [0,1] interval and uses it as a gate signal. and To reset the gate weight, and To update the gate weights, the matrix dimensions are (hidden layer dimension × input feature dimension) and (hidden layer dimension × hidden layer dimension). In this embodiment, the hidden layer dimension is set to 128.
[0085] Based on the reset gate Hidden state of history Filter the data and then combine it with the current input. ,pass Activation function calculates candidate hidden states :
[0086]
[0087] in The function maps the calculation result to the interval [-1, 1]. , The weights are defined in the same dimensions as the weight matrices for the reset and update gates.
[0088] According to the updated gate Controlling candidate hidden state With historical hidden state The fusion ratio is used to calculate the hidden state at the current time step. :
[0089]
[0090] in, The closer it is to 1, the more the current hidden state depends on the candidate hidden state. (i.e., the current input features); The closer it is to 0, the more the current hidden state retains the historical hidden state. (i.e., historical time series information), this mechanism effectively solves the gradient vanishing or gradient exploding problem in long time series data processing;
[0091] Results and Model Validation:
[0092] The hidden state output by the gated loop unit at the last time step The input is fed into the decoder, which converts the hidden state into a time series format through a fully connected layer to obtain the monthly groundwater level simulation results of the unconfined aquifer, the second confined aquifer, and the third confined aquifer in the Huairou emergency water source area.
[0093] Approximately 70% of the data was used as the training set, and the remaining data was used as the validation set. Root mean square error (RMSE) and coefficient of determination (R²) were used as evaluation metrics for model accuracy. The model was compared and validated with the traditional GRU model and the Attention-GRU model. The results are shown in Table 1 below:
[0094]
[0095] Table 1
[0096] Based on the data in Table 1, we can see that:
[0097] For unconfined aquifers: the RMSE of the method of this invention is 0.29m, and the R² is 0.83; the RMSE of the GRU model is 1.06m, and the R² is 0.70; the RMSE of the Attention-GRU model is 0.42m, and the R² is 0.75.
[0098] Second confined aquifer: The RMSE of the method of this invention is 0.31m, and the R² is 0.98; the RMSE of the GRU model is 1.13m, and the R² is 0.86; the RMSE of the Attention-GRU model is 0.95m, and the R² is 0.96.
[0099] For the third confined aquifer: the RMSE of the method of this invention is 1.36m and the R² is 0.88; the RMSE of the GRU model is 1.40m and the R² is 0.65; and the RMSE of the Attention-GRU model is 1.37m and the R² is 0.72.
[0100] The verification results show that the method of the present invention effectively improves the accuracy and reliability of groundwater level simulation by integrating the hysteresis response characteristics of multi-aquifer systems with the multi-head self-attention mechanism, and can be effectively applied to groundwater resource regulation and management.
[0101] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.
Claims
1. A groundwater level simulation method based on a multi-aquifer system lag response characteristic, characterized by, Includes the following steps: Step 1: Extraction of hysteresis response features: Based on the cross wavelet transform method, calculate the hysteresis time between the unconfined aquifer level and precipitation, as well as the hysteresis time between the confined aquifer level and the overlying aquifer level. Step 2: Model Input Data Preparation: Prepare an input dataset containing hysteresis information for the target simulated aquifer. Step 3: Feature recognition of multi-head self-attention mechanism: Input the input dataset prepared in step 2 into the multi-head self-attention mechanism network. The multi-head self-attention mechanism network performs linear mapping on the input features to obtain the query matrix, key matrix and value matrix, calculates the self-attention of each attention head, and concatenates and linearly transforms the outputs of all attention heads. Then, it enhances the feature expression through residual connection and nonlinear activation function to obtain the feature-enhanced output state. Step 4: Timing simulation of the gated loop unit: The output state of the feature enhancement obtained in Step 3 is used as input to the gated loop unit for timing simulation. The gated loop unit controls the retention and updating of historical information through update gate and reset gate, and outputs the hidden state. Step 5, Output Results: The hidden state output by the gated loop unit is converted into a time series format by the decoder to obtain the groundwater level simulation results of the target aquifer. Step two, in preparing the input dataset containing hysteresis information, specifically includes: For unconfined aquifers, the input dataset includes: precipitation, groundwater extraction from the aquifer, the groundwater level of the aquifer in the previous month, and precipitation data adjusted based on the lag time obtained in step one. For confined aquifers, the input dataset includes: precipitation, groundwater extraction from the aquifer, groundwater level of the aquifer in the previous month, and groundwater level data of the overlying aquifer adjusted based on the lag time obtained in step one. In step one, Morlet wavelet is used as the wavelet function of the cross wavelet transform. The lag time is determined by calculating the product of the wavelet transform coefficients of the two time series and their complex conjugates, and based on the obtained phase difference.
2. The groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system according to claim 1, characterized in that, In step three, the operation of the multi-head self-attention mechanism network includes: Perform linear mapping on the input features E to generate the query matrix. Key matrix Sum matrix The specific formula is as follows: in To focus on the number of heads, , and These are the corresponding linear mapping matrices; For each attention head, calculate the self-attention output using the following formula: in The dimension of a single attention head; The outputs of all attention heads are concatenated, and a linear transformation is applied to obtain the multi-head self-attention output, as shown in the formula: in For the output of each attention head, It is a linear transformation matrix; The multi-head self-attention output is added to the original input feature E through a residual network and then processed by a nonlinear activation function to obtain the final feature-enhanced output state, as shown in the formula: in It is a non-linear activation function.
3. The groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system according to claim 1, characterized in that, In step four, the operation of the gated loop unit includes: At each time step t, receive the current input. The hidden state of the previous time step ; Calculate the reset gate and the update gate ; in It's the sigmoid function, which transforms data to a value in the range of 0-1, serving as a gating signal. and To reset the gate weight, and To update the gate weights; Calculate candidate hidden state The formula is: in For activation function, , As weight; Calculate the hidden state at the current time step. The formula is: 。 4. The groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system according to claim 1, characterized in that, The precipitation data in Step 1 are monthly data with a spatial resolution of 1 km. The groundwater level data are monthly monitoring data of unconfined aquifer, second confined aquifer, and third confined aquifer. The groundwater extraction data are monthly data of the first to third aquifers in different districts and counties.
5. The groundwater level simulation method based on the hysteresis response characteristics of a multi-aquifer system according to claim 1, characterized in that, In step four, when performing time-series simulation of the gated cyclic unit, the dataset is divided into a training set and a validation set, with 70% of the data used as the training set and 30% of the data used as the validation set.