A CNN-LSTM Model Based on a Meta-Learning Algorithm to Predict Groundwater Level in the Middle and Lower Reaches of the Heihe River, China

Abstract

In this study, a deep learning model is proposed to predict groundwater levels. The model is able to accurately complete the prediction task even when the data utilized are insufficient. The hybrid model that we have developed, CNN-LSTM-ML, uses a combined network structure of convolutional neural networks (CNN) and long short-term memory (LSTM) network to extract the time dependence of groundwater level on meteorological factors, and uses a meta-learning algorithm framework to ensure the network’s performance under sample conditions. The study predicts groundwater levels from 66 observation wells in the middle and lower reaches of the Heihe River in arid regions and compares them with other data-driven models. Experiments show that the CNN-LSTM-ML model outperforms other models in terms of prediction accuracy in both the short term (1 month) and long term (12 months). Under the condition that the training data are reduced by 50%, the MAE of the proposed model is 33.6% lower than that of LSTM. The results of ablation experiments show that CNN-LSTM-ML is 26.5% better than the RMSE of the original CNN-LSTM structure. The model provides an effective method for groundwater level prediction and contributes to the sustainable management of water resources in arid regions.

1. Introduction

Groundwater is the world’s largest distributed reservoir of fresh water [1], and bears an important responsibility to meet the daily water consumption of 2.5 billion people around the world [2]. In addition, more than half of the world’s food irrigation water comes from groundwater resources [3]. In arid and semi-arid regions with limited rainfall, groundwater resources are especially important for economic and agricultural development [4,5]. In recent years, due to climate change, geological change, and increasing population, the global demand for groundwater resources has also increased sharply [6,7,8]. At the same time, groundwater reserves have rapidly decreased [9], while policies to manage groundwater have been neglected. The reduction of groundwater resources will bring a series of problems, including reduced food production and environmental degradation [10,11]. Therefore, the accurate prediction of groundwater level is crucial for relevant government departments to protect the sustainable development of groundwater resources [12].

In fact, the depth of groundwater level changes over time and is affected by various factors. These factors include meteorological variables, such as precipitation, evapotranspiration, and snow melt [13,14,15,16,17,18,19], as well as surface water around groundwater monitoring sites [20], geological hazards [21], and human factors [22,23]. Therefore, predicting the groundwater level is a highly dynamic and complex problem.

A model that can accurately predict the groundwater level needs to consider time series information on the existing level and the factors that can change the level. Groundwater models based on modern geophysical methods are used to study certain hydrogeological properties of aquifers through electrical conductivity [24,25,26]. These methods are expensive because they usually require manual surveying. With the development of satellite technology and measurement tools, technologies based on remote sensing (RS) and geographic information systems (GIS) have played a huge role in groundwater research [27,28,29]. These methods have obvious advantages in collecting data and drawing graphics and provide technical support for the monitoring of groundwater levels. However, due to the real-time nature of RS and GIS methods, their predictions of groundwater level rely on human judgment, with imperfect results.

By contrast, data-driven models, driven by data sensor technology, have played a greater role in water resource management [30,31,32]. As one of the classic methods for groundwater level prediction, the multiple linear regression method is widely used [33]. The prediction process takes into account the influence of various factors such as rainfall, temperature and river water level on the groundwater level [34], and the constructed model has a high degree of interpretability. However, the multiple regression model performs poorly in predicting the groundwater level under the condition of a small sample size. In addition, if there are too many input variables, problems such as multicollinearity are prone to occur, which will interfere the performance of the model. An autoregressive integrated moving average (ARIMA) model, as a classic time series forecasting method, performs well in the short-term forecasting of groundwater levels [35,36,37]. However, as with multiple regression, the predictive performance of ARIMA is greatly reduced if the amount of data available for modeling is small. In addition, the process of setting appropriate parameters for ARIMA is complicated and its long-term predictions can lack specificity.

With the rapid development of computer technology, deep learning models represented by a large number of neurons have begun to be used in the field of groundwater level prediction. These include multilayer perceptron (MLP) models, also known as artificial neural networks (ANNs). Because these models have high predictive accuracy and do not need to clarify the physical relationship among variables, they have many applications in groundwater prediction [38,39,40]. Recurrent neural networks (RNNs) have natural advantages in time series forecasting and are therefore used for groundwater level forecasting [41]. However, because RNNs are prone to gradient disappearance and gradient explosion problems, a long short-term memory (LSTM) network is usually used to replace RNN. LSTM fully considers the temporal dependence of the input sequence, and has outstanding performance in extracting temporal information. Therefore, many models use LSTM for groundwater level prediction [42,43,44,45,46]. Although LSTM has obvious advantages in long time series prediction, it has insufficient ability to extract the characteristics of influencing factors, so it has certain limitations in extracting the characteristics of groundwater dynamic changes under multiple influences. Convolutional neural networks (CNNs) shine in the field of deep learning due to their powerful feature extraction capabilities. Naturally, there are also studies applying CNN to groundwater level prediction modeling [47,48]. However, after all, CNN is a neighborhood feature extractor, which will have some impact on the predictive accuracy of long-term series. A long- and short-term time-series network (LSTNet) [49] is suitable for multivariate time series forecasting, and CNN and LSTM are also added to the network structure, which is better for time series forecasting with periodic changes. Correspondingly, when the amount of data is small, or the periodicity of sample data is weak, the universality is poor.

On the one hand, the groundwater level prediction model based on deep learning requires large-scale training samples to train the model [50,51], and one of the problems of groundwater prediction is the lack of sufficient data [52], so it is difficult to accumulate enough training samples for the convergence of deep learning model parameters. On the other hand, the existing groundwater level prediction models are weak in modeling the multivariate influencing factors at different times, and most of them cannot simulate the dynamic changes of the groundwater level in time. As an emerging method in the field of deep learning, meta-learning has shown strong performance in model generalization ability, self-learning ability, and few-shot classification and regression problems [53,54,55]. Meta-learning can handle problems with fewer training samples, but currently lacks a network structure that matches the target task, which coincides with the purpose of this study. The model-agnostic meta-learning algorithm (MAML) [56] can finally determine the initial parameters of the model for the target prediction task by meta-training the training samples. This means that the model can learn by itself even under the condition of insufficient sample size, thereby quickly adapting to new prediction tasks. In the end, better performance can still be obtained under the condition of a small number of target samples.

For the above problems and backgrounds, this study proposes a deep learning hybrid model, convolutional neural networks and long short-term memory networks based on meta-learning (CNN-LSTM-ML), to accurately predict groundwater levels. The main contributions of this study are as follows:

(1)

Aiming to resolve the problem that the existing methods are not effective for long-term prediction of groundwater level, and the insufficient modeling ability of various influencing factors. In this paper, a deep learning network structure of hybrid CNN-LSTM is designed, in which the CNN module can effectively extract the multivariate features that affect the groundwater level, and the LSTM module has natural advantages for long-term time series prediction. Therefore, the network structure can effectively solve the above problems.

(2)

For deep learning models, large-scale training data are usually needed as a support, and the real groundwater prediction problem usually lacks sufficient samples. In this study, a meta-learning algorithm is added to the CNN-LSTM network structure, so that the model can train a meta-learner with fewer training samples to complete the task of groundwater level prediction. To the authors’ knowledge, there is no single paper that applies meta-learning to groundwater level prediction, and our study provides support for the expansion of meta-learning algorithm applications.

(3)

To verify the performance of the CNN-LSTM-ML, this study is conducted on a real groundwater level dataset in the middle and lower reaches of the Heihe River. Experimental results show that in short-term prediction (1 month), the MAE of CNN-LSTM-ML is 11.7% lower than that of the multiple regression method. In long-term prediction (12 months), the RMSE of CNN-LSTM-ML is 5% lower than that of LSTM. At the same time, the model can still maintain a high prediction accuracy even when the training samples are reduced. All in all, the proposed model can accurately predict groundwater levels, which can help relevant government departments manage water resources and make evidence-based decisions.

2. Materials and Methods

2.1. Groundwater Level Prediction Process

In this study, the groundwater level prediction process is divided into three main stages: data processing, modeling prediction and model evaluation. The specific flow chart is shown in Figure 1.

Data processing is conducted to identify outliers and impute missing data, to ensure the data integrity of the input model. Data normalization is carried out to eliminate the influence of different dimensions on model predictions. To build the model, we use the meta-learning framework and add the CNN and LSTM networks to extract the spatial features and temporal dependencies of the data. The specific details of the model will be explained in the following sections. In the final experiment, scientific evaluation indicators are used to evaluate the model.

2.2. Study Area and Data Processing

2.2.1. Study Area and Data

Located in an arid region of northwestern China, the Heihe is the second-largest inland river in China and is representative of all inland rivers [57]. The Heihe River is located within 37.7°–42.7° N, 97.1°–102.0° E. It originates in the Qilian Mountains of Qinghai Province and flows northward through Gansu Province to the Inner Mongolia Autonomous Region, with a total length of about 821 km and a total area of about 1,432,000 km². There are glaciers, permafrost, forests, and deserts in the Heihe River Basin (HRB) [58], which affect the production and lives of millions of people in these arid areas.

As shown in Figure 2, this study is carried out in the middle and lower reaches of the HRB, which covers Zhangye, Linze, Gaotai, Jinta in Gansu Province, and EJinaqi in the Inner Mongolia Autonomous Region. The simulated water level data of 66 observation wells from 1 January 2003 to 1 December 2012 were selected in the above five study areas [59]. The data are all monthly data. Seven meteorological factors, including evaporation, precipitation, air pressure, relative humidity, sunshine time, temperature and wind speed in the study area, are selected for the study. Historical meteorological data are collected from local meteorological observation sites, supported by local meteorological bureaus. This paper uses the daily data from 1 January 2003 to 31 December 2012. Many of the dates have missing data that require post-processing.

2.2.2. Missing Data Processing

Due to failure of data acquisition equipment, unstable data transmission networks, failure of storage equipment, and possible human factors, the data used for this type of research are usually incomplete. Missing data can lead to biased parameter estimates, missing key location information, increased standard errors, and poor generalizability of the results [60]. Due to the small amount of data collected in this experiment, the main focus of this study is the prediction problem under the condition of few samples, so the existence of missing values will have a great impact on the model in the experiment and the comparison model. In the process of dealing with missing values in time series, methods such as mean imputation, linear imputation, multiple imputation, and k-nearest neighbor (KNN) imputation are usually used [61,62,63].

For the missing values in the meteorological data collected in this study, we firstly obtained data from five national-level surface meteorological observation stations provided by the China Meteorological Data Service Centre (CMDC) (Zhangye 39°05′ N, 100°17′ E; Gaotai 39°22′ N, 99°50′ E; Linze 39°09′ N, 100°10′ E; Jinta 40°00′ N, 98°53′ E; EJinaqi 41°57′ N, 101°04′ E). If there is still a single missing data, we will adopt the method of mean interpolation. The specific implementation is shown in the following formula:

$$x_{NA}=\frac{x_{t-1}+x_{t+1}}{2}$$

(1)

It is assumed that there is a missing value $x_{NA}$ at the current time $t$, and that $x_{t-1}$ and $x_{t+1}$ represent the complete data at the previous time and the time after the missing data, respectively. For the case where multiple consecutive data are missing in the time series, we use linear interpolation to solve this problem.

2.2.3. Outlier Data Processing

We use the boxplot method to analyze outliers. QL is the lower quartile, meaning that 25% of all observations is less than QL. QU is the upper quartile, meaning that one in four observations is greater than QU. IQR is the interquartile range, which is the difference between the upper quartile QU and the lower quartile QL, that is, IQR = QU − QL. Outliers are usually defined as values less than QL − 1.5IQR or greater than QU + 1.5IQR [64]. Outliers can be treated as missing values.

Due to the large number of observation wells, we select one observation well from each of the Zhangye and Jinta study areas for demonstration. The results of boxplot analysis of outliers are shown in Figure 3. Among them, there is an abnormal value point in an observation well in the Zhangye research area, and it is necessary to analyze and de-cide whether to eliminate the interpolation or keep it. In addition, since the original meteorological data is daily data, but the groundwater level data is monthly data, the monthly average of the meteorological data is calculated. After the above processing, data that can be used for experiments are obtained.

2.2.4. Data Normalization

Before the data are put into some of the deep learning models used in this article, the data need to be normalized. The purpose is to eliminate the influence on the model of the feature attribute dimension in the data, so as to accelerate the model convergence and improve the model’s predictive accuracy. In the experiments, we adopt the max-min normalization method:

$$x^{\prime }=\frac{x-x_{min}}{x_{max}-x_{min}}$$

(2)

where $x^{\prime }$ is the normalized data, $x$ represents the current observation data, and $x_{max}$ and $x_{min}$ represent the maximum and minimum values in the current observation data, respectively.

2.3. Convolutional Neural Network (CNN)

CNN is a type of neural network that includes convolutional computations and is specialized for processing similar gridded data. The concept of convolutional networks was first proposed in 1989 to solve computer vision problems [65]. Due to the weight sharing and local connection characteristics of convolutional neural networks, they have made great contributions in deep learning and have been widely used in various fields [66,67,68]. The main purpose of the convolutional neural network is to extract the feature information of the input, which is mainly done by the convolutional layer. Since the feature information extracted by the convolutional layer is linear, the input data usually has nonlinear characteristics. Therefore, it is necessary to introduce a nonlinear function to solve this problem. This nonlinear function is the activation function in the convolutional neural network. In the specific application of machine learning, the input data is usually not a one-dimensional array. Here, two-dimensional discrete convolution is used as an example to show the convolution operation, and the activation function uses the Relu function as an example. The specific operation is as formula:

$$H\left(i,j\right)=\left(I\ast K\right)\left(i,j\right)={\displaystyle \sum }_m{\displaystyle \sum }_{n}K\left(m,n\right)I\left(i+m,j+n\right)$$

(3)

$$f\left(x\right)=max\left(0,x\right)$$

(4)

where $H\left(i,j\right)$ represents a specific position in the feature map after convolution. $I$ and $K\left(m,n\right)$ represent the size of the input array and the convolution kernel, respectively. $x$ is the feature map after the convolution operation, which is usually a tensor. According to the above formula, it can be seen that the convolution kernel can slide on the input array. Each time it slides to a position, its corresponding array elements are multiplied and summed. Finally, after completing the sliding input array, the feature map matrix will be obtained. Each element on the feature map represents the feature information extracted by the convolution operation. Finally, the feature map is put into the activation function to give it nonlinear characteristics, which strengthens the expressive ability of the entire network and plays a key role in data fitting.

Because this research involves many meteorological factors that affect the groundwater level, the primary challenge of our work involves the feature engineering of multivariate time series. Therefore, it is natural to consider using convolutional neural networks for feature extraction. The entire feature engineering process is shown in Figure 4.

The collected multivariate time series is regarded as a two-dimensional array, the corresponding filters are used for feature extraction and integration, and the activation function is used for nonlinear fitting. The final feature map stores the feature information that affects the groundwater level.

2.4. Long-Short Term Memory (LSTM) Network

LSTM is a neural network for learning long-term dependencies from time series data. LSTM was first proposed in 1997 [69], and currently has many applications in the field of time series forecasting [70,71,72]. In fact, LSTM is a variant of the recurrent neural network, which solves the problem of gradient disappearance in the traditional RNN training process, and thus can preserve long-term information. The internal structure of LSTM is shown in Figure 5.

LSTM is a special RNN. The biggest difference from traditional RNN is that it introduces the concept of memory cells. Due to the existence of memory cells, LSTM has the ability to memorize long-term information in the input sequence. The representation of memory cells in Figure 4 is $C_t$. We can see that the state of the memory cell can directly pass through the cyclic chain, where the LSTM chooses to update or forget the cell state through various gates. Specifically, there are three types of LSTM inputs at time $t$, namely the input $x_t$ at the current time, the hidden layer state $h_{t-1}$ and the cell state $C_{t-1}$ at the previous time. The forget gate $f_t$ and the input gate $i_t$ are used to control which information should be added to the cell state, and the output gate $O_t$ is used to calculate the current hidden layer state $h_t$. The final LSTM outputs are the current cell state $C_t$ and the hidden layer state $h_t$. The formula of the specific gating unit is as follows:

$$f_t=\sigma \left(W_f\left[h_{t-1},X_t\right]+b_f\right)$$

(5)

$$i_t=\sigma \left(W_i\left[h_{t-1},X_t\right]+b_i\right)$$

(6)

$$O_t=\sigma \left(W_o\left[h_{t-1},X_t\right]+b_o\right)$$

(7)

where $\left[h_{t-1},X_t\right]$ indicates that the hidden layer and the sequence data are jointly input to the network. $W_f$, $W_i$, $W_o$ and $b_f$, $b_i$, $b_o$ are the corresponding parameters in the network. Specifically, $W_f$, $W_i$, and $W_o$ are the weighted coefficient matrices corresponding to each gate, and $b_f$, $b_i$, and $b_o$ correspond to the offsets. $\sigma$ is the sigmoid activation function of each gate. The process of updating the cell state is as follows:

$${\tilde{C}}_t=\mathit{tanh}(W_c\left[h_{t-1},X_t\right]+b_c)$$

(8)

$$C_t=f_t\ast C_{t-1}+i_t\ast {\tilde{C}}_t$$

(9)

where $\ast$ represents the multiplication of corresponding elements. $\mathit{tanh}$ is the hyperbolic tangent activation function. The updated cell $C_t$ consists of the remaining value in the cell $C_{t-1}$ filtered by the forget gate and the value of the newly added cell ${\tilde{C}}_t$ after the input gate, respectively. Finally, we can use the updated cells to generate new hidden states through the output gate:

$$h_t=O_t\ast \mathit{tanh}\left(C_t\right)$$

(10)

The role of the output gate is to ensure that only truly valid information is included in the new hidden state. It is worth mentioning that, in order to avoid the disappearance of the gradient information of LSTM through backpropagation, the memory cells are updated by means of addition operations, which allows LSTM to save the cell state for a long time.

2.5. Meta-Learning Algorithm

At present, the monitoring of groundwater level mainly relies on two methods: instrument collection and manual measurement. The collected data are mainly monthly data, and many of them have missing values. Therefore, the amount of data collected for research is generally limited. At present, most deep learning models require large-scale datasets to train the model. When there are few data points, the model usually fails [73]. To address this issue, a meta-learning algorithm was introduced in this study. The meta-learning proposal is aimed at the lack of generalization performance of traditional neural network models and poor adaptability to new types of tasks. It is hoped that the model will acquire a “learning to learn” ability, so that it can acquire existing “knowledge” and quickly learn new tasks.

Considering the structure of our overall model and the training process, the model-agnostic meta-learning (MAML) algorithm [56] was finally chosen. The purpose is to train a deep learning model with better predictions and stronger generalization ability with limited data. We will describe the content of MAML in detail below.

The key idea of MAML is to find a better initial parameter in the process of training the model. Under such conditions, the model can have stronger performance on new tasks after updating the parameters through fewer gradient steps. Suppose there is a distribution $p\left(\mathcal{T}\right)$ for a set of tasks and a model $f_{\theta }$ influenced by the parameter $\theta$. Then, a fraction of tasks ${\mathcal{T}}_i$ is drawn from the task distribution. For each new task ${\mathcal{T}}_i$, its parameters are transformed from the initialized $\theta$ to ${\theta }_i^{\prime }$. In MAML, the updated model parameters ${\theta }_i^{\prime }$ are computed by doing one or more gradient descents on task ${\mathcal{T}}_i$. When considering a gradient calculation:

$${\theta }_i^{\prime }=\theta -\alpha {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)$$

(11)

where $\alpha$ is a hyperparameter, $L_{{\mathcal{T}}_i}$ represents the loss of task ${\mathcal{T}}_i$, and $\nabla$ represents the gradient operation. After a gradient update, the tasks extracted by the algorithm have optimal parameters ${\theta }_i^{\prime }$. Our ultimate goal is for $\theta$ to perform well on any task ${\mathcal{T}}_i$ of the task distribution $p\left(\mathcal{T}\right)$. Therefore, the goal of MAML can be expressed as:

$$\underset{\theta }{\mathit{min}}{\displaystyle \sum }_{{\mathcal{T}}_i~p\left(\mathcal{T}\right)}{L}_{T_i}\left(f_{{\theta }_i^{\prime }}\right)={\displaystyle \sum }_{{\mathcal{T}}_i~p\left(T\right)}{L}_{{\mathcal{T}}_i}\left(f_{\theta -\alpha {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)}\right)$$

(12)

In other words, the goal is to find the best model parameter $\theta$ to minimize the sum of the losses of all tasks after the gradient update. The updating process of the model parameters $\theta$ is as follows:

$$\theta =\theta -\beta {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{{\theta }_i^{\prime }}\right)$$

(13)

where $\beta$ is a hyperparameter. The above steps are also known as meta-optimization or meta-updating.

2.6. CNN-LSTM-ML Prediction Model

2.6.1. Problem Definition

Changes in groundwater levels are usually affected by a variety of factors, including changes in climate, changes in geological and topographical conditions, and anthropogenic influences. The above-mentioned statistical data usually change with the passage of time in different geographical locations. Usually, the change of groundwater level in a region is directly reflected by the water level data of each monitoring point. Due to problems, such as outdated monitoring equipment and insufficient human resources in arid areas, the amount of available data that can be collected is limited. Therefore, in this study, we regard the prediction and modeling process of groundwater levels in arid areas as a multivariate time series prediction problem under the condition of few samples. By predicting the groundwater level of each monitoring well, this method provides a dynamic simulation of groundwater. Here are some symbolic representations of the prediction process:

Definition 1

(Monitoring wells). In this study, the monitoring wells in the data are regarded as $S=\left\{s_1,s_2,\dots ,s_n\right\}$, where $n$ represents the total number of monitoring wells in the study area.

Definition 2

(Relevant features that affect groundwater level). As mentioned above, the meteorological factors that affect the prediction of groundwater level include evaporation, precipitation, pressure, humidity, sunshine time, average temperature and wind speed. Considering a factor influencing each monitoring well $s$ in a time interval $t$ as $X_t$. Then, all the characteristics of a monitoring well $s$ in a time interval $t$ can be ${\mathit{X}}_t=\left\{X_t^1,\dots ,X_t^k\right\}$, where $k$ is the total number of features.

Definition 3

(Groundwater level situation). In this study, the groundwater level of each monitoring well $s$ in a time interval $t$ is represented as $y^t$, which is the prediction target of the model.

Definition 4

(Groundwater level training/prediction task). A groundwater level training/prediction task ${\mathcal{T}}_i$ is composed of a set $D=\left\{{\mathit{X}}_i^t,y_i^t\right\}$ in a continuous time interval $T$. In this study, each task ${\mathcal{T}}_i$ will be divided into training data $D_i^{train}$ and test data $D_i^{test}$, namely ${\mathcal{T}}_i=\left\{\left({\mathit{X}}_i^1,y_i^1\right),\dots ,\left({\mathit{X}}_i^T,y_i^T\right)\right\}=\left\{D_i^{train},D_i^{test}\right\}$.

Usually, to predict the groundwater level of the specified target well $s$, all historical well groundwater level data and meteorological factor data will be put into a model $f_{\varphi }$ affected by the parameter $\varphi$ for training. Then, the trained model is used to make predictions. However, since the goal of this paper is to use a meta-learning algorithm to solve this problem, taking the prediction of a time step as an example, the problem is transformed as follows. We divide the historical groundwater level data and related meteorological feature data into $\tau$ tasks, assuming that all tasks are sampled from the same distribution, i.e., ${\mathcal{T}}_i~p\left(\mathcal{T}\right)$. Then the whole meta-learning process is divided into a meta-training stage and a fine-tune testing stage. In the meta-training stage, our goal is to train so that, under the condition of insufficient data, there is still a model $f_{\theta }$ with high predictive accuracy. In the prediction fine-tune testing stage, our goal can be regarded as putting the trained model $f_{\theta }$ into the target data of few samples $D^{train}=\left\{\left({\mathit{X}}^1,y^1\right),\dots ,\left({\mathit{X}}^{t-1},y^{t-1}\right)\right\}$. The result is a fine-tuned model $f_{{\theta }^{\ast }}$ to complete the prediction task:

$${\widehat{y}}^t=f_{{\theta }^{\ast }}\left(D^{test}\right)$$

(14)

where ${\widehat{y}}^t$ represents the final prediction result

2.6.2. Model Structure

This subsection introduces the network structure designed by CNN-LSTM-ML and the role of each module in capturing the complex groundwater level situation. First, the algorithm samples a batch of tasks ${\mathcal{T}}_i$ from all groundwater level training/prediction tasks $p\left(\mathcal{T}\right)$. Then, the extracted task ${\mathcal{T}}_i$ is divided into training data $D_i^{train}$ and test data $D_i^{test}$. The main task of the training data is to find the optimal parameters for each task through the network structure, while the test data is used to find the optimal parameters of the entire model. Both the training data and the test data go through the same network structure. The input data of the training data will first go through the CNN. CNN is used to capture the multiple influences of each meteorological element on the groundwater level, and, at the same time, it can integrate the spatial relationship between the data, which is convenient for extracting features. The complete sequence is fed into the LSTM layer. Since the observation data is essentially time series data, the feature maps output by the CNN are then fed into the LSTM to learn the long-term dependencies of the sequence. After that, the input of LSTM is put into the fully connected layer to output the prediction result; gradient descent is used to calculate and minimize the loss. Finally, the optimal parameters for the task are obtained. For the input data of the test data, the input network structure is unchanged, but the parameters used by the network are the optimal parameters for the task. On this basis, the gradient descent is performed again to minimize the loss. Finally, we get the relative optimal parameters for the entire model. Figure 6 shows the structure of the training phase of CNN-LSTM-ML.

It is worth noting that the process in the fine-tune testing stage is roughly the same as that in the training stage. The main difference is that the fine-tune testing stage does not need to initialize the model parameters, but rather uses the trained model parameters. Second, the fine-tune testing stage does not have a second gradient update, but instead uses the results of one gradient update to update the model parameters. Subsequent subsections will detail the training and testing algorithm flow.

2.6.3. CNN-LSTM-ML Meta-Training and Fine-Tune Testing

The definition of groundwater level prediction in the field of machine learning is basically a regression problem. We hope to use meta-learning algorithms to solve the problem of small samples in the regression. The goal is to train a statistical function of the groundwater table depth, and ultimately make predictions with fewer data points. Since it is a regression problem, we will use the mean square error loss function, which is described as follows:

$$L_{{\mathcal{T}}_i}\left(f_{\theta }\right)={\displaystyle \sum }_{{\mathit{X}}^j,y^j}\Vert f_{\theta }\left({\mathit{X}}^j\right)-y^j{\Vert }_2^2$$

(15)

where ${\mathit{X}}^j$ and $y^j$ represent input samples and labels sampled from task ${\mathcal{T}}_i$, respectively. ${\Vert ·\Vert }_2$ represents the second norm.

The CNN-LSTM-ML model designed in this study is divided into two processes: meta-training and fine-tune testing (See Algorithms 1 and 2 for details). The purpose of the meta-training stage is to find the optimal initial parameters, so that the model can still make accurate predictions in the case of a small amount of data.

Algorithm 1 CNN-LSTM-ML Meta-Training

Require: Distribution of tasks $p\left(\mathcal{T}\right)$, Step size hyperparameters $\alpha$,$\beta$ 1: Random Initialization parameters $\theta$ 2: While not done do 3: Sample a batch of tasks from $p\left(\mathcal{T}\right)$ i.e.,$\left({\mathcal{T}}_1,{\mathcal{T}}_2,···,{\mathcal{T}}_i\right)\sim p\left(\mathcal{T}\right)$ 4: for all ${\mathcal{T}}_i$ do 5:  Set up training set $D_i^{train}$ for each task in ${\mathcal{T}}_i$ 6:  Calculate ${\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)$ using $D_i^{train}$ and $L_{{\mathcal{T}}_i}$ in Equation (15) 7:  Calculate adapted parameters with gradient descent: ${\theta }_i^{’}=\theta -\alpha {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)$ 8:  Set up test set $D_i^{test}$ for each task in ${\mathcal{T}}_i$ for the meta-update 9: end for 10:  Update $\theta \leftarrow \theta -\beta {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{{\theta }_i^{’}}\right)$ using each $D_i^{test}$ and $L_{{\mathcal{T}}_i}$ in Equation (15) 11: end while

The purpose of the fine-tune testing stage is to optimize the model parameters trained in the meta-training stage on new tasks, and finally make predictions on the test set. The specific process is shown in Algorithm 2:

Algorithm 2 CNN-LSTM-ML Fine-tune Testing

Require: Distribution of tasks $p\left(\mathcal{T}\right)$, Step size hyperparameters $\gamma$ 1: Well-trained parameters $\theta$ 2: While not done do 3: Sample a task ${\mathcal{T}}_i$ from $p\left(\mathcal{T}\right)$ i.e., ${\mathcal{T}}_i\sim p\left(\mathcal{T}\right)$ 4: Set up training set $D_i^{train}$ for each task in ${\mathcal{T}}_i$ 5: Calculate ${\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)$ using $D_i^{train}$ and $L_{{\mathcal{T}}_i}$ in Equation (15) 6: Calculate adapted parameters with gradient descent: ${\theta }_i^{*}=\theta -\gamma {\nabla }_{\theta }{L}_{{\mathcal{T}}_i}\left(f_{\theta }\right)$ 7: Set up test set $D_i^{test}$ for each task in ${\mathcal{T}}_i$ for the prediction 8: Calculate prediction values $\widehat{y}=f_{{\theta }_i^{*}}\left(D_i^{test}\right)$ 9: end while

2.6.4. Model Evaluation

In order to correctly describe the predictive performance of CNN-LSTM-ML, this study uses two commonly used model accuracy evaluation indicators, namely the root mean square error (RMSE) and the mean absolute error (MAE) [74]. The calculation is as follows:

$$\mathrm{MAE}=\frac{1}{N}{\displaystyle \sum }_{n=1}^N\left|\widehat{y_n}-y_n\right|$$

(16)

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\left(\widehat{y_n}-y_n\right)}^2}$$

(17)

where $\widehat{y_n}$ represents the value of the groundwater level predicted by the model, and $y_n$ represents the observed value of the actual groundwater level. $N$ is the data size that needs to be predicted.

3. Results

3.1. Meteorological Factor Time Series

After data processing, the meteorological data of a station in Gaotai County are shown as an example in Figure 7. It can be seen that the annual average temperature is 8.5 °C, the monthly average maximum temperature is 23.5 °C (July), the monthly average minimum temperature is −10.2 °C (January), and the maximum temperature has a rising trend. The average monthly evapotranspiration is greater in spring and summer each year, and the maximum monthly average evaporation is 5 mm. The precipitation in the arid study area is sparse. There is no precipitation for one to two months each year, and the annual average rainfall is about 24 mm. In ten years, the maximum sunshine duration in the study area reached 11.7 h, and the monthly average sunshine duration also reached 8.7 h. The wind speed is relatively stable every year, and the windy weather will occur from time to time in the spring. In addition, the number of months with relative humidity above 60% is relatively small, and most of them are concentrated between 40% and 60%.

In general, the region has a dry climate with less rainfall, high evaporation, large temperature differences between day and night, and long sunshine hours. Therefore, groundwater resources are of great significance to the lives of people in arid areas.

3.2. Prediction Performance

In this section, we use different models to predict the groundwater level of all monitoring wells in the middle and lower reaches of the Heihe River. In the model validation phase, all models including the CNN-LSTM-ML were evaluated by MAE and RMSE. The total evaluation result is the average of the evaluation values of all wells.

Table 1. Predictive accuracy of groundwater level.

	Month	CNN-LSTM-ML	LSTNet	LSTM	MLP	ARIMA	Multiple Regression
	1	0.095	0.123	0.116	0.109	0.154	0.212
	2	0.111	0.154	0.135	0.143	0.176	0.205
MAE	4	0.105	0.141	0.139	0.154	0.209	0.249
	6	0.124	0.152	0.158	0.172	0.223	0.269
	8	0.122	0.168	0.147	0.221	0.238	0.331
	10	0.136	0.187	0.165	0.193	0.244	0.348
	12	0.142	0.201	0.171	0.256	0.251	0.357
	1	0.135	0.172	0.163	0.155	0.21	0.389
	2	0.154	0.202	0.187	0.196	0.248	0.376
RMSE	4	0.148	0.195	0.191	0.237	0.267	0.462
	6	0.172	0.194	0.206	0.266	0.292	0.505
	8	0.175	0.226	0.201	0.301	0.304	0.571
	10	0.183	0.245	0.224	0.298	0.315	0.603
	12	0.196	0.287	0.246	0.429	0.337	0.642

shows the details of the prediction performance.

The prediction results in Table 1 show that the CNN-LSTM-ML model proposed in this study achieves better results in comparison with other models. Especially in the short-term prediction process, the two statistical indicators of CNN-LSTM-ML, MAE and RMSE, are small. As the predicted time step increases, the errors of all models increase, but it is obvious from Table 1 that the error of CNN-LSTM-ML grows slowly, and the predictive performance is the best. This is because the model can perform well in the extraction of multiple meteorological features. It also reflects the adaptability and robustness of the meta-learning algorithm to new tasks. Compared with the classic multiple linear regression method, the MAE and RMSE of the model proposed in this paper are reduced by 0.117 and 0.254, respectively, in the prediction effect of the time step of one month. In the prediction results of twelve months, the MAE and RMSE respectively decrease by 21.5% and 44.6%. Compared with other deep learning models, the MAE of CNN-LSTM-ML is 0.014, 0.021, and 0.028 lower than that of MLP, LSTM, and LSTNet, respectively, in the performance of the prediction step size of one month; and the long-term prediction effect (12 months), RMSE is 23.3%, 5%, and 9.1% better than MLP, LSTM, and LSTNet, respectively. In addition, it can be seen from Table 1 that the long-term prediction ability (12 months) of ARIMA is slightly stronger than that of MLP, but the short-term prediction ability (1 month) of MLP is better than that of other deep learning models, and the prediction performance of CNN-LSTM-ML shows constant stability.

Furthermore, we also evaluate the performance of the CNN-LSTM-ML model by examining the comparative graphical representations. Line charts allow us to better understand the fit between observed and predicted data. The visualization is shown in Figure 8.

Since Table 1 shows the average value of the predictive performance of each model on different monitoring wells, in order to ensure that the mean value is not affected by extreme values, we randomly select monitoring well data for prediction. As can be seen from Figure 8, compared with other models, CNN-LSTM-ML predicts and fits closer to well observations. In addition, during the period from 1 January 2012 to 1 June 2012, the groundwater depth of the observation well shows a significant downward trend. This is predicted by every model. However, it is obvious that the downward trend predicted by CNN-LSTM-ML is closer to the measured value. The multiple linear regression method and ARIMA predict an upward trend in the process, which does not match the actual situation. After 1 June 2012, the depth of groundwater level deepened, and CNN-LSTM-ML does a good job of predicting this change, but the LSTM model appears to predict the time of the change too early, so the error increased. Compared with CNN-LSTM-ML, the prediction errors of other models all increase with the increase of the number of months, and the rate of increase exceeds that of the model proposed in this paper. It can also be seen from Figure 8 that, in terms of the overall prediction and fitting effect, the performance of ARIMA and MLP at some time points is slightly better than that of LSTM, which may be related to the amount of training data. It also shows that the deep learning model using the meta-learning algorithm has strong robustness in few-shot regression.

It is worth mentioning that, in order to verify the predictive ability of the CNN-LSTM-ML model proposed in this paper in the case of insufficient sample size, we performed a reduction operation on the basis of the existing training data volume. By training on datasets reduced by different proportions, the following validation results are obtained, as shown in Figure 9.

Figure 9 shows that, after reducing the amount of training data, the CNN-LSTM-ML model proposed in this paper still achieves better results compared with other models. After reducing the training data by 10%, the MAE and RMSE of CNN-LSTM-ML are reduced by approximately 30% and 17%, respectively, compared to the ARIMA model. When the training data is reduced by 30%, the predictive performance of the CNN-LSTM-ML model is close to that of other models, except for the ARIMA model, but still maintains the best predictive performance. After reducing the training data by 50%, the MAE and RMSE of the CNN-LSTM-ML model are reduced by 0.145 and 0.147, respectively, compared to the multiple linear regression model. The results show that the CNN-LSTM-ML model proposed in this study also has strong predictive ability under the condition of insufficient data. This is related to the meta-learning algorithm we adopted, and also shows that the model can extract other features. It can also be seen from Figure 9 that, after the data volume is reduced by 10%, the predictive ability of the ARIMA model is greatly reduced. It is worth mentioning that under the condition that the training data set is reduced by 10%, the MAE of CNN-LSTM-ML only increases by 0.069 compared with the previous normal amount of training data, and the RMSE increases by 0.21. The MAE and RMSE of LSTM rise by 0.114 and 0.215, respectively. When the training data set is reduced by 50%, CNN-LSTM-ML shows strong robustness, MAE increases by 60%, and RMSE increases by 70% compared with the previous normal amount of training data. However, this figure in LSTM reflects that MAE and RMSE rise by 90.8% and 94.5%, respectively. The performance of the ARIMA is mediocre, with MAE and RMSE rising by 2.612 and 2.538, respectively. However, the MLP and multiple linear regression models show better prediction performance when the amount of training data is reduced.

In addition to the above experiments, this paper also conducts ablation experiments to demonstrate the effectiveness of the meta-learning algorithm that we adopted for the predictive effect of the model. In this experiment, we remove the meta-learning training architecture, use the CNN-LSTM network structure alone for prediction, and compare the results with CNN-LSTM-ML. The experimental results are shown in Figure 10 and Figure 11.

From the results of the ablation experiments, it is clear that the model with the meta-learning algorithm can provide better predictive performance than the model without the meta-learning algorithm, especially under the condition of few samples. From the predictive performance shown in Figure 10, the predictive effect of the CNN-LSTM-ML model is always better than that of the CNN-LSTM in both short-term and long-term prediction. In this experiment, we randomly select a monitoring well, and only use the data of this well for training and prediction. The experimental results are shown in Figure 11. Given the small number of samples available for training, the predicted and observed values of CNN-LSTM-ML are still very close. Compared with CNN-LSTM, the fitting effect of CNN-LSTM is not ideal, which once again confirms that the deep learning model needs a lot of data support. In addition, under the data condition of only one well, it can be seen from the extreme point of 1 September 2021 in Figure 11 that the predictive effect of CNN-LSTM-ML for rapidly changing groundwater levels needs to be strengthened. This is a topic for future study.

4. Discussion

In this study, we propose a method that uses a meta-learning algorithm as the training and prediction framework. CNN extracts feature information and LSTM takes the multivariate feature mixture time series output by CNN as the model input. The model extracts features of relevant groundwater levels and meteorological factors through CNN and calculates long-term dependencies in time series through a LSTM layer. The entire network framework is based on meta-learning, and the model still has strong predictive ability under the condition of a small sample size. In comparative experiments with other baseline models, CNN-LSTM-ML outperforms other models in both short-term and long-term prediction effects. Compared with that of the traditional multiple linear regression model, the prediction performance of CNN-LSTM-ML has been greatly improved. We believe that the result is due to the shortcomings of multiple linear regression in nonlinear fitting, which is also consistent with research [33]. Compared with that of the multivariate time series prediction model LSTNet in deep learning, the predictions of the CNN-LSTM-ML proposed in this paper are outstanding. This optimized result is potentially based on the advantages of the meta-learning algorithm in the regression of limited data [56]. It is possible that LSTNet is more suitable to be applied for time series data with fixed patterns [49], and some of the data involved in this article do not have periodic characteristics. This can also be observed in Table 1, which shows that the predictive performance of LSTM is better than that of LSTNet at some time steps.

Interestingly, the predictive power of the models for observed wells varies widely. This is reflected in the fact that deep learning models are not necessarily better than traditional time series forecasting models. For example, ARIMA performs even slightly better than LSTM at certain time points, and the prediction performance of LSTM suitable for predicting longer time steps may be inferior to MLP on some observation well data. We believe that this may be related to the size of the training data, which can be proved in experiments with reduced datasets. After reducing the training dataset, all model predictions are affected, but CNN-LSTM-ML is more robust. Among them, the impact on ARIMA is tremendous. This consequence may be due to the fact that in the process of reducing the data set, the trend and related seasonal characteristics in the data are also eliminated, so that ARIMA cannot capture this change well, resulting in a declination in forecasting ability. Furthermore, the predictive performance of the deep learning models LSTM and LSTNet is also greatly reduced if the amount of trainable data is greatly reduced. Under the condition that the amount of training data is reduced by 50%, the traditional multiple linear regression model and MLP work better. This shows that the deep learning model is very dependent on the size of the training set. If the available data are in-sufficient, the performance may be similar to that of the traditional model. This is consistent with the findings of [50,51].

The deep learning model is similar to a “black box” structure [75], in which the model has strong predictive ability and weak interpretability. The impact of each meteorological factor on the groundwater level cannot be explained fully from the perspective of the CNN-LSTM-ML. From the perspective of model optimization, attention mechanisms have been widely used in the field of deep learning [76,77,78,79,80,81,82]. An attention mechanism is proposed to simulate the process of human learning by focusing on the most important part of the information. In the case of limited computing power, computing resources can be allocated to more important computing tasks, thereby improving model calculation efficiency. In follow-up research, we will consider adding an attention mechanism to the model. On the one hand, it can improve the calculation speed. On the other hand, the weight assigned by an attention mechanism to each meteorological factor can help explain the changes of groundwater level and improve the model’s interpretability. In addition, according to the requirements of the CNN-LSTM-ML model, the dataset is divided into training data and test data in both the meta-training stage and the fine-tune testing stage. In the process of parameter updating, it undergoes multiple gradient updates. This takes up a lot of computing resources. We will consider optimizing the entire process in future model updates.

In addition to the meteorological factors that affect the groundwater level, the level is also affected by factors, such as landforms, aquifer lithology, and human activities. Subsequent research should consider collecting data on these other factors to facilitate model updating and improve predictive accuracy. In addition, this study uses monthly data, while more fine-grained multivariate time series data could be used in subsequent studies to conduct real-time simulations of groundwater level dynamics.

5. Conclusions

This study proposes a CNN-LSTM-ML model based on a meta-learning algorithm framework, a network structure combining CNN and LSTM, for predicting groundwater levels in the middle and lower reaches of the Heihe River. The structure of CNN+LSTM in this study is designed to solve the problems of the complex nonlinear relationship among multiple meteorological factors and the difficulty of long-time series dependence in the prediction of groundwater level. Because groundwater level monitoring data are limited, this paper uses a meta-learning algorithm to cause the deep learning model to have better initialization parameters that quickly converge, and thereby solves the problem of prediction in the face of limited data. Experimental results show that our proposed CNN-LSTM-ML outperforms other models in both short-term and long-term prediction effects. After reducing the samples available for training, CNN-LSTM-ML still achieves the best prediction performance. Compared with the model with a single CNN-LSTM structure, the CNN-LSTM-ML after adding the meta-learning algorithm improves the accuracy of groundwater level prediction based on the existing sample size, which shows the effectiveness of the meta-learning algorithm. However, the CNN-LSTM-ML model still has room for optimization in terms of interpretability, expansion of multivariate influencing factors, and gradient updates. Future research should start from the above directions. In summary, this paper demonstrates the successful application of the CNN-LSTM-ML model in groundwater level prediction in the middle and lower reaches of the Heihe River. This approach can be extended to other limited-sample time series forecasting with multivariate effects. The model has significance for water resources management in arid regions.