Underground water level change prediction method and system based on machine learning

Abstract

The invention discloses an underground water level change prediction method and system based on machine learning, and relates to the technical field of hydrogeology, and the method comprises the steps: carrying out the standardization processing of multi-source time series data, obtaining a standardized multivariable time series data matrix, and constructing a supervised learning sample set; inputting the underground water level sequence in the supervised learning sample set into a physical information guided variational mode decomposition network; decomposing the underground water level sequence into K intrinsic mode component sequences and residual term sequences through a loss function of physical driving consistency constraint; and for each of the K intrinsic mode component sequences, dynamically assembling a differentiable simulator from the differentiable simplified physical simulator basic library, and carrying out cooperative training by taking approximation to each intrinsic mode component sequence and overall reconstruction of the original water level sequence as targets to obtain K completely trained assembled differentiable simulators. According to the method, a reliable visual prediction result is generated through multi-simulator collaborative deduction and uncertainty quantization.

Description

Technical Field

[0001] This invention relates to the field of hydrogeology, and in particular to a method and system for predicting groundwater level changes based on machine learning. Background Technology

[0002] The scientific management of groundwater resources is inseparable from the accurate prediction of water level changes. Groundwater prediction methods based on machine learning have become one of the mainstream technical approaches. By collecting historical water level and driving data such as meteorology and mining, a time-series dataset is constructed. By using models such as long short-term memory networks and recurrent neural networks, the complex nonlinear relationship between water level and multiple driving factors is learned. High-dimensional and nonlinear mapping patterns are directly extracted from the data without the need to pre-set complex hydrogeological parameters. Therefore, it exhibits excellent fitting performance under good data conditions and has become an effective supplement to traditional physical models.

[0003] Existing purely data-driven prediction methods have limitations, particularly in the lack of physical mechanisms in the models. Current technologies typically treat groundwater level sequences as a whole signal for end-to-end learning, failing to explicitly model and separate the independent contributions of different physical processes within the model. This results in weak interpretability. Furthermore, due to the absence of fundamental physical constraints such as mass conservation and Darcy's law, the extrapolation predictions of the model often lack physical rationality and are difficult to guarantee when the prediction scenario exceeds the historical data distribution. Existing technologies are essentially data fitting tools, and they are insufficient in terms of physical consistency of the model, interpretability of the mechanism, and robust prediction capabilities for unknown scenarios. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a machine learning-based method for predicting groundwater level changes, which solves the problems of weak interpretability and insufficient reliability of supersample extrapolation caused by the lack of physical mechanisms in existing pure data-driven models.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: In a first aspect, the present invention provides a groundwater level change prediction method based on machine learning, which includes standardizing multi-source time series data to obtain a standardized multivariate time series data matrix and constructing a supervised learning sample set. The groundwater level sequence from the supervised learning sample set is input into a physical information-guided variational mode decomposition network; The groundwater level sequence is decomposed into K intrinsic mode component sequences and residual term sequences by using a loss function with physical driving consistency constraints. For each of the K intrinsic modal component sequences, a differentiable simulator is dynamically assembled from the basic library of differentiable simplified physical simulators. The goal is to approximate each intrinsic modal component sequence and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. The predicted sequence of driving factors for future periods is input into K well-trained assembled differentiable simulators for forward state extrapolation to obtain the predicted sequence of groundwater level changes for future periods. The prediction sequence of groundwater level changes in future periods is de-standardized and uncertainty is quantified to generate a visual prediction report.

[0007] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: standardizing multi-source time-series data to obtain a standardized multivariate time-series data matrix: The time-series data of groundwater level observation, precipitation, evaporation and extraction are time-stamped to generate time-aligned multi-source time-series data, and then filled in using a linear interpolation method to generate filled multi-source time-series data. Descriptive statistical analysis was performed on the imputed multi-source time series data to generate statistical characteristics of the mean and standard deviation of each data series. The mean and standard deviation of each data sequence are used to perform Z-score standardization on the imputed multivariate time series data to generate a standardized multivariate time series data matrix.

[0008] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the construction of a supervised learning sample set includes the following steps: The historical window length is determined to define the time span of the input features, and sliding window sampling is performed on the standardized multivariate time series data matrix to generate data slices; The standardized values ​​of each variable at each time point within the window are extracted from each data slice to form a basic feature set. Statistical features of each variable within the window are then extracted to generate derived features. The basic feature set and derived features corresponding to each data slice are merged to form a feature vector that comprehensively represents the state of historical periods, thereby obtaining the prediction step size and defining the future time point of the prediction target. For each data slice corresponding to a historical period, the standardized groundwater level observation value after the end of each feature window is located from the standardized multivariate time series data matrix; The feature vectors of each time period are paired with standardized groundwater level observations to construct a supervised learning sample set.

[0009] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: inputting the groundwater level sequence from the supervised learning sample set into a physical information-guided variational mode decomposition network: The standardized groundwater level observation sequence is formed by extracting all feature vectors corresponding to the supervised learning sample set. By analyzing the correspondence between precipitation and groundwater recharge processes, as well as between extraction and groundwater discharge processes, a set of physical driving sequences corresponding to hydrological processes is defined, and a variational mode decomposition network guided by physical information is constructed.

[0010] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the groundwater level sequence is decomposed into K intrinsic mode component sequences and residual term sequences through a loss function with physical driving consistency constraints, including the following steps: The signal reconstruction loss term in the loss function based on physical-driven consistency constraints is constrained to ensure that the sum of the K intrinsic mode component sequences and the residual term sequence is consistent with the standardized groundwater level observation sequence. The energy concentration of each intrinsic mode component sequence in the frequency domain is constrained by the modal bandwidth loss term in the loss function of the physical-driven consistency constraint. The structural differences between intrinsic modal component sequences and unspecified physical driving sequences in the time-frequency distribution are quantified by quantifying the structural consistency between each intrinsic modal component sequence and the pre-specified physical driving sequence in the time-frequency distribution. The alignment loss term separates the physical orientation by maximizing the time-frequency similarity measure between the target intrinsic mode component sequence and the specified physical driving sequence, and minimizing the time-frequency similarity measure with the physical driving sequence. During training, the parameters of the physical information-guided variational mode decomposition network are updated by jointly optimizing the signal reconstruction loss term, modal bandwidth loss term, and alignment loss term through the backpropagation algorithm. Based on a variational mode decomposition network guided by trained physical information, the standardized groundwater level observation sequence is decomposed into K intrinsic mode component sequences and residual term sequences that are strongly correlated with different physical driving processes.

[0011] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: For each of the K intrinsic mode component sequences, a differentiable simulator is dynamically assembled from a differentiable simplified physical simulator library. Based on the fundamental physical laws of groundwater flow, a basic library of a differentiable simplified physical simulator is constructed. The dynamic characteristics of each intrinsic mode component sequence are analyzed, and the feature vector of each intrinsic mode component sequence is extracted. Based on a neural network architecture, a simulator selection and parameterization network is set up. The simulator selection and parameterization network outputs the selection weight vector and parameter adjustment vector corresponding to each intrinsic mode component sequence. By using a weighted combination of differentiable simplified physics simulators in the basic library of differentiable simplified physics simulators, a preliminary assembled differentiable simulator is generated. The parameter adjustment vector is injected into the initially assembled differentiable simulator to adjust the key internal parameters of the initially assembled differentiable simulator and generate the final differentiable simulator. For each of the K intrinsic mode component sequences, feature vector extraction, simulator selection and parameterized network processing, weighted combination and parameter adjustment are repeatedly performed to dynamically assemble K differentiable simulators.

[0012] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: Co-training is performed with the goal of approximating the sequences of each intrinsic modal component and reconstructing the original water level sequence as a whole, to obtain K fully trained assembled differentiable simulators. Based on the multi-task learning framework, the total loss function of collaborative training is defined and constructed as a comprehensive index that jointly measures the independent prediction accuracy of each simulator and the overall prediction consistency. Each differentiable simulator and its corresponding dominant physics driving sequence generate an eigenmode component sequence by performing state transitions and response outputs through the simplified physical laws inherent within it. The global sequence reconstruction loss measures the fidelity of the normalized groundwater level observation sequence in terms of both overall waveform and numerical value after the intrinsic modal component sequences are superimposed and combined. During the training iteration, the gradient is calculated based on the total loss function using the backpropagation algorithm, and the parameters of the simulator selection and parameterization network, as well as the parameters inside the differentiable simulator, are updated. Set the training termination condition. When the total loss function value reaches the preset number of iterations, the training process is completed, and K assembled differentiable simulators with fixed parameters are obtained.

[0013] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: inputting the predicted sequence of driving factors for future periods into K well-trained assembled differentiable simulators for forward state extrapolation to obtain the predicted sequence of groundwater level changes for future periods. Based on the acquired driving factor forecast sequence, the precipitation forecast sequence in the driving factor forecast sequence is guided to a fully trained assembled differentiable simulator that simulates the effect of rainfall on groundwater recharge. Based on the current dynamic state of memory, combined with the input precipitation forecast sequence, a water level component sequence reflecting the contribution of future rainfall replenishment is generated; The extraction volume forecast sequence in the driving factor forecast sequence is guided to a fully trained modular differentiable simulator that simulates the effect of extraction on groundwater discharge. Based on the current memorized dynamic state and combined with the input extraction volume forecast sequence, a water level component sequence reflecting the impact of future extraction and discharge is generated. After receiving a continuous sequence of driving factor forecasts, the simulator iteratively updates its internal state and generates the corresponding water level component sequence step by step. All water level component sequences generated by each simulator throughout the entire future prediction period are superimposed point by point along the time dimension with the residual water level sequence estimated by trend extrapolation, and integrated into a complete standardized groundwater level change prediction sequence for the future period.

[0014] As a preferred embodiment of the machine learning-based groundwater level change prediction method of the present invention, the method includes the following steps: performing destandardization and uncertainty quantification on the predicted groundwater level change sequence for future periods to generate a visual prediction report: Perform an inverse Z-score transformation on the standardized groundwater level change prediction sequence to restore the standardized groundwater level change prediction sequence to a groundwater level change prediction sequence with actual physical dimensions. Apply random perturbations that conform to prior knowledge to the key parameters of the fully trained assembled differentiable simulator to generate multiple sets of parameter settings; Using a fully trained assembly-type differentiable simulator with multiple sets of parameter settings, multiple independent forward state extrapolations are performed on the future driving factor prediction sequence to generate multiple groundwater level change prediction sequences, forming a groundwater level change prediction sequence set. At each future time point in the groundwater change prediction sequence set, specific quantiles of the distribution of all predicted values ​​are extracted to obtain the prediction interval for the future time point. The groundwater level change prediction sequence with actual physical dimensions is used as the median, and the prediction interval is used as the confidence band to draw a time series curve. Add legends, axis labels, and forecast step size descriptions to the time series curve to generate a visual forecast report.

[0015] Secondly, the present invention provides a groundwater level change prediction system based on machine learning, including a construction module for standardizing multi-source time series data to obtain a standardized multivariate time series data matrix and constructing a supervised learning sample set. The decomposition network module inputs the groundwater level sequence from the supervised learning sample set into the physical information-guided variational mode decomposition network; The component sequence module decomposes the groundwater level sequence into K intrinsic mode component sequences and residual term sequences through a loss function with physical driving consistency constraints. The training module dynamically assembles differentiable simulators from the basic library of differentiable simplified physical simulators for each of the K intrinsic modal component sequences. The goal is to approximate each intrinsic modal component sequence and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. The prediction sequence module inputs the predicted sequence of driving factors for future periods into K well-trained assembled differentiable simulators to perform forward state extrapolation, thereby obtaining the predicted sequence of groundwater level changes for future periods. The forecast report module performs destandardization and uncertainty quantification on the predicted sequence of groundwater level changes for future periods to generate a visual forecast report.

[0016] The beneficial effects of this invention are as follows: By utilizing a variational mode decomposition network guided by physical information and introducing physical-driven consistency constraints into the loss function, the groundwater level sequence is decomposed into multiple intrinsic mode components strongly correlated with specific physical processes. This enables the tracing and decoupling of the physical causes of water level changes, enhancing the model's mechanistic interpretability. For each mode component with a clear physical meaning, a dedicated differentiable simulator is dynamically assembled from the basic library of differentiable simplified physical simulators. Through collaborative training, these simulators accurately fit their respective components while jointly reconstructing the overall water level sequence, thus constructing an interpretable prediction engine that combines data fitting capabilities with physical regularity. This improves the model's extrapolation robustness under unknown scenarios. Through multi-simulator collaborative deduction and uncertainty quantification, reliable visual prediction results are generated. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This is a flowchart of a machine learning-based method for predicting groundwater level changes.

[0019] Figure 2 This is a schematic diagram of a machine learning-based groundwater level change prediction system. Detailed Implementation

[0020] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0021] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0022] Secondly, the term "one embodiment" or "example" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the invention. The appearance of an embodiment in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that mutually excludes other embodiments.

[0023] Reference Figures 1-2 This is one embodiment of the present invention, which provides a groundwater level change prediction method based on machine learning, including the following steps: S1. Standardize the multi-source time series data to obtain a standardized multivariate time series data matrix.

[0024] S1.1. Timestamp alignment is performed on the multi-source time-series data of groundwater level observation, precipitation, evaporation and extraction to generate time-aligned multi-source time-series data, and linear interpolation is used to fill in the gaps to generate filled multi-source time-series data.

[0025] Furthermore, all input data sequences are resampled using a unified time base frequency, such as daily values, to ensure that each data sequence has identical and continuous time point indices on the time axis. For missing values ​​that appear after resampling, a linear interpolation method based on adjacent time point data is used to fill them in, thereby generating complete and aligned time-aligned multi-source time series data in the time dimension.

[0026] S1.2 Perform descriptive statistical analysis on the imputed multi-source time series data to generate statistical characteristics of the mean and standard deviation of each data sequence.

[0027] Furthermore, each independent data sequence in the imputed multi-source time-series data matrix is ​​processed. For each data sequence, the center position of the numerical distribution is evaluated, which is defined as the arithmetic mean of all data points in that column. The dispersion of the data sequence is also evaluated. The dispersion is measured by the square of the difference between each data point and the aforementioned arithmetic mean, the average of the squares is calculated, and the square root of the average is taken. This yields a mean representing the average level of the values ​​for each data sequence and a standard deviation representing the magnitude of the fluctuations in the values. Together, these constitute the statistical characteristics of the mean and standard deviation that describe the core distribution characteristics of the data sequence.

[0028] S1.3. Perform Z-score standardization on the imputed multivariate time series data using the mean and standard deviation of each data sequence to generate a standardized multivariate time series data matrix.

[0029] Furthermore, for each specific value in the imputed multivariate time series data, the mean, corresponding to the center position of its data column, is subtracted. Then, the difference is divided by the standard deviation, which is the quantification value of the dispersion of its data column. Through this operation, all data columns are transformed into a common dimension centered at zero and measured by unit fluctuation, ultimately generating a standardized multivariate time series data matrix.

[0030] S2. Construct a supervised learning sample set.

[0031] S2.1 Determine the historical window length to define the time span of the input features, perform sliding window sampling on the standardized multivariate time series data matrix, and generate data slices.

[0032] Furthermore, constructing a supervised learning sample set requires determining a history window length. The history window length defines the number of past time steps reviewed to construct the input features of the model. On a standardized multivariate time series data matrix, using this history window length as the size, starting from the beginning of the sequence, slide forward one time point at a time in chronological order. Each slide extracts a continuous sub-block of all variable data containing the history window length of time steps. Each such sub-block is a data slice.

[0033] Specifically, continuous time-series data is reconstructed into state segments with clear time boundaries. The setting of the historical window length is not arbitrary, but needs to cover the memory duration and inertial effects of key physical processes in the groundwater level system. For example, for shallow groundwater mainly recharged by precipitation infiltration, the water level response to rainfall may be completed within several days, and the historical window length needs to cover the complete response cycle from the start of rainfall to the peak water level. For deep confined water recharged by long-distance lateral runoff or slow overflow, water level changes may be driven by hydrological and meteorological conditions from months or even years ago, and the required historical window length needs to be significantly extended to capture this long-range correlation. Sliding window sampling is performed on a standardized multivariate time-series data matrix to generate data slices, realizing the transformation from non-stationary time series to stationary state samples.

[0034] It should be noted that the sliding mechanism ensures a smooth transition and partial overlap of information between adjacent states in consecutive time periods. This allows the generated data slices to independently represent the complete state of the system during a specific time period while maintaining the continuity of states in the sequence. This provides machine learning models with sufficient and structured contextual information to learn the dynamic evolution of the system, enabling the models to identify the complete evolution trajectory of different hydrological scenarios, such as the continuous decline in water level during drought, the rapid rise in water level during wet season, and the reversal of water level trends caused by increased mining.

[0035] S2.2 Extract the standardized values ​​of each variable at each time point within the window from each data slice to form a basic feature set. Extract the statistical features of each variable within the window to generate derived features.

[0036] Furthermore, for each data slice, the standardized values ​​of each variable contained within the data slice at each time point within the historical window length are extracted. The standardized values ​​are arranged in order of variable and time to form a basic feature set representing the original observation state of the historical period. For each variable sequence in the data slice, several statistics are obtained within the historical window length. The slope of the linear trend of the value change over time within the window is used to characterize its direction and rate of change, and the range or standard deviation of the value within the window is used to characterize its fluctuation range. The statistics together constitute the derived features describing the behavior pattern of the variable within the window period.

[0037] Specifically, by extracting a basic feature set from each data slice and calculating derived features, feature engineering embodies a multi-layered information mining approach that delves into process mechanisms from surface observations. The basic feature set directly retains the standardized values ​​of variables such as groundwater level observations, precipitation, evaporation, and extraction at each time point within the historical window. This is the most authentic record of the system's instantaneous state. The dynamics of groundwater level are the result of the integration and interaction of various hydrological fluxes over time. Instantaneous values ​​at a single point are insufficient to characterize the cumulative effects, trends, and fluctuation patterns within the window. Statistical features of each variable within the window are extracted separately to generate derived features. It should be noted that the sum of precipitation sequences within a historical window can characterize the total potential recharge for that period; the mean of the extraction sequence can reflect the average extraction intensity; the slope of the linear trend of the water level sequence can quantify the rate of rise or fall of the water level; and the standard deviation of the water level sequence can measure its volatility or the severity of external disturbances. Derived features quantify the behavioral patterns of driving factors and response variables at specific time scales from different dimensions, transforming implicit physical processes into explicit, quantifiable feature indicators. Merging the basic feature set and derived features essentially constructs a comprehensive feature representation that includes both the original high-resolution time-series details and the condensed process statistics. This feature construction method enables the model to learn using both local details and global statistical information simultaneously.

[0038] S2.3. Merge the basic feature set and derived features corresponding to each data slice to form a feature vector that comprehensively represents the state of historical periods, obtain the prediction step size, and define the future time point of the prediction target.

[0039] Furthermore, the values ​​of all elements in the basic feature set and the values ​​of all elements in the derived features are connected in a predetermined order to form a feature vector with higher dimension and richer information. The feature vector comprehensively represents the state of the corresponding historical period. A prediction step size is defined, which indicates the number of time steps between the future time point to be predicted and the end time of the current historical window, thus clarifying the position of the prediction target on the time axis.

[0040] Specifically, merging the basic feature set and derived features corresponding to each data slice to form a comprehensive feature vector is an integration and sublimation of feature information. The comprehensive feature vector integrates the details of the original signal and the generalization of higher-order statistics, forming a feature representation that comprehensively and multi-scale describes the system state over historical periods. It obtains the prediction step size and defines the future time point of the prediction target, serving as a key bridge connecting historical states and future predictions. The determination of the prediction step size needs to be based on hydrogeological conditions and prediction tasks. Its essence is to define a reasonable time delay between cause and effect. For example, for water level prediction in the middle of irrigation collection, a short prediction step size may be needed to capture diurnal scale changes; for assessing the impact of seasonal climate change on groundwater storage, a prediction step size of several months to quarters may be needed. The definition of the prediction step size gives the sample construction a clear time orientation, ensuring that the constructed supervised learning problem is based on the state-response relationship that conforms to the physical time delay law. This allows the model to internalize the specific response time corresponding to different physical processes during learning, thereby improving the accuracy of prediction in the time dimension.

[0041] S2.4 For each historical period corresponding to each data slice, locate the standardized groundwater level observation value after the end of each feature window from the standardized multivariate time series data matrix.

[0042] Furthermore, for each historical time period corresponding to a data slice, based on the defined prediction step size, the specific moment is located in the standardized multivariate time series data matrix, which is a number of time points after the end of the feature window corresponding to that historical time period, separated by the prediction step size. The standardized groundwater level observation value at that moment is then read. Specifically, for each historical period corresponding to each data slice, the standardized groundwater level observation value after the end of each feature window is located from the standardized multivariate time series data matrix. This achieves precise causal alignment between features and labels on the time axis. Based on the prediction step size, the future time point corresponding to the step size after the end of the feature window is found in the standardized multivariate time series data matrix, and the standardized groundwater level observation value at that point is read. This ensures that each historical state feature vector uniquely corresponds to a future water level observation result with causal probability in time series. The strict alignment mechanism abandons the practice of arbitrarily combining historical data with data at a certain future time point, and forces the model to learn a time series mapping relationship that conforms to the actual hydrological response.

[0043] It should be noted that, for example, for a historical window following a heavy rainfall event, the prediction target should correspond to the future water level after the rainfall infiltration and replenishment effect has been fully realized; for a historical window in a period of continuous high-intensity mining, the prediction target should reflect the water level drop that may be caused by subsequent continuous mining, providing the model with physically consistent training samples, which is an important basis for improving the reliability and interpretability of the prediction results.

[0044] S2.5 Pair the feature vectors of each time period with standardized groundwater level observations to construct a supervised learning sample set.

[0045] Furthermore, the feature vector generated for each historical period is combined with the standardized groundwater level observation value that is located and read and corresponds to a future moment of the same historical period. Each such combination constitutes a complete supervised learning sample. By traversing all possible historical period data slices, a supervised learning sample set containing multiple samples is finally constructed.

[0046] Specifically, feature vectors from each time period are paired with standardized groundwater level observations to construct a supervised learning sample set. This completes the final transformation from time-series data to a machine learning-processable format. Each sample consists of a feature vector comprehensively representing the historical state and a future water level value with a clear temporal causal relationship. By traversing all generateable data slices and pairing them, a supervised learning sample set containing a large number of diverse state-response instances is formed. The generated sample set encodes the evolutionary patterns of the groundwater system towards the future under various historical drivers and state combinations. By training on this sample set, it is possible to learn the complex nonlinear mapping of how the system state will evolve from any given historical state under known or predicted future driving scenarios. The application of physically-informed variational mode decomposition networks and differentiable simulators for refined prediction provides high-quality, high-information training data, which is a key prerequisite for achieving high-precision, physically interpretable groundwater level change prediction.

[0047] S3. Input the groundwater level sequence from the supervised learning sample set into the physical information-guided variational mode decomposition network.

[0048] S3.1 Extract the standardized groundwater level observation sequence corresponding to all feature vectors from the supervised learning sample set to form a standardized groundwater level observation sequence.

[0049] Furthermore, the process iterates through each sample in the supervised learning sample set. Each sample consists of a feature vector and a corresponding label for the future standardized groundwater level observation. The feature vector portion of all samples is temporarily set aside, and the labels for the future standardized groundwater level observation corresponding to all samples are extracted and arranged according to the chronological order of the samples in the supervised learning sample set. This process reorganizes these scattered future observations, which were originally the prediction targets, into a continuous and complete sequence of standardized groundwater level observations in the time dimension.

[0050] It's important to note that the structural reorganization of sample information to restore a panoramic temporal perspective of groundwater dynamic evolution is crucial. The supervised learning sample set is constructed for future prediction; its samples are discrete and oriented towards specific future points in time. To deeply understand and decompose the physical mechanisms driving water level changes, it's essential to examine the evolution patterns of the water level sequence itself from a holistic, continuous temporal perspective. By extracting and reorganizing the complete standardized groundwater level observation sequence from the supervised learning sample set in reverse chronological order, the discrete samples prepared for prediction are restored to continuous signals for understanding the mechanisms. This allows subsequent physical information-guided decomposition to act on the complete historical water level evolution trajectory, rather than isolated future prediction points. For example, the model can analyze the frequency component evolution of the water level sequence during a complete flood event, from rise to peak to recession; it can also analyze the low-frequency decay trend of the water level sequence under long-term drought conditions. This processing of continuous sequences provides a foundation for separating modal components corresponding to different physical processes from the time-frequency domain, serving as a key bridge connecting data-driven prediction and physical mechanism explanation. S3.2 By analyzing the correspondence between precipitation and groundwater recharge processes and between extraction and groundwater discharge processes, a set of physical driving sequences corresponding to hydrological processes is defined, and a variational mode decomposition network guided by physical information is constructed.

[0051] Furthermore, by analyzing the correspondence between precipitation and groundwater recharge processes, and between extraction and groundwater discharge processes, a set of physical driving sequences corresponding to hydrological processes is defined. A physical information-guided variational mode decomposition network is constructed. Based on hydrogeological principles, the main physical processes affecting groundwater level dynamics are identified, and external driving variables that can characterize these processes are determined. For example, the precipitation sequence is defined as the dominant physical driving sequence characterizing the groundwater recharge process, and the extraction sequence is defined as the dominant physical driving sequence characterizing the groundwater discharge process. Based on these defined physical driving sequences, a physical information-guided variational mode decomposition network is constructed, which includes a master encoder for receiving and encoding standardized groundwater level observation sequences, a conditional encoder for receiving and encoding the defined physical driving sequences, and a mode decoder for fusing and decomposing the encoded information into physically oriented mode components. Through this structure, external physical driving information is embedded as a guiding signal into the variational mode decomposition process.

[0052] Specifically, by analyzing the correspondence between precipitation and groundwater recharge processes, as well as between extraction and groundwater discharge processes, a set of physical driving sequences corresponding to hydrological processes is defined. A variational mode decomposition network guided by physical information is constructed. The core innovation of this method lies in the fusion of physical information. Traditional signal decomposition methods, such as standard variational mode decomposition, are unsupervised. Their decomposition is performed automatically based on the mathematical characteristics of the signal (such as frequency and bandwidth), and the decomposed components lack clear physical meaning. This method transforms knowledge in the field of hydrogeology into structured guiding signals that can be used by machine learning models. Physical driving sequences corresponding to hydrological processes are defined. For example, the precipitation sequence is explicitly defined as the physical driving sequence characterizing the precipitation infiltration and recharge process, and the extraction sequence is defined as the physical driving sequence characterizing the artificial extraction and discharge process. These are used as supervision signals or reference templates for the decomposition process.

[0053] It should be noted that the physical information-guided variational mode decomposition network receives standardized groundwater level observation sequences to be decomposed, and also receives these predefined physical driving sequences in parallel. The training objective is to learn to decompose the water level sequence into several components, and to force each component to be highly correlated with a specified physical driving sequence in terms of time-frequency characteristics, while being uncorrelated with other driving sequences. The prior knowledge that the decomposition should conform to the physical mechanism is encoded as a constraint in the machine learning loss function. For example, the network is guided to decompose a component that is synchronized with precipitation events in the time domain and reflects the rapid infiltration response characteristics in the frequency domain, and strongly correlates it with the precipitation physical driving sequence. Another component is decomposed to reflect a long-term, stable change trend in the time domain and is low-frequency in the frequency domain, and strongly correlated with the extraction volume physical driving sequence. The physical information-guided variational mode decomposition network achieves a deep integration of data-driven decomposition and physical mechanism interpretation. Each intrinsic mode component decomposed naturally has a physical process label, which fundamentally improves the interpretability of the model.

[0054] S4. Using a loss function with physical-driven consistency constraints, the groundwater level sequence is decomposed into K intrinsic mode component sequences and residual term sequences.

[0055] S4.1. In the loss function based on physical-driven consistency constraints, the signal reconstruction loss term is constrained so that the sum of the K intrinsic mode component sequences and the residual term sequence is consistent with the standardized groundwater level observation sequence.

[0056] Furthermore, during the training process of the variational mode decomposition network guided by physical information, the network first outputs K intrinsic mode component sequences and a residual term sequence. These K intrinsic mode component sequences and the residual term sequence are added point by point in the time dimension to obtain the estimated value of the reconstructed groundwater level sequence. The difference between the estimated value of the reconstructed groundwater level sequence and the original standardized groundwater level observation sequence at each time point is used as the sum of the squares of these differences as the value of the signal reconstruction loss term. The role of the loss term is to ensure the completeness of the decomposition, that is, all decomposed components can completely reconstruct the original signal without distortion.

[0057] Specifically, the cornerstone of ensuring the integrity of the decomposition results is that the sum of all output components of the variational mode decomposition network guided by physical information must accurately reproduce the original input signal. While allowing the network to freely explore the decomposition space to match the physical drive, strict mathematical constraints are imposed to ensure that the decomposition process is a lossless transformation. The variational mode decomposition network guided by physical information can flexibly decompose the standardized groundwater level observation sequence into a variety of possible component combinations, but only those combinations whose sum of components is completely consistent with the original sequence or has a very small error are recognized. This avoids the root cause bias in subsequent predictions due to the loss of information in the decomposition itself. For example, when decomposing a water level sequence that simultaneously includes rapid rainfall response and slow mining trend, the signal reconstruction loss term ensures that the weights and forms of these two processes are reproduced regardless of how the network allocates them.

[0058] S4.2. By using the modal bandwidth loss term in the loss function of the physical-driven consistency constraint, the energy concentration of each intrinsic mode component sequence in the frequency domain is constrained.

[0059] Furthermore, the modal bandwidth loss term in the loss function of the physical-driven consistency constraint is used to constrain the energy concentration of each intrinsic modal component sequence in the frequency domain. For each intrinsic modal component sequence output by the network, the analytic signal of the component is constructed by Hilbert transform to obtain its instantaneous frequency. The dispersion of the spectral energy of the intrinsic modal component sequence around its estimated center frequency is obtained. This is achieved by demodulating the analytic signal of the component by multiplying it by the complex exponent of its center frequency, and then calculating the sum of the squares of the time derivatives of the demodulated signal. The summation is performed on K intrinsic modal component sequences respectively, and the result is used as the value of the modal bandwidth loss term. The role of the loss term is to make each intrinsic modal component sequence as compact and narrow-band as possible in the frequency domain, ensuring the sparsity of the decomposition in the frequency domain and the separation between components.

[0060] Specifically, frequency domain constraints ensure that each intrinsic mode component sequence decomposed has a clearly defined and separable oscillation mode in physical terms. This is derived from the idealized definition of the intrinsic mode function in signal processing theory. A narrow-band amplitude-frequency modulated signal is used, and the instantaneous frequency change rate (i.e., bandwidth) of each intrinsic mode component sequence after demodulation at the center frequency is minimized. The modal bandwidth loss term forces the network to find the components with the most concentrated energy in the frequency domain, which mathematically guarantees the sparsity and orthogonality of the decomposition results, so that different components do not overlap as much as possible in the frequency domain.

[0061] It should be noted that, for groundwater level sequences, this means that fluctuations associated with different physical processes (such as high-frequency rainfall infiltration fluctuations, medium-frequency evaporation impact fluctuations, and low-frequency mining or climate trends) can be automatically separated into different intrinsic mode component sequences. For example, the energy of a rapid rise and fall in water level caused by a rainstorm event will be mainly concentrated in a relatively high-frequency component, while the energy of a slow decline in water level caused by long-term seasonal mining will be concentrated in a low-frequency component.

[0062] S4.3. By quantifying the structural consistency of each intrinsic modal component sequence with the pre-specified physical driving sequence in the time-frequency distribution, the structural differences between the intrinsic modal component sequence and the unspecified physical driving sequence in the time-frequency distribution are quantified.

[0063] Furthermore, by quantifying the structural consistency of each intrinsic modal component sequence with a pre-specified physical driving sequence in the time-frequency distribution, and quantifying the structural differences between the intrinsic modal component sequence and the unspecified physical driving sequence in the time-frequency distribution, for each intrinsic modal component sequence output by the network, the similarity measure in the joint time-frequency domain with a pre-specified dominant physical driving sequence representing a specific hydrophysical process, and the average similarity measure in the joint time-frequency domain between the intrinsic modal component sequence and all other unspecified physical driving sequences, the quantification of structural consistency usually adopts methods such as spectral coherence to measure the degree of matching between the two sequences in the time-frequency energy distribution.

[0064] Specifically, this is achieved by calculating and comparing the similarity metrics of the intrinsic mode component (IMC) sequence with different physical driving sequences in the joint time-frequency domain. For each IMC sequence output by the network, it is necessary to calculate its structural consistency metric in time-frequency distribution with a pre-specified dominant physical driving sequence representing a specific hydrophysical process (such as a precipitation sequence), and the average similarity metric in time-frequency distribution between the IMC sequence and all other unspecified physical driving sequences (such as extraction or evaporation sequences). The time-frequency similarity metrics used here, such as spectral coherence, can assess the stability of phase synchronization between two signals at different frequency components, thereby capturing the degree of matching between the two in time-frequency energy distribution. This metric considers both temporal and frequency domain correlations, and can reveal deeper dynamic relationships between signals than a simple time-domain correlation coefficient. For example, a water level fluctuation component truly related to a precipitation process should not only show a rise in water level after a rainfall event (time domain correlation), but its frequency characteristics (such as response speed and decay period) should also match the spectral characteristics of the rainfall input (frequency domain correlation). S4.4 The alignment loss term separates the physical orientation by maximizing the time-frequency similarity measure between the target intrinsic mode component sequence and the specified physical driving sequence, and minimizing the time-frequency similarity measure with the physical driving sequence.

[0065] Furthermore, the alignment loss term separates physical orientation by maximizing the time-frequency similarity metric between the target intrinsic modal component sequence and the specified physical driving sequence, and minimizing the time-frequency similarity metric with the physical driving sequence. The negative value of the time-frequency similarity metric between each intrinsic modal component sequence and its specified dominant physical driving sequence is taken as an attraction term, and the average time-frequency similarity metric between each intrinsic modal component sequence and the non-specified physical driving sequence is taken as a repulsion term. The sum of these two terms over all K intrinsic modal component sequences is taken as the value of the alignment loss term. This loss term is minimized through the optimization process, thereby forcing each intrinsic modal component sequence to move closer to its specified physical driving sequence and away from other physical driving sequences in terms of time-frequency structure.

[0066] Specifically, the alignment loss term constructs a competitive learning mechanism to drive the physical information-guided variational mode decomposition network to clearly decouple the physical processes. The similarity measure is symbolically transformed and combined to construct the final loss term. For each intrinsic mode component sequence, the time-frequency similarity measure between it and the specified dominant physical driving sequence is negatively valued as an attraction term. Minimizing this term means maximizing the similarity. The average time-frequency similarity measure between the intrinsic mode component sequence and the non-specified physical driving sequence is used as a repulsion term. Minimizing this term means minimizing the association with other driving sequences. The attraction and repulsion terms for all K intrinsic mode component sequences are summed to form the alignment loss term.

[0067] It should be noted that, simulating the separation idea of ​​grouping similar things together, the network sets clear belonging targets and exclusion objects for each intrinsic mode component sequence at the loss function level. In order to minimize the total alignment loss, the network will continuously adjust the parameters so that each output component becomes more and more like its specified driving sequence in terms of time-frequency characteristics, and less and less like other driving sequences. For example, the network will be forced to push the part of the water level sequence that is synchronized with the precipitation event and has a fast response-attenuation spectrum characteristic into the component aligned with the precipitation sequence, and pull the part of the water level sequence that is synchronized with mining activities and has a slow change and low frequency trend characteristic into the component aligned with the mining volume sequence, and ensure that these two components are distinguishable from each other in terms of time-frequency characteristics. By using external physical driving sequences as templates or anchors, the unsupervised decomposition process is guided to produce components with clear physical labels, which fundamentally solves the pain point of unclear physical meaning in traditional signal decomposition methods.

[0068] S4.5 During training, the parameters of the physical information-guided variational mode decomposition network are updated by jointly optimizing the signal reconstruction loss term, modal bandwidth loss term, and alignment loss term through the backpropagation algorithm.

[0069] Furthermore, during training, the signal reconstruction loss term, modal bandwidth loss term, and alignment loss term are jointly optimized using the backpropagation algorithm to update the parameters of the physically-guided variational mode decomposition network. The signal reconstruction loss term, modal bandwidth loss term, and alignment loss term are linearly weighted and summed according to preset weight coefficients to form the total loss function of the physical-driven consistency constraint. The gradient of the total loss function with respect to all trainable parameters in the physically-guided variational mode decomposition network is calculated using the backpropagation algorithm. Then, gradient descent or its variant optimization algorithm is used to update these parameters along the gradient in the opposite direction. By iterating this process, the total loss function is continuously reduced.

[0070] Specifically, the physical information-guided variational mode decomposition network must find the optimal solution under multiple constraints: it must accurately reconstruct the original signal (fidelity), ensure the compact frequency band of each component (mathematical rationality), and satisfy the alignment of the components with the specified physical drive in the time-frequency structure (physical interpretability). This joint optimization avoids the bias that may be caused by optimizing any one objective alone. For example, optimizing only the alignment loss may lead to severe decomposition distortion and failure to reconstruct the original signal; optimizing only the reconstruction loss and bandwidth loss will degenerate into standard VMD, losing physical meaning. By pre-setting weight coefficients (such as assigning appropriate weights to the alignment loss), a balance can be achieved between mathematical fidelity and physical orientation. The gradient descent algorithm updates the network parameters along the negative gradient direction of the total loss function, which essentially guides the network parameter space to evolve in a direction that can simultaneously satisfy the above three constraints. After sufficient training, the parameters of the physical information-guided variational mode decomposition network will stabilize in a certain state, so that its forward propagation can naturally decompose the input water level sequence into components that are both mathematically sound and have clear physical attributes.

[0071] The expression for the signal reconstruction loss term is:

[0072] in, For signal reconstruction loss, For time The reconstructed groundwater level sequence estimate, For time The original standardized groundwater level sequence is as follows. For time.

[0073] The expression for the modal bandwidth loss term is: ; in, For modal bandwidth loss, for, for, The symbol is for partial differentials. For time The partial derivatives, For the Dirac function, The imaginary unit, Pi For time Each intrinsic mode component For the first The center frequencies of each intrinsic mode component.

[0074] The alignment loss term is expressed as follows: ; in, To align the loss terms, for, For the first Each intrinsic mode component has a pre-defined dominant physical driving sequence. The total number of physical drive sequences. For physical drive sequence indexing, For the first A physical drive sequence.

[0075] S4.6. Based on the trained physical information-guided variational mode decomposition network, the standardized groundwater level observation sequence is decomposed into K intrinsic mode component sequences and residual term sequences that are strongly correlated with different physical driving processes.

[0076] Furthermore, based on the trained physical information-guided variational mode decomposition network, the standardized groundwater level observation sequence is decomposed into K intrinsic mode component sequences and residual term sequences, each strongly correlated with different physical driving processes. The standardized groundwater level observation sequence to be decomposed, along with a predefined set of physical driving sequences, is input into the parameter-optimized physical information-guided variational mode decomposition network. The network is then propagated forward, and the output of the network is K intrinsic mode component sequences and one residual term sequence. The features of each intrinsic mode component sequence in the time-frequency domain are strongly correlated with its corresponding specified physical driving sequence, thereby decoupling the original water level sequence into multiple sub-components with clear physical meaning.

[0077] Specifically, once network training is complete, the complex multi-constraint optimization process is condensed into a set of fixed network parameters. The decomposition operation is simplified into a single efficient forward propagation. The standardized groundwater level observation sequence to be decomposed, along with the predefined physical driving sequences (such as precipitation and extraction sequences), are input into the network. The network then automatically outputs K intrinsic mode component sequences and a residual sequence. The decomposition result is now physical: each intrinsic mode component sequence is not only mathematically a component of the original signal, but is also physically clearly labeled as strongly correlated with a specific driving process. For example, one component in the output may have a peak highly synchronized with precipitation events and a spectral characteristic similar to the precipitation sequence; this component is interpreted as a water level fluctuation induced by precipitation recharge. Another component may exhibit a slow linear downward trend, with a change similar to the cumulative extraction curve; this component is interpreted as a trend-driven downward component caused by extraction. The residual sequence may contain some unmodeled noise or very slow background changes.

[0078] S5. For each of the K intrinsic mode component sequences, dynamically assemble a differentiable simulator from the basic library of differentiable simplified physics simulators.

[0079] S5.1 Based on the fundamental physical laws of groundwater flow, construct a basic library of a differentiable simplified physical simulator, analyze the dynamic characteristics of each intrinsic mode component sequence, and extract the feature vector of each intrinsic mode component sequence.

[0080] Furthermore, based on the principles of groundwater dynamics, a series of simplified mathematical models, expressed in the form of differential equations or difference equations, are constructed to describe different typical hydrophysical processes. For example, a differentiable simplified physical simulator is constructed to describe the rapid recharge response using linear reservoir equations; a differentiable simplified physical simulator is constructed to describe solute transport or slow runoff processes using one-dimensional diffusion equations; and a differentiable simplified physical simulator is constructed to describe the recession process using exponential decay functions. Gradient backpropagation is allowed, and statistical characteristics are described for each intrinsic mode component sequence output by the variational mode decomposition network guided by physical information. These characteristics include, but are not limited to, the autocorrelation time of the sequence to describe its memory, the main peak frequency of the power spectrum of the sequence to describe its dominant fluctuation period, and the cross-correlation coefficient between the sequence and the corresponding physical driving sequence to describe its correlation strength. The statistics are combined into a multidimensional vector, which is the feature vector of each intrinsic mode component sequence.

[0081] Specifically, a simplified model library (such as linear reservoirs, diffusion equations, and attenuation models) built based on the fundamental physical laws of groundwater flow (such as Darcy's law and mass conservation) encapsulates the core mathematical essence of hydrological response at different spatiotemporal scales. It analyzes the dynamic characteristics of each intrinsic mode component sequence and extracts feature vectors. In essence, it transforms the clearly directional signal components obtained from the previous physical decomposition into a series of quantifiable and comparable mathematical descriptors. For an intrinsic mode component sequence identified as strongly correlated with rainfall, its feature vector may contain a short autocorrelation time (reflecting rapid response), a high dominant frequency (reflecting frequent fluctuations), and a high cross-correlation coefficient with the rainfall sequence. For a component strongly correlated with mining, its feature vector may exhibit a long autocorrelation time (reflecting persistent impact), a low-frequency dominant power, and a high trend. This transforms the physically ambiguous time series curve into a characteristic fingerprint that can characterize its dynamic mode.

[0082] S5.2. Based on the neural network architecture, set up a simulator selection and parameterization network. The simulator selection and parameterization network outputs the selection weight vector and parameter adjustment vector corresponding to each intrinsic mode component sequence.

[0083] Furthermore, a simulator selection and parameterization network is set up based on a neural network architecture. The simulator selection and parameterization network outputs the selection weight vector and parameter adjustment vector corresponding to each intrinsic mode component sequence. A neural network with a multilayer perceptron structure is constructed as the simulator selection and parameterization network. The input layer dimension matches the dimension of the feature vector of the intrinsic mode component sequence. After nonlinear transformation of several hidden layers, the network branches into two output heads. The first output head outputs a selection weight vector with a length equal to the number of simulators in the base library of differentiable simplified physics simulators. It is normalized by the Softmax function, and the value of each element represents the probability or intensity weight of selecting the corresponding simulator from the base library. The second output head outputs a parameter adjustment vector, the dimension of which matches the number of key physical parameters that need to be adjusted for the selected differentiable simplified physics simulator. When the feature vector of a certain intrinsic mode component sequence is input into the simulator selection and parameterization network, the network outputs the selection weight vector and parameter adjustment vector corresponding to that intrinsic mode component sequence after forward propagation.

[0084] Specifically, leveraging the universal approximation capability of neural networks, a complex mapping from component feature fingerprints to the optimal simulator configuration is learned. The simulator selection and parameterization network receives feature vectors representing the dynamic modes of the components. Its output consists of two parts: a selection weight vector and a parameter adjustment vector. The selection weight vector enables soft selection of the basic simulator library. Instead of rigidly specifying a model, it provides a probability distribution, allowing the advantages of multiple basic models to be utilized in a weighted combination. This enhances the flexibility and expressiveness of the solution, enabling it to handle real-world mixed processes that cannot be perfectly described by a single simplified model. The parameter adjustment vector enables personalized fine-tuning of the physical parameters of the selected model. It recognizes that even if the physical processes are essentially the same, key parameters (such as hydraulic conductivity and specific yield) differ under different locations and aquifer conditions. By learning this adjustment amount through the neural network, the assembled simulator can both inherit the mechanistic framework of the physical model and accurately adapt to the data characteristics of specific components.

[0085] S5.3. Use the selected weight vector to perform a weighted combination of the differentiable simplified physics simulators in the basic library of differentiable simplified physics simulators to generate a preliminary assembled differentiable simulator.

[0086] Furthermore, the selection weight vector of the simulator and the output of the parameterized network is weighted and summed with all differentiable simplified physical simulators in the basic library of differentiable simplified physical simulators. The output of each differentiable simplified physical simulator in the basic library under the current input is multiplied by the weight coefficient in its corresponding selection weight vector. Then, all the weighted outputs are summed to obtain the output of the initially assembled differentiable simulator. This process is essentially a dynamic and soft selection and combination from the basic library to create a composite simulator that is most suitable for simulating the dynamic characteristics of the current intrinsic mode component sequence.

[0087] S5.4. Inject the parameter adjustment vector into the initially assembled differentiable simulator, adjust the internal key parameters of the initially assembled differentiable simulator, and generate the final differentiable simulator.

[0088] Furthermore, key parameters with clear physical meanings are identified from each of the differentiable simplified physical simulators participating in the weighted combination, such as the storage and release coefficient in the linear reservoir model or the diffusion coefficient in the diffusion equation. The parameter adjustment vectors selected by the simulator and output by the parameterized network are added to these identified key parameters according to a preset mapping relationship, thereby performing personalized fine-tuning of the parameters of each component unit in the initially assembled differentiable simulator. The initially assembled differentiable simulator after parameter adjustment is the final differentiable simulator dynamically customized and generated for the current intrinsic mode component sequence.

[0089] S5.5 For each of the K intrinsic mode component sequences, repeatedly perform feature vector extraction, simulator selection and parameterized network processing, weighted combination and parameter adjustment to dynamically assemble K differentiable simulators.

[0090] Furthermore, for each intrinsic mode component sequence output by the variational mode decomposition network guided by physical information, the following operations are performed sequentially: feature vector extraction, feature vector input to simulator selection and parameterization network to obtain the corresponding selection weight vector and parameter adjustment vector, weighted combination of the basic library of differentiable simplified physical simulator using the selection weight vector to generate a preliminary assembled differentiable simulator, and parameter adjustment vector injection into the preliminary assembled differentiable simulator to generate the final differentiable simulator. After traversing all K intrinsic mode component sequences and completing the above process, K final differentiable simulators with different parameters, dynamically assembled for different physical process components, are obtained.

[0091] S6. Co-training is performed with the goal of approximating the sequence of each intrinsic modal component and reconstructing the original water level sequence as a whole, to obtain K fully trained assembled differentiable simulators.

[0092] S6.1. Based on the multi-task learning framework, define the total loss function for collaborative training, and construct the total loss function for collaborative training as a comprehensive index that jointly measures the independent prediction accuracy of each simulator and the overall prediction consistency.

[0093] Furthermore, based on the multi-task learning framework, the total loss function for collaborative training is defined as the method to achieve split optimization. The total loss function is constructed as a comprehensive index that jointly measures the independent prediction accuracy of each simulator and the overall prediction consistency. It sets up dual optimization objectives. The first objective is local fitting loss, which requires each assembled differentiable simulator to independently and accurately predict its corresponding intrinsic mode component sequence. The second objective is global reconstruction loss, which requires that the sum of the outputs of all simulators must accurately restore the original complete groundwater level sequence.

[0094] Specifically, on the one hand, it incentivizes each simulator to focus on learning and simulating the specific physical process it is responsible for (such as rainfall response or mining trends) and pursue local optima; on the other hand, it forces all simulator outputs to be harmoniously superimposed through global reconstruction loss, avoiding mutually exclusive or conflicting outputs from each simulator when pursuing local optima. For example, if one simulator over-predicts the supply amount while another simulator fails to adjust its discharge amount accordingly, it will lead to overall reconstruction distortion. The total loss function is a weighted sum of these two objectives, like a conductor who must ensure that each musician (individual simulator) plays its part well (fits its own components) and that the entire orchestra (all simulators) plays a harmonious overall melody (reconstructs the original sequence).

[0095] S6.2 Each differentiable simulator and its corresponding dominant physics driving sequence generate an eigenmode component sequence by performing state transitions and response outputs through the simplified physical laws inherent within it.

[0096] Furthermore, each differentiable simulator and its corresponding dominant physics driving sequence, through the simplified physical laws inherent within it, executes state transitions and response outputs to generate an eigenmode component sequence. This is the concrete execution and manifestation of physical knowledge in the prediction stage. Each assembled differentiable simulator is essentially a parameterized differential or difference equation solver. It receives the corresponding dominant physics driving sequence as an external input or forcing term, and performs step-by-step state calculations based on the simplified physical laws encoded within it (such as linear storage-drainage relations and diffusion equations), outputting the predicted eigenmode component sequence. The simulator's operation process strictly follows the mathematical expression of physical laws, and the output results can be traced back to specific input drivers and physical equations.

[0097] Specifically, for example, a simulator responsible for rainfall response may internally be a linear reservoir model. It receives the rainfall sequence as input and obtains the response component of the groundwater level through a differential equation representing the infiltration and discharge process. This makes the prediction process no longer a black-box data mapping, but a white-box inference with clear physical causal relationships, which greatly enhances the reliability and extrapolation ability of the model, especially when facing new driving scenarios that have not appeared in historical data.

[0098] S6.3. The global sequence reconstruction loss measures the fidelity of the standardized groundwater level observation sequence in terms of overall waveform and numerical value after the intrinsic modal component sequences are superimposed and combined.

[0099] Furthermore, the global sequence reconstruction loss measures the fidelity of the standardized groundwater level observation sequence after the intrinsic modal component sequences are superimposed and combined in terms of overall waveform and numerical values. This ensures the final effective assembly line and consistency constraints of the model prediction. The intrinsic modal component prediction sequences output by each of the K differentiable simulators are added point by point to obtain the overall reconstruction prediction of the original complete groundwater level sequence. Then, the difference measure (such as mean square error) between the overall reconstruction prediction sequence and the actual standardized groundwater level observation sequence is calculated.

[0100] Specifically, starting from the ultimate application goal (accurately predicting the total water level), a final constraint is imposed on the front-end physical decomposition process and the independent predictions of each physical simulator. This constraint mandates that each simulator, while optimizing its own local objective, must consider the compatibility of its predictions with those of other simulators and its contribution to the final sum. Without this constraint, each simulator might merely pursue overfitting to its own target component, ignoring the overall rationality of its prediction. This could lead to offsetting or amplifying errors when the predicted values ​​of each component are superimposed, making it impossible to accurately reconstruct the total water level. The global reconstruction loss term ensures that the systematic principle that the whole is greater than the sum of its parts is upheld. S6.4 During the training iteration, based on the total loss function, the gradient is calculated using the backpropagation algorithm, and the parameters of the simulator selection and parameterization network, as well as the parameters inside the differentiable simulator, are updated. Furthermore, the training iteration process, which uses the total loss function and backpropagation to calculate gradients to update parameters, is the key engine for achieving end-to-end optimization and fusing data and physics. Its core lies in utilizing differentiable programming, allowing gradients to permeate the entire computational graph: starting from the total loss function, the gradient can propagate back through global summation operations, the physics calculations of each simulator, the simulator's assembly parameters (weight selection and adjustment), all the way back to the simulator's selection and parameterization of the network's weights. The gradient update expression embodies the optimization process of minimizing the total loss function by iteratively fine-tuning all trainable parameters, achieving automation and joint optimization of data-driven calibration of physical model parameters and model structure (weight selection).

[0101] Specifically, for example, gradient signals guide the simulator's selection and parameterization network to adjust their selection weights, giving higher weights to models in the base simulator library that can more effectively reduce the total loss. Simultaneously, gradients also adjust key physical parameters (such as permeability coefficients) within each differentiable simplified physics simulator, making them more accurately reflect the actual properties of the current aquifer. This end-to-end gradient flow seamlessly connects physical constraints (the simulator's equations) with the data fitting objective (total loss), enabling the physical model to learn optimal parameters and combinations from the data, and allowing the data to be interpreted more reasonably through the constraints of the physical model.

[0102] The gradient update expression is: ; in, In the first The set of all trainable parameters at the time of the training iteration In the first The set of all trainable parameters at the time of the training iteration For learning rate, loss function In parameters The gradient at that point.

[0103] S6.5 Set the training termination condition. When the total loss function value reaches the preset number of iterations, the training process is completed, and K assembled differentiable simulators with fixed parameters are obtained.

[0104] Furthermore, setting training termination conditions and obtaining K assembled differentiable simulators with fixed parameters after the conditions are met marks the transition of the model from the training phase to the stable application phase. The training termination condition is usually set to the total loss function no longer decreasing on the validation set or reaching the preset maximum number of iterations. After training is completed, the simulator selection, the weights of the parameterized network, and the internal parameters of each differentiable simulator all reach a stable state. The obtained K assembled differentiable simulators are the finalized prediction components after co-training optimization.

[0105] Specifically, each simulator undergoes rigorous training with dual objectives (local fitting and global consistency), enabling it to independently and accurately simulate its corresponding physical process components, and to work in collaboration with other simulators to accurately reconstruct the total water level. The simulator parameters have been calibrated to the optimal level by data, and its structure (reflected by the selection of weights) is adapted to the characteristics of the groundwater system to be predicted.

[0106] S7. Input the predicted sequence of driving factors for future periods into K well-trained assembled differentiable simulators to perform forward state extrapolation and obtain the predicted sequence of groundwater level changes for future periods.

[0107] S7.1 Based on the acquired driving factor forecast sequence, the precipitation forecast sequence in the driving factor forecast sequence is guided to a fully trained assembled differentiable simulator that simulates the effect of rainfall on groundwater recharge.

[0108] Furthermore, based on the precise data routing established during the training phase, the physical isomorphism matching between prediction-driven and simulator functions is achieved. In the early training phase, the variational mode decomposition network guided by physical information has decomposed the water level sequence into multiple components. Each component is forced to align with a specific physical driving sequence (such as precipitation sequence) in terms of time-frequency features. The co-training process specifically optimizes a simulator to accurately fit this component that is strongly correlated with precipitation.

[0109] Specifically, during the forecast phase, precipitation forecast sequences are directed to an assembled differentiable simulator that has been specially trained to simulate rainfall replenishment processes. This ensures that the driving data is processed by the model that best understands its physical effects, avoiding the confusion and coupling errors that may be caused by inputting mixed driving data into a single model. For example, forecast data for a rainstorm will only enter the simulator responsible for simulating rapid infiltration response, and will not be mistakenly used to affect the simulator simulating long-term mining trends. This ensures clear separation and independent deduction of physical processes at the source of the forecast, which is the foundation for achieving interpretable forecasts.

[0110] S7.2 Based on the current memory of the dynamic state, combined with the input precipitation forecast sequence, generate a water level component sequence that reflects the contribution of future rainfall replenishment.

[0111] Furthermore, based on the current dynamic state of memory and combined with the input precipitation forecast sequence, a future water level component sequence is generated, and a fully trained assembled differentiable simulator is used as a dynamic system model. Its internal state (such as water storage volume and water level value of the previous period) encodes the historical memory of the system up to the current moment. Before the prediction begins, the simulator needs to be driven to the current moment using the most recent historical observation data so that its internal state is synchronized with the current state of the real system.

[0112] Specifically, starting from this dynamic state that aligns with current reality, the simulator performs step-by-step calculations based on the input precipitation forecast sequence and its embedded physical laws (such as the storage and release equations of a linear reservoir) to generate future replenishment contribution components. It respects the path dependence of the physical system; the response of the groundwater level to rainfall depends not only on future rainfall but also strongly on current initial conditions such as soil saturation and groundwater depth. For example, after a prolonged drought, the replenishment generated by the same amount of rainfall is far less than after a wet period. Through state memory-based deduction, the prediction results naturally incorporate the influence of current system conditions, making the predictions closer to physical reality and improving the accuracy of predictions under changing initial conditions.

[0113] S7.3 The extraction volume forecast sequence in the driving factor forecast sequence is guided to a fully trained assembled differentiable simulator that simulates the effect of extraction on groundwater discharge. Based on the current memorized dynamic state, combined with the input extraction volume forecast sequence, a water level component sequence reflecting the impact of future extraction and discharge is generated.

[0114] Furthermore, a routing and extrapolation logic similar to, but independent and parallel to, the precipitation sequence was adopted. The extraction volume forecast sequence was guided to a fully trained, modular, differentiable simulator specifically designed to simulate the effects of extraction and discharge. This simulator, also based on its current memorized dynamic state (reflecting the current extraction response characteristics of the aquifer, such as the unit rate of decline), physically extrapolates future input extraction volume sequences, generating independent water level discharge influence components, thus achieving resolvable prediction scenarios.

[0115] Specifically, because the extraction component is simulated separately, decision-makers can easily assess the independent impact of different extraction management schemes (such as limited extraction and increased extraction) on future water levels without having to rerun the entire complex model. For example, it is possible to quickly compare the differences in water level trend components under two scenarios: maintaining current extraction levels versus reducing extraction after implementing water-saving plans. This provides water resource management with a powerful and transparent if-then scenario analysis tool, making the predictive model not only a forecasting tool but also a decision support system.

[0116] S7.4 After receiving a continuous sequence of driving factor forecasts, the simulator iteratively updates its internal state and generates the corresponding water level component sequence step by step.

[0117] Furthermore, the core inference mechanism of the simulator receives continuous driving sequences and iteratively updates them. Each fully trained assembled differentiable simulator runs as a discrete-time dynamic model. Within each time step, the simulator performs the following operations: reads the driving input values ​​of the current time step; calculates the state update amount and output water level component values ​​of the current time step based on its internal physical equations and the internal state at the end of the previous time step; updates the internal state to the new values ​​after the current time step ends, and iterates and repeats this process over time to generate a water level component sequence covering the entire future prediction period.

[0118] Specifically, it transforms complex long-term prediction problems into a series of manageable, physics-based short-step calculations. It strictly follows the causal temporal logic of physical systems, meaning that the current state is the result of the interaction between the previous state and the current input. This mechanism ensures the stability of the simulation and can naturally capture the dynamic cumulative effects of the process (such as the continuous decline in water storage caused by consecutive droughts) or hysteresis effects (such as the delayed response of precipitation infiltration). Through iterative time-step calculations, the predictive capabilities of the physical simulator can be robustly extended over time.

[0119] S7.5. All water level component sequences generated by each simulator throughout the entire future prediction period are superimposed point by point in the time dimension with the residual water level sequence estimated by trend extrapolation, and integrated into a complete standardized groundwater level change prediction sequence for the future period.

[0120] Furthermore, synthesizing all independently extrapolated water level component sequences with the trend-extrapolated residual water level sequence is the final and crucial step in restoring a complete water level prediction. The water level component sequences generated by each simulator represent the contributions of interpretable physical processes (such as rainfall replenishment and extraction discharge). By using simple time series models (such as linear extrapolation and ARIMA models) to predict the future of the residual sequence obtained from historical decomposition, a residual water level sequence is obtained. This sequence contains unmodeled slow trends, periodic noise, or other systematic biases. The standardized groundwater level change prediction sequence is obtained by algebraically adding the predicted values ​​of all K water level component sequences to the predicted value of the residual water level sequence at each future time point. This adheres to the linear superposition assumption that the whole equals the sum of its parts, which is the theoretical cornerstone of the entire decomposition-simulation-synthesis framework.

[0121] Specifically, it provides a complete and unbiased overall forecast because all identified physical process contributions and unexplained residual variations are taken into account, maintaining an interpretable framework for the forecast results. Decision-makers not only receive the predicted total water level but can also clearly trace how much of that total value is attributable to anticipated rainfall, how much to planned mining, and how much to other trends or uncertainties. This transparent forecast composition is invaluable for understanding the forecast results, assessing the impact weights of different factors, and conducting risk management, transforming the model output from a single numerical value into an information-rich basis for decision-making.

[0122] S8. Perform destandardization and uncertainty quantification on the predicted sequence of groundwater level changes in future periods to generate a visual prediction report.

[0123] S8.1 Perform inverse Z-score transformation on the standardized groundwater level change prediction sequence to restore the standardized groundwater level change prediction sequence to a groundwater level change prediction sequence with actual physical dimensions.

[0124] Furthermore, the statistical characteristics of the mean and standard deviation calculated during the standardization process of the original groundwater level observation sequence are utilized. Specifically, each value in the standardized groundwater level change prediction sequence is multiplied by its corresponding standard deviation statistical characteristic, and then the corresponding mean statistical characteristic is added. Through this inverse linear transformation, the previously normalized prediction values, centered at zero and on a unit variance scale, are restored to a groundwater level change prediction sequence with actual physical units (such as meters).

[0125] Specifically, the model's internal mathematical calculations are converted back into numerical values ​​with real-world physical meaning that end users (such as water managers and geologists) can intuitively understand and directly apply. This completes the final step in the conversion from model space to the physical world, enabling predicted values ​​to be directly compared and evaluated with historical observation records and management thresholds (such as warning water levels), which is a prerequisite for the practical application of prediction results.

[0126] S8.2 Apply random perturbations that conform to prior knowledge to the key parameters of the fully trained assembled differentiable simulator to generate multiple sets of parameter settings.

[0127] Furthermore, applying random perturbations consistent with prior knowledge to the key parameters of the fully trained assembled differentiable simulator demonstrates the integration of Bayesian thinking into deterministic physical models. Key parameters refer to those parameters that have clear physical meaning in the simulator and significantly affect the simulation results, such as the decay time constant in a linear reservoir model and the diffusion coefficient in a diffusion model. The random perturbations consistent with prior knowledge are not arbitrary but rather reasonable estimates of parameter uncertainties based on domain knowledge. For example, the prior uncertainty of the permeability coefficient might be considered to follow a log-normal distribution; the uncertainty range of the specific yield might be set based on empirical lithological values.

[0128] Specifically, by independently sampling from these pre-defined probability distributions, a set of random perturbation increments (such as additive or multiplicative perturbations) relative to their post-training determined values ​​is generated for each key parameter, thus forming a new parameter setting. This independent sampling process is repeated multiple times to generate multiple sets of parameter settings. This acknowledges and quantifies the inherent uncertainty of the physical model parameters. Due to the spatial heterogeneity of aquifer properties and the finiteness of observational data, the optimal parameters obtained through training are essentially an estimate of the true parameters and have confidence intervals. By applying reasonable perturbations, multiple equally likely model instances are generated.

[0129] S8.3. Using a fully trained assembly-type differentiable simulator with multiple sets of parameter settings, perform multiple independent forward state extrapolations on the future driving factor prediction sequence to generate multiple groundwater level change prediction sequences, forming a groundwater level change prediction sequence set.

[0130] Furthermore, the specific method for using simulators with multiple sets of parameter settings to perform multiple independent simulations is integrated Monte Carlo simulation. For each set of parameter settings, these parameters are injected into the corresponding positions of K fully trained assembled differentiable simulators, thereby temporarily constructing K simulator instances with fine-tuned parameters. Using these simulator instances with fine-tuned parameters, the forward state inference process is completely re-executed: the predicted sequence of driving factors for future periods is input, and each simulator iteratively calculates based on its current memorized dynamic state to generate a set (K components) of water level component sequences. After synthesis, a complete groundwater level change prediction sequence under this set of parameter settings is obtained. This inference process is repeated independently and completely for each different set of parameter settings. After traversing all parameter sets, a set of multiple prediction sequences is obtained, namely the groundwater level change prediction sequence set, where each sequence represents a reasonable inference of the future water level under a possible parameter reality.

[0131] S8.4 At each future time point in the groundwater change prediction sequence set, extract the specific quantiles of the distribution of all predicted values ​​to obtain the prediction interval for the future time point. Use the groundwater level change prediction sequence with actual physical dimensions as the median and the prediction interval as the confidence band to draw a time series curve.

[0132] Furthermore, to provide an intuitive visualization of the uncertainty in predictions, each generated groundwater level change prediction sequence corresponds to a series of future time points. For any given future time point, the predicted values ​​of all prediction sequences at that time point are extracted to form a sample set. This sample set is then sorted, and values ​​corresponding to specific probability quantiles (e.g., the 5th and 95th percentiles) are identified. These two quantile values ​​constitute a prediction interval for that future time point, such as a 90% confidence interval. Based on the current model and parameter knowledge, there is a 90% probability that the actual future water level value will fall within this interval. The same calculation is performed for all future time points to obtain a series of prediction intervals that change over time.

[0133] Specifically, when plotting time series graphs, the predicted sequence of groundwater level changes with actual physical dimensions generated from the original optimal parameters is plotted as the median or best estimate line. The upper and lower bounds of the prediction interval for each time point are filled with different colors or line types (such as light-colored bands) to form a confidence band around the median. This compresses and refines the abstract, multidimensional prediction uncertainty (reflected in the divergence of multiple prediction trajectories) into an easily understandable graphical language—the confidence interval. Decision-makers can clearly see the overall trend of the prediction (median) and the range of uncertainty around this trend (confidence band). For example, they can clearly identify whether the uncertainty is increasing sharply at a certain stage of the prediction period, which may correspond to an increase in the uncertainty of the driving factor prediction itself, or a weakening of the model's extrapolation ability during that period.

[0134] S8.5 Add legends, axis labels, and forecast step size descriptions to the time series curve to generate a visual forecast report.

[0135] Furthermore, adding legends, axis labels, and forecast step size descriptions to time series graphs to generate the final visualized forecast report is the final step in information integration and effective communication. Legends clearly distinguish the median forecast line, confidence bands, and different scenario lines (such as different mining schemes). Axis labels must clearly indicate physical quantities and units (such as time (year-month), groundwater level (meters)). The forecast step size description clearly indicates the start time of the forecast, the forecast duration, and the future driving scenario on which the forecast is based. Integrating these elements into one or a set of charts constitutes the visualized forecast report.

[0136] This example also provides a machine learning-based groundwater level change prediction system, including: a construction module that standardizes multi-source time series data to obtain a standardized multivariate time series data matrix and constructs a supervised learning sample set; The decomposition network module inputs the groundwater level sequence from the supervised learning sample set into the physical information-guided variational mode decomposition network; The component sequence module decomposes the groundwater level sequence into K intrinsic mode component sequences and residual term sequences through a loss function with physical driving consistency constraints. The training module dynamically assembles differentiable simulators from the basic library of differentiable simplified physical simulators for each of the K intrinsic modal component sequences. The goal is to approximate each intrinsic modal component sequence and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. The prediction sequence module inputs the predicted sequence of driving factors for future periods into K well-trained assembled differentiable simulators to perform forward state extrapolation, thereby obtaining the predicted sequence of groundwater level changes for future periods. The forecast report module performs destandardization and uncertainty quantification on the predicted sequence of groundwater level changes for future periods to generate a visual forecast report.

[0137] In summary, this invention utilizes a variational mode decomposition network guided by physical information, introducing a physics-driven consistency constraint into the loss function to decompose the groundwater level sequence into multiple intrinsic mode components strongly correlated with specific physical processes. This enables the tracing and decoupling of the physical causes of water level changes, enhancing the model's mechanistic interpretability. For each mode component with a clear physical meaning, a dedicated differentiable simulator is dynamically assembled from a library of differentiable simplified physical simulators. Through collaborative training, these simulators accurately fit their respective components while jointly reconstructing the overall water level sequence, thus constructing an interpretable prediction engine that combines data fitting capabilities with physical regularity. This improves the model's extrapolation robustness under unknown scenarios. Through multi-simulator collaborative deduction and uncertainty quantification, reliable visual prediction results are generated. It should be noted that the above examples are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred examples, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A groundwater level change prediction method based on machine learning, characterized in that: This includes standardizing multi-source time series data to obtain a standardized multivariate time series data matrix and constructing a supervised learning sample set; The groundwater level sequence from the supervised learning sample set is input into a physical information-guided variational mode decomposition network; The groundwater level sequence is decomposed into K intrinsic mode component sequences and residual term sequences by using a loss function with physical driving consistency constraints. For each of the K intrinsic modal component sequences, a differentiable simulator is dynamically assembled from the basic library of differentiable simplified physical simulators. The goal is to approximate each intrinsic modal component sequence and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. The predicted sequence of driving factors for future periods is input into K well-trained assembled differentiable simulators for forward state extrapolation to obtain the predicted sequence of groundwater level changes for future periods. The prediction sequence of groundwater level changes in future periods is de-standardized and uncertainty is quantified to generate a visual prediction report.

2. The groundwater level change prediction method based on machine learning as described in claim 1, characterized in that: Standardizing multi-source time series data to obtain a standardized multivariate time series data matrix includes the following steps: The time-series data of groundwater level observation, precipitation, evaporation and extraction are time-stamped to generate time-aligned multi-source time-series data, and then filled in using a linear interpolation method to generate filled multi-source time-series data. Descriptive statistical analysis was performed on the imputed multi-source time series data to generate statistical characteristics of the mean and standard deviation of each data series. The mean and standard deviation of each data sequence are used to perform Z-score standardization on the imputed multivariate time series data to generate a standardized multivariate time series data matrix.

3. The groundwater level change prediction method based on machine learning as described in claim 2, characterized in that: Constructing a supervised learning sample set includes the following steps: The historical window length is determined to define the time span of the input features, and sliding window sampling is performed on the standardized multivariate time series data matrix to generate data slices; The standardized values ​​of each variable at each time point within the window are extracted from each data slice to form a basic feature set. Statistical features of each variable within the window are then extracted to generate derived features. The basic feature set and derived features corresponding to each data slice are merged to form a feature vector that comprehensively represents the state of historical periods, thereby obtaining the prediction step size and defining the future time point of the prediction target. For each data slice corresponding to a historical period, the standardized groundwater level observation value after the end of each feature window is located from the standardized multivariate time series data matrix; The feature vectors of each time period are paired with standardized groundwater level observations to construct a supervised learning sample set.

4. The groundwater level change prediction method based on machine learning as described in claim 3, characterized in that: The groundwater level sequence from the supervised learning sample set is input into a physically-guided variational mode decomposition network, including the following steps: The standardized groundwater level observation sequence is formed by extracting all feature vectors corresponding to the supervised learning sample set. By analyzing the correspondence between precipitation and groundwater recharge processes, as well as between extraction and groundwater discharge processes, a set of physical driving sequences corresponding to hydrological processes is defined, and a variational mode decomposition network guided by physical information is constructed.

5. The groundwater level change prediction method based on machine learning as described in claim 4, characterized in that: The groundwater level sequence is decomposed into K intrinsic mode component sequences and residual term sequences using a loss function driven by physical consistency constraints, including the following steps: The signal reconstruction loss term in the loss function based on physical-driven consistency constraints is constrained to ensure that the sum of the K intrinsic mode component sequences and the residual term sequence is consistent with the standardized groundwater level observation sequence. The energy concentration of each intrinsic mode component sequence in the frequency domain is constrained by the modal bandwidth loss term in the loss function of the physical-driven consistency constraint. The structural differences between intrinsic modal component sequences and unspecified physical driving sequences in the time-frequency distribution are quantified by quantifying the structural consistency between each intrinsic modal component sequence and the pre-specified physical driving sequence in the time-frequency distribution. The alignment loss term separates the physical orientation by maximizing the time-frequency similarity measure between the target intrinsic mode component sequence and the specified physical driving sequence, and minimizing the time-frequency similarity measure with the physical driving sequence. During training, the parameters of the physical information-guided variational mode decomposition network are updated by jointly optimizing the signal reconstruction loss term, modal bandwidth loss term, and alignment loss term through the backpropagation algorithm. Based on a variational mode decomposition network guided by trained physical information, the standardized groundwater level observation sequence is decomposed into K intrinsic mode component sequences and residual term sequences that are strongly correlated with different physical driving processes.

6. The groundwater level change prediction method based on machine learning as described in claim 5, characterized in that: For each of the K eigenmode component sequences, a differentiable simulator is dynamically assembled from the base library of differentiable simplified physics simulators, including the following steps: Based on the fundamental physical laws of groundwater flow, a basic library of a differentiable simplified physical simulator is constructed. The dynamic characteristics of each intrinsic mode component sequence are analyzed, and the feature vector of each intrinsic mode component sequence is extracted. Based on a neural network architecture, a simulator selection and parameterization network is set up. The simulator selection and parameterization network outputs the selection weight vector and parameter adjustment vector corresponding to each intrinsic mode component sequence. By using a weighted combination of differentiable simplified physics simulators in the basic library of differentiable simplified physics simulators, a preliminary assembled differentiable simulator is generated. The parameter adjustment vector is injected into the initially assembled differentiable simulator to adjust the key internal parameters of the initially assembled differentiable simulator and generate the final differentiable simulator. For each of the K intrinsic mode component sequences, feature vector extraction, simulator selection and parameterized network processing, weighted combination and parameter adjustment are repeatedly performed to dynamically assemble K differentiable simulators.

7. The groundwater level change prediction method based on machine learning as described in claim 6, characterized in that: The goal is to approximate the sequences of each intrinsic mode component and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. This includes the following steps: Based on the multi-task learning framework, the total loss function of collaborative training is defined and constructed as a comprehensive index that jointly measures the independent prediction accuracy of each simulator and the overall prediction consistency. Each differentiable simulator and its corresponding dominant physics driving sequence generate an eigenmode component sequence by performing state transitions and response outputs through the simplified physical laws inherent within it. The global sequence reconstruction loss measures the fidelity of the normalized groundwater level observation sequence in terms of both overall waveform and numerical value after the intrinsic modal component sequences are superimposed and combined. During the training iteration, the gradient is calculated based on the total loss function using the backpropagation algorithm, and the parameters of the simulator selection and parameterization network, as well as the parameters inside the differentiable simulator, are updated. Set the training termination condition. When the total loss function value reaches the preset number of iterations, the training process is completed, and K assembled differentiable simulators with fixed parameters are obtained.

8. The groundwater level change prediction method based on machine learning as described in claim 7, characterized in that: The predicted sequence of driving factors for future periods is input into K well-trained assembled differentiable simulators for forward state extrapolation to obtain the predicted sequence of groundwater level changes for future periods, including the following steps: Based on the acquired driving factor forecast sequence, the precipitation forecast sequence in the driving factor forecast sequence is guided to a fully trained assembled differentiable simulator that simulates the effect of rainfall on groundwater recharge. Based on the current dynamic state of memory, combined with the input precipitation forecast sequence, a water level component sequence reflecting the contribution of future rainfall replenishment is generated; The extraction volume forecast sequence in the driving factor forecast sequence is guided to a fully trained modular differentiable simulator that simulates the effect of extraction on groundwater discharge. Based on the current memorized dynamic state and combined with the input extraction volume forecast sequence, a water level component sequence reflecting the impact of future extraction and discharge is generated. After receiving a continuous sequence of driving factor forecasts, the simulator iteratively updates its internal state and generates the corresponding water level component sequence step by step. All water level component sequences generated by each simulator throughout the entire future prediction period are superimposed point by point along the time dimension with the residual water level sequence estimated by trend extrapolation, and integrated into a complete standardized groundwater level change prediction sequence for the future period.

9. The groundwater level change prediction method based on machine learning as described in claim 8, characterized in that: The prediction series of groundwater level changes for future periods is de-standardized and uncertainty quantified to generate a visual prediction report, including the following steps: Perform an inverse Z-score transformation on the standardized groundwater level change prediction sequence to restore the standardized groundwater level change prediction sequence to a groundwater level change prediction sequence with actual physical dimensions. Apply random perturbations that conform to prior knowledge to the key parameters of the fully trained assembled differentiable simulator to generate multiple sets of parameter settings; Using a fully trained assembly-type differentiable simulator with multiple sets of parameter settings, multiple independent forward state extrapolations are performed on the future driving factor prediction sequence to generate multiple groundwater level change prediction sequences, forming a groundwater level change prediction sequence set. At each future time point in the groundwater change prediction sequence set, specific quantiles of the distribution of all predicted values ​​are extracted to obtain the prediction interval for the future time point. The groundwater level change prediction sequence with actual physical dimensions is used as the median, and the prediction interval is used as the confidence band to draw a time series curve. Add legends, axis labels, and forecast step size descriptions to the time series curve to generate a visual forecast report.

10. A machine learning-based groundwater level change prediction system, based on the machine learning-based groundwater level change prediction method according to any one of claims 1 to 9, characterized in that: This includes a construction module that standardizes multi-source time-series data to obtain a standardized multivariate time-series data matrix and constructs a supervised learning sample set. The decomposition network module inputs the groundwater level sequence from the supervised learning sample set into the physical information-guided variational mode decomposition network; The component sequence module decomposes the groundwater level sequence into K intrinsic mode component sequences and residual term sequences through a loss function with physical driving consistency constraints. The training module dynamically assembles differentiable simulators from the basic library of differentiable simplified physical simulators for each of the K intrinsic modal component sequences. The goal is to approximate each intrinsic modal component sequence and reconstruct the original water level sequence as a whole through collaborative training, resulting in K fully trained assembled differentiable simulators. The prediction sequence module inputs the predicted sequence of driving factors for future periods into K well-trained assembled differentiable simulators to perform forward state extrapolation, thereby obtaining the predicted sequence of groundwater level changes for future periods. The forecast report module performs destandardization and uncertainty quantification on the predicted sequence of groundwater level changes for future periods to generate a visual forecast report.

Images