A groundwater level multi-scale prediction method for meteorological observation sparse areas

Abstract

This invention belongs to the field of hydrology and water resources technology, specifically disclosing a multi-scale prediction method for groundwater levels in areas with sparse meteorological observations. The invention preprocesses historical groundwater level time-series data and constructs a lag candidate feature library, employing the HSIC-Lasso algorithm to adaptively select the optimal input features. It then uses LPSO-Adaptive-VMD to achieve adaptive multi-scale decomposition of the groundwater level sequence, dividing it into low-frequency trend components and mid-to-high-frequency fluctuation components. The N-BEATS and BiGRU models are used to predict the two types of components in parallel. Finally, the LightGBM model is used to learn the nonlinear coupling relationship between the multi-scale components, completing adaptive fusion reconstruction and outputting the prediction results. This invention uses only historical groundwater level sequences as the sole input, requiring no meteorological driving data, and has the advantages of strong generalization, high accuracy, and good robustness. It can be directly extended to dynamic groundwater level prediction in remote grasslands, deserts, and other data-scarce areas.

Description

Technical Field

[0001] This invention belongs to the field of hydrology, water resources, and eco-hydrology technology, and specifically discloses a multi-scale prediction method for groundwater levels in areas with sparse meteorological observations. It involves a combination of machine learning and hydrological time-series analysis techniques, and is particularly suitable for high-precision dynamic prediction of groundwater levels in arid and semi-arid grasslands and deserts where meteorological observation stations are sparse and hydrogeological data is scarce. Background Technology

[0002] Groundwater is a core water resource for maintaining ecosystem stability and ensuring agricultural and livestock production and residents' lives in arid and semi-arid regions. Accurate dynamic prediction of groundwater levels is the core foundation for regional water resource regulation, grassland ecological protection, and drought risk prevention and control.

[0003] Existing groundwater level prediction technologies are mainly divided into two categories: traditional physical process numerical models and data-driven models. Traditional physical process numerical models, based on hydrogeological mechanisms, can provide mechanistic explanations. However, when applied at a regional scale, they face challenges such as difficulty in obtaining aquifer parameters, uncertain boundary conditions, and high model calibration costs. This is especially true in grassland hinterlands where monitoring wells are sparse and hydrogeological data is incomplete, making model construction and long-term operation extremely difficult and hindering large-scale deployment. Data-driven models, including machine learning and deep learning methods, are widely used for hydrological time series prediction due to their strong nonlinear fitting capabilities. Current mainstream solutions often use meteorological factors such as precipitation, temperature, and evapotranspiration along with groundwater level sequences as model inputs. While this can improve accuracy in areas with complete data, in grassland areas with sparse meteorological observation stations, precipitation exhibits strong spatiotemporal heterogeneity, and observations of high-intensity recharge events are lacking. High-resolution meteorological data is difficult to obtain stably, leading to a sharp decline in model generalization ability and making practical deployment impossible.

[0004] Self-driven prediction models relying solely on historical groundwater level sequences can circumvent the core pain point of missing meteorological data, but their existing technology still suffers from three major drawbacks: the construction of input lag features depends on human experience, and fixed windows or full-lag inputs easily introduce a large number of redundant features, reducing model training efficiency and cross-well point generalization ability; groundwater levels exhibit significant multi-scale and non-stationarity, and a single model often struggles to simultaneously account for long-term trends, seasonal cycles, and short-term disturbances; prediction results after multi-scale decomposition generally use linear superposition reconstruction, which cannot characterize the cross-scale nonlinear coupling relationship between different frequency band components, resulting in systematic biases and poor prediction consistency during critical stages such as rapid rises and falls in water levels and peak and valley extremes.

[0005] In summary, to address the aforementioned issues, there is an urgent need for a prediction framework that can combine multi-scale characterization, adaptive input selection, and cross-scale nonlinear fusion capabilities under the constraint of "low data availability". Summary of the Invention

[0006] The purpose of this invention is to overcome the above-mentioned defects of the prior art and provide a multi-scale prediction method for groundwater level in areas with sparse meteorological observations. This method solves the technical problems of high dependence on meteorological driving data, redundancy of lag features, low accuracy of non-stationary sequence modeling, and large deviation of linear reconstruction system in the prior art, and achieves high-precision, high-generalization, and high-robust prediction based solely on historical groundwater level sequences.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A multi-scale method for predicting groundwater levels in areas with sparse meteorological observations includes the following steps: Step S1: Data Preprocessing and Construction of Candidate Lag Feature Library: Quality control and preprocessing are performed on historical 15-day groundwater level (GWL) time-series data. Based on autocorrelation (ACF) and partial autocorrelation (PACF) analysis, the upper limit of the lag order is determined, and a GWL library covering a 1-year hydrological memory scale is constructed. (t-1) ~GWL (t-24) The lag candidate feature library is used, where t is the prediction time, and the training set, validation set, and test set are divided in chronological order. Step S2: HSIC-Lasso Adaptive Lag Feature Selection: The HSIC-Lasso algorithm is used to capture the nonlinear dependency between candidate lag terms and the prediction target through the radial basis function (RBF) kernel. A sparse optimization objective function is constructed and solved. The optimal regularization parameter and the number of features are determined through 5-fold rolling time-series cross-validation. Redundant features are removed to obtain the optimal input feature set. Step S3: Adaptive multi-scale decomposition of groundwater level time series is performed using LPSO-Adaptive-VMD: an adaptive variational mode decomposition framework optimized by linear variational particle swarm optimization (LPSO) is constructed, and the mode number K and bandwidth penalty coefficient α of variational mode decomposition (VMD) are jointly optimized. A comprehensive objective function is constructed based on reconstruction error, mode orthogonality, and entropy complexity to achieve data-driven adaptive decomposition. The intrinsic mode function (IMF) components are obtained by sorting them from low to high center frequency and are divided into low-frequency trend components and medium-to-high-frequency fluctuation components. Step S4: Scale-based differentiated parallel prediction: For low-frequency trend components, an N-BEATS deep learning model is used to model and capture the long-term evolution of groundwater level; for medium- and high-frequency fluctuation components, a bidirectional gated cyclic unit (BiGRU) model is used to model and capture annual cycle fluctuations and short-term nonlinear disturbances; the two models are trained in parallel and output independent prediction results for each IMF component. Step S5: LightGBM Nonlinear Ensemble Reconstruction: The prediction results of each IMF component are used as meta-features. An out-of-sample meta-feature construction strategy is adopted to avoid information leakage in ensemble learning. The data is input into the LightGBM gradient boosting tree model to learn the nonlinear coupling relationship between multi-scale components, complete the adaptive fusion reconstruction, and output the final groundwater level prediction value.

[0008] Furthermore, in step S1, data preprocessing specifically includes: identifying and removing outliers using the 3σ principle; repairing short-term missing data using linear interpolation; and not using smoothing operations across time boundaries; performing 0-1 normalization based on training set statistics, with the normalization parameters only applied to the validation and test sets; and dividing the training, validation, and test sets into a time sequence ratio of 7:1.5:1.5.

[0009] Furthermore, in step S2, the accelerated proximal gradient method is used to solve the sparse optimization objective function, with a convergence rate of O(1 / t). 2 The optimal regularization parameter λ and the number of features n are determined based on the criterion of minimizing the RMSE or maximizing the NSE of the validation set, and features with non-zero weights are retained as the final input feature set.

[0010] Furthermore, in step S3, the search range for the mode number K is 3~10, and the search range for the bandwidth penalty coefficient α is 500~3500; a comprehensive objective function is constructed for the VMD decomposition results of different parameter combinations. The reconstruction error is the mean square error between the superposition of the decomposed components and the original sequence. The mode orthogonality index is used to suppress frequency band overlap, and the entropy complexity index is used to suppress pseudo-modes.

[0011] Furthermore, in step S3, LPSO replaces the fixed inertia term by introducing linear variation weights. The inertia weights decrease linearly with iteration, making the algorithm biased towards global search in the early stage and focused on local convergence in the later stage. IMF1-IMF2 are divided into low-frequency trend components, and IMF3 and above are divided into medium- and high-frequency fluctuation components.

[0012] Furthermore, in step S4, the N-BEATS model is configured as 3 stacked modules, each stack containing 4 fully connected layers with 256 hidden neurons per layer, using the ReLU activation function and the Adam optimizer, with an initial learning rate of 1×10⁻⁶. -3 The batch size is 32, the maximum number of iterations is 200, and an early stopping strategy is introduced to terminate training when the performance on the validation set does not improve for 20 consecutive iterations.

[0013] Furthermore, in step S4, the BiGRU model is configured as a two-layer bidirectional GRU structure with 64 hidden units per layer. During training, a dropout of 0.20 is applied, and the AdamW optimizer is used with an initial learning rate of 5×10⁻⁶. -4 The batch size is 32, and an early stopping strategy is introduced to terminate training when the validation set error has not improved for 15 consecutive training cycles.

[0014] Furthermore, in step S5, the strategy for constructing meta-features for temporal samples is as follows: after training the N-BEATS and BiGRU sub-models on the training set, rolling temporal cross-validation is performed on the training set to generate one-step advance prediction outputs of each sub-model on the corresponding validation fold, and these are concatenated to form a meta-feature matrix.

[0015] Furthermore, in step S5, after training the N-BEATS and BiGRU sub-models on the training set, rolling temporal cross-validation is performed on the training set to generate one-step-ahead prediction outputs for each sub-model on the corresponding validation fold, which are then concatenated to form a meta-feature matrix; the LightGBM model uses a regression mode to construct the ensemble tree.

[0016] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention uses only historical groundwater level sequences as the sole input, without requiring any meteorological driving data such as precipitation or temperature. This solves the core pain point that models cannot be deployed in areas with sparse meteorological observations, and can be directly extended to remote grasslands, deserts, and other areas with scarce data.

[0017] 2. The HSIC-Lasso algorithm is used to adaptively filter nonlinear hysteresis features and eliminate redundant inputs. Compared with traditional partial correlation and partial mutual information methods, it improves NSE, has smaller performance dispersion across well points, and remains stable under different hydrological backgrounds.

[0018] 3. By using LPSO-Adaptive-VMD, the VMD parameters are fully optimized in a joint manner, avoiding the problems of over-decomposition / under-decomposition caused by manual parameter selection. The decomposition structure fits the real hydrological scale characteristics of the groundwater level sequence, and the physical interpretability and decomposition stability are significantly improved.

[0019] 4. To address the different characteristics of low-frequency trends and medium-to-high-frequency fluctuations, we employ adapted N-BEATS and BiGRU models respectively. Compared to a single model, this approach effectively balances the accurate prediction of both long-term trends and short-term disturbances.

[0020] 5. Breaking through the inherent limitations of traditional linear superposition, LightGBM nonlinear ensemble captures cross-scale coupling relationships, significantly improving prediction bias during peak-valley extremes and rapid transition stages. Compared with linear reconstruction, the phase consistency between the predicted sequence and the observed sequence is greatly improved. Attached Figure Description

[0021] Figure 1 This is a diagram illustrating the method steps of the present invention; Figure 2 This is an overview of the study area and a spatial distribution map of monitoring wells and meteorological stations for this invention. Figure 3 This is a comparison chart of the hysteresis feature screening results of three representative grassland monitoring wells of this invention using the HSIC-Lasso, PMI, and PCR methods. Figure 4 The box plot shows the predictive performance comparison of the three feature screening methods of the present invention in 27 monitoring wells, where (a) is the NSE index and (b) is the RMSE index. Figure 5 This is a time series diagram of the LPSO-Adaptive-VMD decomposition results of representative monitoring wells in the Hulunbuir meadow grassland of this invention; Figure 6 This is a comparison chart of the IMF-level prediction performance of the low-frequency components and mid-to-high-frequency components in the three types of grassland areas in this invention. Figure 7 This is a time-series comparison diagram of the predicted and observed groundwater levels reconstructed by linear superposition according to the present invention. Detailed Implementation

[0022] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0023] Numerous 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 therefore the invention is not limited to the specific embodiments disclosed below.

[0024] Example The complete implementation steps of the method of the present invention are as follows: Figure 1-7 : This embodiment uses groundwater level time-series data (15 days each year from 2000 to 2020) from 27 groundwater monitoring wells in three types of grassland areas in the east-west trending zone of Inner Mongolia: Hulunbuir meadow steppe, Ordos desert steppe, and Xilingol typical steppe. The method of this invention is applied to these areas. Overview of the implementation areas: The eastern Hulunbuir region has a cold-temperate continental monsoon climate with an average annual temperature of -3.0 to 1.5℃ and precipitation of 350 to 450 mm. The soil is mainly black calcareous soil and meadow soil, and groundwater is primarily replenished by snowmelt and summer precipitation. The central Xilingol steppe is a temperate semi-arid zone with an temperature of 2.0 to 5.0℃ and precipitation of 250 to 350 mm. Evaporation is much higher than precipitation. The soil is mainly chestnut calcareous soil and sandy soil, and groundwater replenishment mainly comes from precipitation infiltration and local slope runoff. The Ordos desert steppe in western China is located in arid-semi-arid transitional zones, with temperatures ranging from 6.0 to 9.0℃, precipitation from 200 to 300 mm with significant interannual variability, and evaporation reaching 2000 to 2600 mm. Aeolian sandy soils and calcareous gray soils are widespread, and groundwater is primarily replenished by stormwater runoff and alluvial fans at the foot of mountains. All three types of grassland face common challenges, including sparse meteorological observation stations, uneven spatial and temporal distribution of precipitation, and difficulties in obtaining high-resolution meteorological inputs for data-driven models.

[0025] The specific steps of the method proposed in this invention are as follows: Step S1: Data Preprocessing and Construction of Candidate Lag Feature Library: Quality control and preprocessing are performed on historical 15-day groundwater level (GWL) time-series data. Based on autocorrelation (ACF) and partial autocorrelation (PACF) analysis, the upper limit of the lag order is determined, and a GWL library covering a 1-year hydrological memory scale is constructed. (t-1) ~GWL (t-24) A lag candidate feature library is used, where t is the prediction time, and the library is divided into training, validation, and test sets in chronological order; the details are as follows: 1.1 Representative wells for each type of grassland were selected to display the decomposition results and time series comparison diagram, while the remaining wells were used for spatial statistics and generalization evaluation. After quality control, the overall missing rate of the raw data was less than 2%. Outliers were identified and removed using the 3σ principle. Short-term missing segments were repaired using linear interpolation. Smoothing operations across time series boundaries were not used to avoid introducing potential future information and changing the high-frequency structure.

[0026] 1.2 0-1 normalization is used to eliminate the influence of burial depth differences and ensure the comparability of different well points. The minimum and maximum values ​​of normalization are calculated only from the training set data and applied to the validation and test sets with the same parameters.

[0027] 1.3 Based on autocorrelation (ACF) and partial autocorrelation (PACF) analysis, the upper limit of the lag order was determined to be 24, corresponding to 24 15-day periods, covering a hydrological annual cycle of approximately 360 days, and a GWL was constructed. (t-1)~GWL (t-24) The lagging candidate feature library.

[0028] 1.4 HSIC-Lasso is used to filter lag features in the candidate library to obtain the final input lag set T={τ1,…,τn}. All data are arranged in chronological order, with 70% as the training set, 15% as the validation set, and 15% as the test set. Random partitioning is not used to ensure temporal integrity.

[0029] Step S2: HSIC-Lasso Adaptive Lag Feature Selection: The HSIC-Lasso algorithm is used to capture the nonlinear dependency between candidate lag terms and the prediction target through the radial basis function (RBF) kernel. A sparse optimization objective function is constructed and solved. The optimal regularization parameter and the number of features are determined through 5-fold rolling cross-validation. Redundant features are removed to obtain the optimal input feature set; the details are as follows: 2.1 Input and output variables are normalized to 0-1 as in step 1.2; the kernel matrix is ​​centered using a centering matrix during calculation; a radial basis function (RBF) kernel is used. (1) K(x i ,x j The kernel function value between the i-th sample and the j-th sample; σ represents the kernel width parameter, which controls the rate of similarity decay; γ is the kernel parameter.

[0030] Construct the kernel matrix for each candidate hysteresis feature The kernel matrix L of the target variable is centered using a centering matrix.

[0031] 2.2 Optimization Solution: Construct the HSIC-Lasso sparse optimization objective function and solve it quickly using the Accelerated Proximal Gradient (APG) method. In each iteration, the update steps are as follows: (2) Where S λ / L For soft-thresholding, a (t) and a (t+1) This represents the value of the variable at the t-th iteration. f(a) is the Lipschitz constant. (t) () indicates that the objective function f(a) is at the current point a. (t) The gradient at point . The convergence rate of this algorithm is . It exhibits high computational efficiency and stability under sparse high-dimensional sample conditions.

[0032] 2.3 Parameter Tuning and Feature Selection: The weight vector obtained by solving HSIC-Lasso. Based on this, select the option that satisfies... The input variables are used as the selected features; when the number of non-zero features is greater than the expected size, they are selected according to... Sort the features from largest to smallest and retain the top n features. This embodiment uses a 5-fold rolling temporal cross-validation, where each fold is trained using only the earlier time period and validated using the next time period. The optimal regularization parameter λ and the number of features n are determined based on the minimum RMSE (or maximum NSE) of the validation set. Features with non-zero weights are retained as the final input feature set.

[0033] In this embodiment, the number of features retained for the three types of grassland representative well points is 6 to 8, which is significantly lower than the 24 features in the candidate library, and the redundancy is greatly reduced.

[0034] Step S3: Adaptive multi-scale decomposition of groundwater level time series is performed using LPSO-Adaptive-VMD to obtain IMF components sorted by center frequency, and divided into low-frequency trend components and mid-to-high-frequency components according to frequency level. Groundwater level sequences exhibit typical seasonal variations and multi-scale recharge response superposition properties. The coexistence of long-term trends and high-frequency disturbances makes them significantly non-stationary. Traditional Variational Mode Decomposition (VMD) can decompose a signal into a set of intrinsic mode functions (IMFs) with finite bandwidth, but its decomposition quality is strongly dependent on the number of modes. K With bandwidth penalty coefficient α Manual specification is crucial. Improper parameter settings can easily lead to over-decomposition (modal redundancy) or under-decomposition (overlapping trends and disturbances), thus weakening the effectiveness of subsequent prediction models. Therefore, this invention constructs an adaptive parameter optimization VMD framework based on particle swarm optimization (LPSO-Adaptive-VMD). This framework utilizes a multi-index objective function to drive automatic parameter search, ensuring that the decomposition structure is determined by data characteristics rather than manual setting, thereby enhancing the stability and physical interpretability of multi-scale groundwater level analysis.

[0035] The specific steps are as follows: 3.1 Parameter Search Space Setting: To reduce the subjectivity of parameter setting, this invention introduces Linear-Mutation Particle Swarm Optimization (LPSO) to achieve joint optimization of K and α using two variables; the search range for the mode number K is 3~10, and the search range for the bandwidth penalty coefficient α is 500~3500, which have been adjusted through preliminary experiments; each particle is encoded as follows: This is used to generate different VMD decomposition candidate structures.

[0036] 3.2 Construction of the comprehensive objective function; A comprehensive objective function is constructed for the VMD decomposition results with different parameter combinations. The reconstruction error is the mean square error between the superposition of the decomposed components and the original sequence. The modal orthogonality index is used to suppress frequency band overlap, and the entropy complexity index is used to suppress pseudo-modes. (3) in: Reconstruction error: (4) Modal orthogonality (suppressing overlap): (5) Entropy complexity (suppressing pseudo-modalities): (6) Weight ω 1, ω 2, ω The model is automatically tuned based on its performance on the validation set.

[0037] 3.3 LPSO Optimization: LPSO avoids premature convergence of standard PSO by introducing linearly variable weights to replace the fixed inertia term, thus improving optimization stability and speed. The particle velocity and position are updated as follows: (7) (8) The inertia weight decreases linearly with iteration: (9) This approach initially favors global search, then focuses on local convergence in the later stages. The final optimal parameters are: (10) 3.4 Decomposition and Component Division: After obtaining the VMD decomposition results under the optimal parameters, the IMF components are sorted from low to high center frequency. IMF1-IMF2 are divided into low-frequency trend components (LOW), and IMF3 and above are divided into medium- and high-frequency fluctuation components (HIGH).

[0038] Step S4: Scale-based differentiated parallel prediction: For low-frequency trend components, an N-BEATS deep learning model is used to model and capture the long-term evolution of groundwater level; for medium- and high-frequency fluctuation components, a bidirectional gated cyclic unit (BiGRU) model is used to model and capture annual cycle fluctuations and short-term nonlinear disturbances; the two models are trained in parallel and output independent prediction results for each IMF component. 4.1 Low-frequency component N-BEATS prediction: The model is configured with 3 stacked modules, each stack containing 4 fully connected layers with 256 hidden neurons per layer, and ReLU is used as the activation function to enhance the model's non-linear expressive power; the Adam optimizer is used for model training with an initial learning rate of 1×10⁻⁶. -3 The batch size is 32, the maximum number of iterations is 200, and an early stopping strategy is introduced to prevent the model from overfitting. Training is terminated when the performance on the validation set does not improve for 20 consecutive iterations, and the prediction results of the low-frequency components are output.

[0039] 4.2 BiGRU Prediction for Mid-to-High Frequency Components: The model is configured as a two-layer bidirectional GRU structure with 64 hidden units per layer. A dropout of 0.20 is added during training to enhance the network's generalization ability. The BiGRU output is mapped to the corresponding future IMF prediction value through a fully connected layer. The AdamW optimizer is used for training with an initial learning rate of 5×10⁻⁶. -4 The batch size is 32. To avoid overfitting, an early stopping mechanism is also introduced. When the validation set error does not improve within 15 consecutive training cycles, the training is terminated early, and the prediction results of each mid-to-high frequency IMF component are output.

[0040] The results of steps 3 and 4 in this embodiment are analyzed as follows: The optimal decomposition order of the three representative well points of the three types of grasslands differs significantly: the optimal K of the three representative monitoring wells of Hulunbuir meadow grassland, Ordos desert grassland and Xilingol typical grassland are 9, 11 and 14 respectively. This difference reflects that the energy spectrum distribution and resolvable time scale of the groundwater level sequence in the three places are significantly different: as the hydrothermal gradient transitions from humid to arid, intra-annual fluctuations and short-term disturbances account for a higher proportion and are more abundant in the sequence, thus requiring more modes to characterize its multi-scale structure.

[0041] In the optimal decomposition results, each IMF is sorted from low to high according to its center frequency, exhibiting a clear scale hierarchy: IMF1-IMF2 are multi-year scale gradual components, characterizing the long-term background evolution and interannual fluctuations of groundwater level; IMF3 reflects significant annual cycle oscillations; higher-order IMFs (IMF4+) mainly correspond to rapid fluctuations within the year to shorter time scales, and exhibit amplitude modulation and stage-specific enhancement characteristics in some stages. Based on this center frequency structure, this paper defines IMF1-IMF2 as low-frequency components (LOW), and IMF3 and above as high-frequency components (HIGH), and uses N-BEATS and BiGRU for subscale prediction respectively.

[0042] The IMF-level predictive performance assessment results show that, overall, the low-frequency component exhibits significantly higher and more stable predictability in all three monitoring wells: the average NSE for LOW (IMF1-IMF2) is 0.765-0.766, and the average RMSE is 0.124-0.138m. In contrast, the prediction difficulty of the high-frequency component increases significantly: the average NSE for HIGH (IMF3+) decreases to 0.669-0.690, and the average RMSE increases to 0.187-0.203m. Furthermore, with the increase of the IMF order, there is an overall degenerate trend of decreasing NSE and increasing RMSE, reflecting that the short-term disturbance components are more random and have weaker time persistence. In terms of regional comparison, the performance degradation of high-order IMFs in the Ordos desert steppe is the most obvious (high-order IMF NSE close to 0.64, RMSE up to 0.24 m), consistent with its strong evapotranspiration constraints, unstable precipitation timing, and more pulsed recharge events, making high-frequency fluctuations more difficult to predict stably. The typical grassland in Xilingol is affected by both precipitation fluctuations and freeze-thaw processes, resulting in a larger number of high-frequency modes (K=14), and its high-frequency error exhibits a more complex structure as the frequency changes. In contrast, in the Hulunbuir meadow steppe, snowmelt and relatively regular seasonal precipitation give water level fluctuations stronger continuity, so the high-frequency components remain at a relatively high level overall.

[0043] Based on the linear superposition and reconstruction of the prediction results from each IMF sub-prediction, the time series prediction effect of groundwater levels from the three monitoring wells is as follows: Figure 7As shown in the figure, the two vertical dashed lines divide the samples into training, validation, and test sets; the small window in the upper right corner is a magnified comparison of the entire test period, used to more clearly demonstrate the time-series tracking capability of the extrapolation stage. Overall, the linear superposition of the prediction results of each IMF sub-prediction can reproduce the long-term evolution trend and dominant cyclical fluctuations of groundwater level relatively well. The reconstruction accuracy of groundwater in the three types of grasslands during the test period is as follows: Hulunbuir meadow grassland NSE=0.766, RMSE=0.227m; Ordos desert grassland NSE=0.752, RMSE=0.254m; Xilingol typical grassland NSE=0.761, RMSE=0.135m. From the magnified results of the test period, it can be seen that during the rapid decline / rise of water level and the local peak-valley change stages, there are still identifiable deviations between the predicted curve and the observed curve, and these deviations are more obvious in the Ordos and Xilingol regions where high-frequency disturbances are stronger.

[0044] Step 5: LightGBM Nonlinear Integration Reconstruction: A nonlinear reconstruction model is constructed using Light Gradient Boosting Machine (LightGBM), as shown in Equation (11). All IMF prediction results are used as input features. The nonlinear combination relationship between multi-scale components is learned through gradient boosting tree, thereby generating the final groundwater level prediction. (11) in: The predicted value of the k-th IMF; : The nonlinear mapping function learned by LightGBM.

[0045] The specific steps are as follows: 5.1 Training LightGBM using out-of-sample meta-feature construction strategy: After training the N-BEATS and BiGRU sub-models on the training set, rolling temporal cross-validation is performed on the training set to generate one-step-ahead prediction outputs of each sub-model on the corresponding validation fold. These are then concatenated to form a meta-feature matrix, avoiding overfitting and information leakage caused by "training both sub-models and ensemblers on the same data".

[0046] 5.2 LightGBM was trained using the out-of-sample meta-features and the real GWL as supervisory signals. LightGBM used a regression mode to construct the ensemble tree. The core parameters included: a learning rate of 0.05 to balance the model's convergence speed and generalization ability; a maximum tree depth of 6 to avoid overfitting; 31 leaf nodes to enhance the model's expressive power; sampling rates (bagging_fraction and feature_fraction) of 0.8 to improve model stability; and 500 iteration trees to obtain sufficient fitting ability.

[0047] 5.3 Prediction Output: With the sub-model and integrator parameters fixed, the test set is input into the model. After feature selection, decomposition, scale prediction, and nonlinear fusion, the final predicted groundwater level is output and reverse-normalized to the original dimension (m) to complete the final prediction.

[0048] The final test set performance of this embodiment: the overall NSE of 27 monitoring wells was 0.871~0.925, and the RMSE was 0.075~0.159m. Among them, the average NSE of the typical grassland in Xilingol reached 0.909, and the average RMSE was as low as 0.089m. The prediction accuracy and consistency are significantly better than existing methods.

[0049] Comparative example: This comparative study employed partial mutual information (PMI) and partial correlation analysis (PCR) to screen candidate lag terms GWL(t-1) to GWL(t-24) from 27 groundwater monitoring wells across three grassland types: Hulunbuir meadow steppe, Ordos desert steppe, and Xilingol typical steppe. Figure 3 This paper presents the screening results of representative monitoring wells in three types of grassland areas using the two methods in this comparative example and the HSIC-Lasso method in the examples, where n represents the number of hysteresis terms ultimately retained. The results show that HSIC-Lasso produces a relatively sparse hysteresis structure in all three types of grassland areas. Figure 3 (n=6-8), and the selected lag distribution exhibits obvious multi-scale characteristics. Hulunbuir mainly shows short-term lags (t-1, t-2), accompanied by several medium- and longer lags (such as t-4, t-7, t-11, t-20), indicating that groundwater level changes are simultaneously influenced by rapid response and seasonal-scale memory. Ordos retains more dispersed lags, including not only short-term lags but also multiple medium-term lags (n=6-8), and the selected lags show obvious multi-scale characteristics. Figure 3 (n=8), reflecting the delayed and slow-release characteristics of the groundwater system response in this area. Xilingol also retains both short-term and long-term lag terms ( Figure 3(n=8), showing the characteristics of short-term response superposition with seasonal-cross-seasonal scale signals. Overall, the hysteresis system obtained by HSIC-Lasso is more sparse and structurally clear, providing input with low redundancy for subsequent prediction models.

[0050] In contrast, the PMI typically retains more lagged terms in the three types of grassland regions ( Figure 3 (n=8-11), and the inclusion lags often appear in clusters in adjacent or nearby lags, suggesting that this method is more likely to retain lag variables with similar time scales simultaneously, thus introducing some redundancy. PCR screening results are more inclined towards linear responses, mainly with short-term lags (t-1, t-2), supplemented by some mid-term and longer-term lags (t-1, t-2). Figure 3 (n=10-11), which has relatively limited characterization of nonlinear dependence and cross-seasonal lag structure.

[0051] The accuracy of the BiGRU prediction model built based on the screening results of three types of variables is as follows: Figure 4 As shown, the three methods exhibited consistent relative performance rankings across different grassland regions: HSIC-Lasso was the best, followed by PMI, and PCR was the weakest. In the Hulunbuir meadow steppe, groundwater primarily responds to short-period fluctuations in spring snowmelt pulses. The short-period lag set identified by HSIC-Lasso resulted in the highest predictive performance (NSE≈0.644, RMSE≈0.365), outperforming PMI (NSE≈0.627) and PCR (NSE≈0.616). In the Ordos desert steppe, hydrological processes are controlled by evapotranspiration, and the "delayed infiltration-deep slow release" mechanism is more pronounced. HSIC-Lasso accurately identified seasonal-mid-season characteristics after GWL(t-10), better characterizing the "slow replenishment" process of the desert steppe compared to the features screened by PMI and PCR (NSE≈0.636). The typical grassland in Xilin Gol is affected by both summer precipitation and spring and autumn freeze-thaw cycles, resulting in a complex lag structure with superimposed short-term and seasonal scales. The "bimodal characteristic" of HSIC-Lasso matches this mechanism perfectly, thus achieving the highest accuracy in this region (NSE≈0.650).

[0052] Looking at the differences between different methods, HSIC-Lasso showed an average improvement in NSE of approximately 0.016-0.034 compared to PMI, and an improvement of approximately 0.027-0.047 compared to PCR. In terms of RMSE, HSIC-Lasso showed a reduction of approximately 0.001-0.029 compared to PMI, and a reduction of approximately 0.011-0.041 compared to PCR. These results were consistent across the three types of grassland areas, indicating that the input features constructed by HSIC-Lasso have good stability and generalization potential under different hydrological backgrounds. Figure 4The distribution characteristics shown indicate that HSIC-Lasso not only performs better at the mean level, but also has relatively smaller performance dispersion at the well point scale, suggesting that this method can more stably support the construction of groundwater level prediction models under different monitoring conditions.

[0053] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the claimed invention. The scope of protection claimed by the appended claims and their equivalents is defined.

Claims

1. A multi-scale prediction method for groundwater levels in areas with sparse meteorological observations, characterized in that, Includes the following steps: Step S1: Data Preprocessing and Construction of Candidate Lag Feature Library: Quality control and preprocessing are performed on historical 15-day groundwater level (GWL) time-series data. Based on autocorrelation (ACF) and partial autocorrelation (PACF) analysis, the upper limit of the lag order is determined, and a GWL library covering a 1-year hydrological memory scale is constructed. (t-1) ~GWL (t-24) The lag candidate feature library is used, where t is the prediction time, and the training set, validation set, and test set are divided in chronological order. Step S2: HSIC-Lasso Adaptive Lag Feature Selection: The HSIC-Lasso algorithm is used to capture the nonlinear dependency between candidate lag terms and the prediction target through the radial basis function (RBF) kernel. A sparse optimization objective function is constructed and solved. The optimal regularization parameter and the number of features are determined through 5-fold rolling time-series cross-validation. Redundant features are removed to obtain the optimal input feature set. Step S3: Adaptive multi-scale decomposition of groundwater level time series is performed using LPSO-Adaptive-VMD: an adaptive variational mode decomposition framework optimized by linear variational particle swarm optimization (LPSO) is constructed, and the mode number K and bandwidth penalty coefficient α of variational mode decomposition (VMD) are jointly optimized. A comprehensive objective function is constructed based on reconstruction error, mode orthogonality, and entropy complexity to achieve data-driven adaptive decomposition. The intrinsic mode function (IMF) components are obtained by sorting them from low to high center frequency and are divided into low-frequency trend components and medium-to-high-frequency fluctuation components. Step S4: Scale-based differentiated parallel prediction: For low-frequency trend components, an N-BEATS deep learning model is used to model and capture the long-term evolution of groundwater levels; for medium- and high-frequency fluctuation components, a bidirectional gated cyclic unit (BiGRU) model is used to model and capture annual cycle fluctuations and short-term nonlinear disturbances. The two models are trained in parallel, and each model outputs an independent prediction result for each IMF component. Step S5: LightGBM Nonlinear Ensemble Reconstruction: The prediction results of each IMF component are used as meta-features. An out-of-sample meta-feature construction strategy is adopted to avoid information leakage in ensemble learning. The data is input into the LightGBM gradient boosting tree model to learn the nonlinear coupling relationship between multi-scale components, complete the adaptive fusion reconstruction, and output the final groundwater level prediction value.

2. The multi-scale prediction method for groundwater levels in areas with sparse meteorological observations according to claim 1, characterized in that, In step S1, data preprocessing specifically includes: identifying and removing outliers using the 3σ principle; repairing short-term missing data using linear interpolation; and not using smoothing operations across time boundaries; performing 0-1 normalization based on training set statistics, with the normalization parameters only applied to the validation and test sets; and dividing the training, validation, and test sets into a time sequence ratio of 7:1.5:1.

5.

3. The multi-scale groundwater level prediction method for areas with sparse meteorological observations as described in claim 1, characterized in that, In step S2, the accelerated proximal gradient method is used to solve the sparse optimization objective function, with a convergence rate of O(1 / t). 2 The optimal regularization parameter λ and the number of features n are determined based on the criterion of minimizing the RMSE or maximizing the NSE of the validation set, and features with non-zero weights are retained as the final input feature set. 4.The groundwater level multi-scale prediction method for meteorological observation sparse areas according to claim 3, characterized in that, In step S3, the search range for the mode number K is 3~10, and the search range for the bandwidth penalty coefficient α is 500~3500; for the VMD decomposition results of different parameter combinations, a comprehensive objective function is constructed. The reconstruction error is the mean square error between the superposition of the decomposed components and the original sequence. The mode orthogonality index is used to suppress frequency band overlap, and the entropy complexity index is used to suppress pseudo-modes.

5. The multi-scale prediction method for groundwater levels in areas with sparse meteorological observations according to claim 1, characterized in that, In step S3, LPSO replaces the fixed inertia term with a linear variation weight. The inertia weight decreases linearly with iteration, making the algorithm biased towards global search in the early stage and focused on local convergence in the later stage. IMF1-IMF2 are divided into low-frequency trend components, and IMF3 and above are divided into medium- and high-frequency fluctuation components.

6. The multi-scale prediction method for groundwater levels in areas with sparse meteorological observations according to claim 1, characterized in that, In step S4, the N-BEATS model is configured as 3 stacked modules, each stack containing 4 fully connected layers with 256 hidden neurons per layer, using the ReLU activation function and the Adam optimizer, with an initial learning rate of 1×10⁻⁶. -3 The batch size is 32, the maximum number of iterations is 200, and an early stopping strategy is introduced to terminate training when the performance on the validation set has not improved for 20 consecutive iterations. 7.The groundwater level multi-scale prediction method for meteorological observation sparse areas according to claim 1, characterized in that, In step S4, the BiGRU model is configured as a two-layer bidirectional GRU structure with 64 hidden units per layer. During training, a dropout of 0.20 is applied, and the AdamW optimizer is used with an initial learning rate of 5×10⁻⁶. -4 The batch size is 32, and an early stopping strategy is introduced to terminate training when the validation set error has not improved for 15 consecutive training cycles. 8.The groundwater level multi-scale prediction method for meteorological observation sparse areas according to claim 1, characterized in that, In step S5, the strategy for constructing meta-features for temporal samples is as follows: after training the N-BEATS and BiGRU sub-models on the training set, rolling temporal cross-validation is performed on the training set to generate one-step advance prediction outputs of each sub-model on the corresponding validation fold, and these are concatenated to form a meta-feature matrix. 9.The groundwater level multi-scale prediction method for meteorological observation sparse areas according to claim 1, characterized in that, In step S5, after training the N-BEATS and BiGRU sub-models on the training set, rolling temporal cross-validation is performed on the training set to generate one-step-ahead prediction outputs of each sub-model on the corresponding validation fold, which are then concatenated to form a meta-feature matrix; the LightGBM model uses a regression mode to construct an ensemble tree.

Images