Karst spring flow prediction method based on bayesian optimization algorithm and deep learning and storage medium
By combining CNN and Bayesian optimization, multi-scale hydrological features of karst spring flow are extracted, and parameters are optimized in the LSTM framework to construct a BI-LSTM coupled model. This solves the problems of insufficient feature capture and inadequate parameter optimization in existing karst spring flow prediction models, and achieves high-precision karst spring flow prediction.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- KUNMING PROSPECTING DESIGN INSTITUTE OF CHINA NONFERROUS METALS INDUSTRY CO LTD
- Filing Date
- 2026-04-14
- Publication Date
- 2026-06-05
Smart Images

Figure CN122153404A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of hydrological forecasting technology, specifically relating to a method for predicting karst spring flow based on Bayesian optimization algorithm and deep learning, as well as a storage medium. Background Technology
[0002] Karst areas are a special type of hydrogeological unit, characterized by the coexistence of fissures and conduits in their aquifers. This results in highly uneven spatial distribution of groundwater, which fluctuates significantly seasonally due to factors such as climate and topography. Karst springs, as the primary discharge mechanism for karst groundwater, directly reflect the response characteristics and dynamic patterns of groundwater in karst areas. Accurate and reliable prediction of karst spring flow is a core prerequisite for the rational development, efficient utilization, and ecological environmental protection of water resources in karst areas, and has significant practical implications for regional water resource management, water conservancy project planning, and ecosystem protection.
[0003] Currently, karst spring flow prediction models are mainly divided into two categories: physical models and data-driven models. Physical models describe the dynamic processes of karst groundwater transport and storage by constructing mathematical equations, and rely on analytical or numerical methods to solve these equations to simulate and predict spring flow. However, the hydraulic characteristics of karst aquifers exhibit strong spatial heterogeneity, and factors such as climate fluctuations and human activities cause significant non-stationarity in spring flow. These characteristics are difficult to fully quantify by physical models, resulting in high construction difficulty, limited applicability, and the need for costly field investigations and parameter determination.
[0004] Compared to physical models, data-driven models do not require in-depth characterization of hydrophysical mechanisms; they can achieve prediction simply by exploring the mapping relationship between input and output data. This gives them advantages such as low modeling cost and high computational efficiency, making them the mainstream research direction for karst spring flow prediction. Time series models are the most widely used data-driven models, including Autoregressive Moving Average (ARMA), Recurrent Neural Network (RNN), and Vector Autoregression (VAR). However, existing time series model-based data-driven prediction methods for karst spring flow still have many technical shortcomings, making it difficult to meet the high-precision prediction requirements of practical engineering. Specific problems include the following: 1. Existing models do not adequately capture the input sequence dependency between karst spring flow and meteorological driving factors. The model parameters are mostly determined by empirical settings or simple optimization methods, lacking targeted parameter optimization strategies for non-stationary and highly fluctuating sequences in karst areas. This results in low simulation accuracy of the peak, trough and trend changes of spring flow.
[0005] 2. Since karst areas are mostly located in remote mountainous areas, meteorological stations are sparsely distributed, and meteorological data such as rainfall are scarce. Furthermore, there is often a mismatch between meteorological data and spring flow data in terms of observation time scale, resulting in insufficient effective training samples for the model, which limits the construction and performance improvement of the prediction model.
[0006] 3. Most models simply establish a direct mapping relationship between meteorological factors and spring flow, failing to fully consider the unique storage and transmission characteristics of groundwater in karst areas. They lack sequence analysis and prediction design based on the hydrogeological characteristics of karst areas, resulting in the models being unable to accurately capture the hydrological response patterns of karst spring flow, especially with large prediction errors under extreme rainfall events.
[0007] 4. Existing studies mostly use single methods (such as single variable screening and simple error correction) to improve model performance. They have not yet formed an integrated optimization system that combines spatial feature extraction, global parameter optimization, and time series deep learning. It is difficult to fully explore the spatiotemporal characteristics of karst spring flow series, and the generalization ability and robustness of the model are poor.
[0008] Due to the practical need for accurate prediction of spring flow in water resource management and ecological protection in karst areas, there is a prominent contradiction with the technical shortcomings of existing methods, such as insufficient feature capture, lack of targeted parameter optimization, poor data matching, and failure to adapt to the characteristics of karst hydrology. Therefore, there is an urgent need for an efficient and accurate prediction method that adapts to the hydrological characteristics of karst areas to fill the technical gap. Summary of the Invention
[0009] To address the technical problems of existing methods, this invention provides a method for predicting karst spring flow based on the extraction of spatial hydrological features using a convolutional neural network (CNN) and combined with Bayesian optimization (Bayes) and long short-term memory (LSTM) networks. It also provides a computer storage medium.
[0010] The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning in this invention is implemented as follows: it includes data acquisition and preprocessing, CNN spatial feature extraction, Bayesian parameter optimization, and karst spring flow prediction steps. The specific contents of each step are as follows: A. Data Acquisition and Preprocessing: A right-angled triangular weir was set up below the target spring in the study area and a water level gauge was placed. The water level in front of the weir was continuously observed over a long period of time. The observed water level in front of the weir was converted into the spring flow rate using the weir flow formula. Rainfall data around the spring was collected, and the spring flow rate data and rainfall data were integrated to form the original dataset. The original dataset was preprocessed to remove outliers, unify the scale, and suppress noise to obtain a standardized dataset. B. CNN Spatial Feature Extraction: The standardized dataset is input into the CNN model. The output features of the convolutional layer are obtained through the convolution operation of the CNN model. Then, multiple sets of convolutional kernels are stacked to obtain comprehensive features. Then, average pooling is performed on the comprehensive features to achieve feature dimensionality reduction and noise suppression. Finally, multi-scale hydrological spatial features in the rainfall-spring flow process are extracted. C. Bayesian parameter optimization: Based on multi-scale hydrological spatial characteristics, a surrogate model is established using a Gaussian process, and the EI function is selected as the acquisition function. A Bayesian optimization framework is constructed through the surrogate model and the acquisition function. The optimal value of the log-likelihood function corresponding to the parameter is solved by the gradient descent method using the Bayesian optimization framework. Through multiple iterations, the global optimization of all parameters is achieved to obtain the optimal parameter combination. D. Karst spring flow prediction: The multi-scale hydrological spatial features extracted in step B are imported into the LSTM framework to construct a BI-LSTM coupled prediction model. Then, the optimal parameter combination obtained in step C is substituted into the BI-LSTM coupled prediction model. Through the synergistic effect of the model's forget gate, input gate, cell state update and output gate, the karst spring flow is predicted.
[0011] Furthermore, in step A, the original dataset is preprocessed using MATLAB packages, including the following sub-steps: A10. Outlier Removal: The 3σ principle outlier detection method is used to identify and remove extreme outlier data points in the original dataset, retaining valid data. A20. Data Normalization: Perform IQR normalization on the dataset after removing extreme outliers; A30. Noise Suppression: For the dataset after IQR normalization, calculate the average value of all data within a preset window before and after each data point, and replace the original value of the corresponding data point with this average value to obtain the standardized dataset.
[0012] Furthermore, in the IQR normalization process, the IQR calculation formula is as follows: IQR = Q 3- Q 1, In the formula: Q 3 represents the data value at the 75th percentile after the original dataset is sorted in ascending order; Q 1 represents the data value at the 25th percentile after the original dataset is sorted in ascending order.
[0013] Furthermore, in step B, the specific mathematical expression and process of the convolution operation, stacking of multiple convolutional kernels, and average pooling operation of the CNN model are as follows: B10. The process data sequence X input to the CNN model is represented as: X = { x1 , x 2 ,..., x T}, x T ∈R d , In the formula: T For time step, d Let R be the input feature dimension, and R be the real number field. B20. The comprehensive feature H obtained after stacking multiple sets of convolutional kernels from the output features of the convolutional layer is represented as follows: , In the formula: Concat It is K Each column vector (or row vector) is placed side-by-side by column (or by channel) to form a matrix; K The number of convolution kernels, H The final output of the convolution kernel; B30. The mathematical form of the average pooling operation is: , In the formula: This represents the output value of the j-th feature channel after average pooling. h j The pooling window size corresponding to the j-th feature channel. pool This is an average pooling operation; r This represents the size of the pooling window, and also the index range of positions within the pooling window; j This is the feature channel index, and also indicates that the current process is the first [number]th [channel] in that channel. j A pooled window, h jr Let be the feature value at position r within the j-th feature channel and the j-th pooling window.
[0014] Furthermore, in step C, the EI function, which serves as the acquisition function in the Bayesian optimization algorithm, has the following mathematical expression: , In the formula: α EI ( x ) is the EI function with respect to the independent variable x The function value at that point is the expected improvement value; Φ is the cumulative distribution function of the standard normal distribution. f Let be the probability density function of the standard normal distribution. f ( x + This represents the best predicted value obtained by the surrogate model so far. m ( x) is the proxy model in the independent variable x The predicted mean at that location, s ( x ) is the proxy model in the independent variable x The predicted standard deviation at point Z; Z is a standardized variable; in, ; The selection logic for the EI function is: when the independent variable... x Predicted mean at m ( x (Better than the current best prediction) f ( x + When this happens, the EI function tends to select that independent variable. x The optimal predicted value is the best objective function value that the surrogate model believes to be so far. Wherein, the optimal objective function value f ( x + It is usually defined as: , In the formula: x i For the first in the training set i One sample point; X train The training dataset contains all evaluated datasets. x i ; m ( x i )for x i The predicted mean at the location; Or when the independent variable x Predicted standard deviation at s ( x When the value is large, the EI function will assign a value to the independent variable. x A certain selection weight.
[0015] Furthermore, in step C, the specific process of using a Bayesian optimization framework to solve for the optimal value of the log-likelihood function corresponding to the parameters through gradient descent is as follows: C10: The probability density functions corresponding to each parameter are: , In the formula: x 1 represents the set of the first t observed data points, i.e. x 1 , x 2 , ..., x t ; t is the number of observed data points, η is the noise term, y tLet be the observed value corresponding to the t-th observation data point. m t For the surrogate model at the t-th observation data point x t The predicted mean at point T is the time step; C20: Due to the set of observation data points x 1:t and corresponding observations f(x 1:t ) Since all parameters are known, the parameter values that maximize the aforementioned probability density function are transformed into solving for the maximum value of the log-likelihood function, which is: , C30: Maximizing the above log-likelihood function is equivalent to minimizing its negative function, i.e., minimizing: , C40: Following the mechanism and process of solving the optimal value of probability density as described above, Bayesian optimization is performed on all other parameters in the BI-LSTM coupled prediction model and the surrogate model one by one, and finally the optimal values of all parameters are obtained, forming the optimal parameter combination.
[0016] Furthermore, in step D, the specific structure, mathematical expression, and function of the forget gate, input gate, cell state update, and output gate of the BI-LSTM coupled prediction model are as follows: Forget gate: Used to selectively filter historical information, deciding whether to discard or retain information from a previous moment; it uses an activation function. s Process the hidden state from the previous moment h t-1 and current input x t Generate the gate value; The mathematical expression is: , In the formula: f t The activation value for the forget gate, with a range of [0,1]; s The Sigmoid activation function has an output range of [0,1]. W f This is the weight matrix; h t-1 The hidden state from the previous moment; x t The input vector at the current time step; b f For the bias term of the forget gate; Input gate: Used to determine whether the current input information is updated to the cell state, providing new information to be stored in the memory unit; includes input gate activation value. i t It is a gated unit based on the activation function σ, which determines h t-1 and x t Whether the information in the memory is updated; and selectively updating the cell state to provide new information to be stored in the memory unit; The mathematical expression for candidate memory cell states is: , The mathematical expression for the input gate activation value is: , In the formula: represents the candidate memory cell state; tanh is the hyperbolic tangent activation function with an output range of [-1, 1]; W C , w i These are the weight matrices for the input gate and the candidate memory cell states, respectively. h t-1 The hidden state from the previous moment; x t This is the input vector at the current moment; i t The input gate activation value, which takes the range [0,1]. bi , b C These are the bias terms for the input gate and the candidate memory cell state, respectively; Cell state update: The cell state is updated by updating the cell state at the previous time step. C t-1 Output with the Forgot Gate f t Multiplication enables selective forgetting; then the input gate output is superimposed. i t With candidate state The product of these factors generates the updated cell state. The mathematical expression is: , In the formula: C t This represents the current state of the cell. Output gate: Used to control the output of cell state, and the final model output information is determined through the gating mechanism; The mathematical expression for the output gate activation value is: , The mathematical expression for the hidden state at the current moment is: , In the formula: ht This represents the hidden state at the current moment, i.e., the output value of the BI-LSTM coupled prediction model; O t The activation value of the output gate, with a value range of [0,1]; W o This is the weight matrix of the output gate; b o This is the bias term for the output gate.
[0017] Furthermore, the present invention also includes a model performance evaluation step, which is performed after step D, using RMSE, MAE, and R... 2 The predictive performance of the BI-LSTM coupled prediction model is comprehensively evaluated using four evaluation metrics, including MAPE, residual analysis plots, Taylor plots, and time series comparison plots.
[0018] Furthermore, the RMSE is the root mean square error, which is used to measure the sample standard deviation of the difference between the predicted value and the observed value; The formula for calculating RMSE is: , The MAE is the Mean Absolute Error, which measures the average absolute error between the predicted and actual values. The formula for calculating MAE is: , The R 2 The coefficient of determination is used to measure how well the model fits the data. R 2 The calculation formula is: , The MAPE stands for Mean Absolute Percentage Error, which measures the average percentage difference between the predicted value and the actual observed value. The formula for calculating MAPE is: , In the formula: m For the number of samples, y i Let i be the true value of the i-th sample. Let be the predicted value for the i-th sample. This is the average of the true values of all samples.
[0019] The computer-readable storage medium of the present invention is implemented as follows: a computer program is stored thereon, which can be executed by one or more processors to implement the karst spring flow prediction method based on Bayesian optimization algorithm and deep learning as described above.
[0020] The present invention has the following beneficial effects: 1. The front end of this invention leverages the powerful spatial feature extraction capabilities of CNNs to extract hydrological spatial features. Through convolution operations and stacking multiple sets of convolutional kernels, it deeply mines multi-scale spatial features in the rainfall-spring flow process. Then, Bayesian optimization is used to automatically iteratively optimize the model parameters. The back end combines the strong ability of LSTM to capture the dependencies of long-term time-series data, thereby constructing a spring flow BI-LSTM coupled prediction model. This model can adapt to the complex nonlinear characteristics of non-stationary and strongly fluctuating karst spring flow, thus significantly improving the response capability and prediction accuracy to the peak, valley and trend changes of spring flow under extreme hydrological events, making up for the shortcomings of insufficient feature capture in traditional models.
[0021] 2. This invention introduces a Bayesian optimization algorithm, constructs a surrogate model using a Gaussian process, employs an EI function as the data acquisition function, and integrates it with a deep learning framework to build an intelligent parameter optimization mechanism. This mechanism enables rapid global optimization of the model's high-dimensional parameters, effectively avoiding getting trapped in local optima. Furthermore, this optimization mechanism tailors parameter strategies to the hydrological characteristics of karst areas, significantly improving the model's generalization ability and robustness under different karst geological backgrounds and hydrological conditions. This solves the problem of insufficient stability and accuracy caused by the lack of specificity in parameter settings in existing models.
[0022] 3. This invention constructs a purely data-driven prediction architecture that does not rely on complex prior geological information such as aquifer structure, fracture network distribution, and groundwater flow field. It can model and perform high-precision predictions solely using long-term monitored spring flow and rainfall time-series data, avoiding the high costs of field exploration in karst areas. Especially in karst areas where the hydrogeological structure is unclear, the evolution of spring flow can be accurately predicted by analyzing the correlation between historical spring flow and rainfall. Furthermore, data preprocessing techniques such as 3σ outlier removal, IQR normalization, and sliding window noise suppression effectively improve sample quality issues caused by scarce meteorological data and spatiotemporal scale mismatches, significantly reducing the model's dependence on data volume and greatly lowering the application threshold and cost.
[0023] 4. This invention constructs an integrated system of "spatial feature extraction - global parameter optimization - temporal deep prediction," realizing full-process optimization from feature mining and parameter calibration to model prediction; combined with RMSE, MAE, and R... 2 The model uses four quantitative indicators (MAPE) and visualization tools such as residual analysis charts and Taylor charts to comprehensively evaluate model performance, ensuring the accuracy and reliability of prediction results and providing scientific data support for water resource management and ecological protection in karst areas.
[0024] In summary, this invention utilizes CNN to extract multi-scale spatial features, Bayesian optimization algorithm to achieve global parameter optimization, and BI-LSTM to capture temporal nonlinear relationships, constructing a purely data-driven integrated prediction system. This system can accurately adapt to the hydrological characteristics of karst areas, improve the prediction accuracy of extreme events and the generalization robustness of the model, and reduce the application threshold and cost, providing reliable technical support for water resource management in karst areas. Attached Figure Description
[0025] Figure 1 This is a flowchart of the karst spring flow prediction method of the present invention; Figure 2 This is a framework diagram of the coupling model in the karst spring flow prediction method of the present invention; Figure 3 This is a fitting effect diagram of the BI-LSTM coupled prediction model in an embodiment of the present invention; In the diagram: 3a is the training set, and 3b is the test set; Figure 4 This is a radar chart showing the model performance of an embodiment of the present invention; Figure 5 This is a residual analysis diagram of an embodiment of the present invention; Figure 6 This is a Taylor diagram of an embodiment of the present invention; Figure 7 This is a comparison chart of the various models and measured water volumes in the embodiments of the present invention. Detailed Implementation
[0026] The present invention will be further described below with reference to the accompanying drawings and embodiments, but this does not limit the present invention in any way. Any changes or improvements made based on the teachings of the present invention shall fall within the protection scope of the present invention.
[0027] like Figure 1 and 2 As shown, this invention presents a karst spring flow prediction method based on Bayesian optimization algorithm and deep learning, including data acquisition and preprocessing, CNN spatial feature extraction, Bayesian parameter optimization, and karst spring flow prediction steps. The specific details of each step are as follows: A. Data Acquisition and Preprocessing: A right-angled triangular weir was set up below the target spring in the study area, and a water level gauge was placed. The water level in front of the weir was continuously observed over a long period of time. The observed water level in front of the weir was converted into spring flow rate using the weir flow formula. Rainfall data from the surrounding area of the spring (i.e., the meteorological station closest to the spring) was collected. The spring flow rate data and rainfall data were integrated to form a raw dataset including spring flow rate, rainfall, and time data. The raw dataset was preprocessed to remove outliers, unify the scale, and suppress noise to obtain a standardized dataset. B. CNN Spatial Feature Extraction: The standardized dataset is input into the CNN (Convolutional Neural Network) model. The convolutional operation of the CNN model is used to obtain the output features of the convolutional layer. Then, multiple sets of convolutional kernels are stacked to obtain the comprehensive features. Then, the average pooling operation is performed on the comprehensive features to achieve feature dimensionality reduction and noise suppression. Finally, the multi-scale hydrological spatial features in the rainfall-spring flow process are extracted. C. Bayesian parameter optimization: Based on multi-scale hydrological spatial characteristics, a surrogate model is established using Gaussian processes (Mockuse et al., 1978). The EI (Expected Improvement) function (Jones, 2001) is selected as the acquisition function. Bayesian optimization (Bayesian Optimization) is constructed through the surrogate model and the acquisition function. For Bayesian formulas, please refer to J Mockus, V Tiesis, and A Zilinskas. The application of Bayesian methods for seeking the extremum. Towards Global Optimization, 2:117–129, 1978. and DR Jones. A taxonomy of global optimization methods based on response surfaces. Journal of Global Optimization, 21(4):345–383. (2001.) Framework; The Bayesian optimization framework is used to solve for the optimal value of the log-likelihood function corresponding to the parameters (such as learning rate, regularization coefficient, etc.) through gradient descent. Through multiple rounds of iteration, global optimization of all parameters is achieved to obtain the optimal parameter combination. D. Karst spring flow prediction: The multi-scale hydrological spatial features extracted in step B are imported into the LSTM (Long Short-Term Memory) framework to construct a BI-LSTM coupled prediction model. Then, the optimal parameter combination obtained in step C is substituted into the BI-LSTM coupled prediction model. Through the synergistic effect of the model's forget gate, input gate, cell state update and output gate, the karst spring flow is predicted.
[0028] It should be noted that the synergistic effect in step A is that after constructing the BI-LSTM coupled model, the model will continuously optimize intelligently during operation to obtain the optimal solution for each parameter; then, using the spring flow dataset as input, with the parameters set as the optimal solution, the spring flow is predicted through the gradient descent algorithm combined with the gating mechanism of LSTM, and various performance indicators are generated.
[0029] In step A, the original dataset is preprocessed using MATLAB packages, including the following sub-steps: A10. Outlier Removal: The 3σ principle outlier detection method is used to identify and remove extreme outlier data points in the original dataset, retaining valid data. A20. Data Normalization: Perform IQR normalization on the dataset after removing extreme outliers; A30. Noise Suppression: For the dataset after IQR normalization, calculate the average value of all data within a preset window before and after each data point, and replace the original value of the corresponding data point with this average value to obtain the standardized dataset.
[0030] In the IQR normalization process, the IQR calculation formula is as follows: IQR = Q 3- Q 1, In the formula: Q 3 represents the data value at the 75th percentile after the original dataset is sorted in ascending order; Q 1 represents the data value at the 25th percentile after the original dataset is sorted in ascending order.
[0031] In step B, the specific mathematical expression and process of the convolution operation, stacking of multiple convolutional kernels, and average pooling operation of the CNN model are as follows: B10. The process data sequence X input to the CNN model is represented as: X = { x 1 , x 2 ,..., x T}, x T ∈R d , In the formula: T For time step, d Let R be the input feature dimension, and R be the real number field. B20. The comprehensive feature H obtained after stacking multiple sets of convolutional kernels from the output features of the convolutional layer is represented as follows: , In the formula: Concat It is K Each column vector (or row vector) is placed side-by-side by column (or by channel) to form a matrix; K The number of convolution kernels, H The final output of the convolution kernel; B30. The mathematical form of the average pooling operation is: , In the formula: This represents the output value of the j-th feature channel after average pooling. h j The pooling window size corresponding to the j-th feature channel. pool This is an average pooling operation; r This represents the size of the pooling window, and also the index range of positions within the pooling window; j This is the feature channel index, and also indicates that the current process is the first [number]th [channel] in that channel. j A pooled window, h jr Let be the feature value at position r within the j-th feature channel and the j-th pooling window.
[0032] In step C, the EI function, which serves as the acquisition function in the Bayesian optimization algorithm, has the following mathematical expression: , In the formula: α EI ( x ) is the EI function with respect to the independent variable x The function value at that point is the expected improvement value; Φ is the cumulative distribution function of the standard normal distribution. f Let be the probability density function of the standard normal distribution. f ( x + This represents the best predicted value obtained by the surrogate model so far. m ( x ) is the proxy model in the independent variable x The predicted mean at that location, s ( x ) is the proxy model in the independent variable x The predicted standard deviation at point Z; Z is a standardized variable; in, ; The selection logic for the EI function is: when the independent variable... x Predicted mean at m ( x (Better than the current best prediction) f ( x + When this happens, the EI function tends to select that independent variable. x The optimal predicted value is the best objective function value that the surrogate model believes to be so far. Wherein, the optimal objective function value f ( x + It is usually defined as: , In the formula: xi For the first in the training set i One sample point; X train The training dataset contains all evaluated datasets. x i ; m ( x i )for x i The predicted mean at the location; Or when the independent variable x Predicted standard deviation at s ( x When the value is large, the EI function will assign a value to the independent variable. x Certain selection weights (the weights are determined adaptively by the data).
[0033] It should be noted that in machine learning, the predicted mean, variance, and standard deviation can be calculated directly through the surrogate model.
[0034] In step C, the specific process of using a Bayesian optimization framework to solve for the optimal value of the log-likelihood function corresponding to the parameters through gradient descent is as follows: C10: The probability density functions corresponding to each parameter are: , In the formula: x 1 represents the set of the first t observed data points, i.e. x 1 , x 2 , ..., x t ; t is the number of observed data points, η is the noise term, y t Let be the observed value corresponding to the t-th observation data point. m t For the surrogate model at the t-th observation data point x t The predicted mean at point T is the time step; C20: Due to the set of observation data points x 1:t and corresponding observations f(x 1:t ) Since all parameters are known, the parameter values that maximize the aforementioned probability density function are transformed into solving for the maximum value of the log-likelihood function, which is: , C30: Maximizing the above log-likelihood function is equivalent to minimizing its negative function, i.e., minimizing: , C40: Following the mechanism and process of solving the optimal value of probability density as described above, Bayesian optimization is performed on all other parameters in the BI-LSTM coupled prediction model and the surrogate model one by one, and finally the optimal values of all parameters are obtained, forming the optimal parameter combination.
[0035] In step D, the specific structure, mathematical expression, and function of the forget gate, input gate, cell state update, and output gate of the BI-LSTM coupled prediction model are as follows: Forget gate: Used to selectively filter historical information, deciding whether to discard or retain information from a previous moment; it uses an activation function. s Process the hidden state from the previous moment h t-1 and current input x t Generate the gate value; The mathematical expression is: , In the formula: f t σ is the activation value of the forget gate, with a value range of [0,1]; σ is the Sigmoid activation function, with an output range of [0,1]. W f This is the weight matrix; h t-1 The hidden state from the previous moment; x t The input vector at the current time step; b f For the bias term of the forget gate; Input gate: Used to determine whether the current input information is updated to the cell state, providing new information to be stored in the memory unit; includes input gate activation value. i t It is a gated unit based on the activation function σ, which determines h t-1 and x t Whether the information in the memory is updated; and selectively updating the cell state to provide new information to be stored in the memory unit; The mathematical expression for candidate memory cell states is: , The mathematical expression for the input gate activation value is: , In the formula: represents the candidate memory cell state; tanh is the hyperbolic tangent activation function with an output range of [-1, 1]; W C , w i These are the weight matrices for the input gate and the candidate memory cell states, respectively.h t-1 The hidden state from the previous moment; x t This is the input vector at the current moment; i t The input gate activation value, which takes the range [0,1]. bi , b C These are the bias terms for the input gate and the candidate memory cell state, respectively; Cell state update: The cell state is updated by updating the cell state at the previous time step. C t-1 Output with the Forgot Gate f t Multiplication enables selective forgetting; then the input gate output is superimposed. i t With candidate state The product of these factors generates the updated cell state. The mathematical expression is: , In the formula: C t This represents the current state of the cell. Output gate: Used to control the output of cell state, and determines the final model output information through gating mechanism; The mathematical expression for the output gate activation value is: , The mathematical expression for the hidden state at the current moment is: , In the formula: h t This represents the hidden state at the current moment, i.e., the output value of the BI-LSTM coupled prediction model; O t The activation value of the output gate, with a value range of [0,1]; W o This is the weight matrix of the output gate; b o This is the bias term for the output gate.
[0036] This invention also includes a model performance evaluation step, which is performed after step D, using RMSE, MAE, and R... 2 The predictive performance of the BI-LSTM coupled prediction model is comprehensively evaluated using four evaluation metrics, including MAPE, residual analysis plots, Taylor plots, and time series comparison plots.
[0037] The RMSE is the root mean square error, which measures the sample standard deviation of the difference between the predicted value and the observed value. It reflects the overall level of prediction error. The smaller the value, the higher the prediction accuracy. The formula for calculating RMSE is: , The MAE is the mean absolute error, which measures the average absolute error between the predicted value and the actual value. It reflects the average level of the prediction error. The smaller the value, the better the prediction effect. The formula for calculating MAE is: , The R 2 The coefficient of determination measures how well the model fits the data. Its value ranges from (-∞, 1]. The closer it is to 1, the better the model fits the data and the stronger the consistency between the predicted and actual values. R 2 The calculation formula is: , The MAPE stands for Mean Absolute Percentage Error, which measures the average percentage difference between the predicted value and the actual observed value. It intuitively reflects the relative magnitude of the prediction error; the smaller the value, the higher the prediction accuracy. The formula for calculating MAPE is: , In the formula: m For the number of samples, y i This represents the true value of the i-th sample (i.e., the actual observed flow rate of the karst spring). This is the predicted value for the i-th sample (i.e., the karst spring flow rate output by the BI-LSTM coupled prediction model). This is the average of the true values of all samples.
[0038] The present invention provides a computer-readable storage medium having a computer program stored thereon, the computer program being executable by one or more processors to implement the aforementioned method for predicting karst spring flow based on Bayesian optimization algorithm and deep learning.
[0039] Example
[0040] This study uses the karst spring flow rate in central Guizhou Province as a case study. This region has a mild climate, abundant rainfall, and its main water supply source is precipitation, with a relatively short groundwater flow path. The specific details are as follows: S100: By setting up a right-angled triangular weir and placing a water level gauge below the target spring in the study area, the water level in front of the weir from September 2021 to April 2024 was obtained. The observed water level in front of the weir was converted into spring flow using the weir flow formula. Rainfall data from the meteorological station closest to the spring was collected, and the spring flow data and rainfall data were integrated to form the original dataset.
[0041] Then, by employing the 3σ outlier detection method, extreme outliers in the original dataset are identified and removed, retaining only valid data. Next, the dataset after removing extreme outliers undergoes IQR normalization to standardize the data scale and eliminate the influence of unit dimensions. Subsequently, for the IQR-normalized dataset, the average value of all data within a preset window before and after each data point is calculated. This average value replaces the original value of the corresponding data point, achieving data smoothing and reducing random noise interference. This enhances the feature representation of the data, improves the stability of the time series, and yields a standardized dataset.
[0042] S200: The standardized dataset is input into the CNN model. The convolutional layer output features are obtained through the convolutional operation of the CNN model. Then, multiple sets of convolutional kernels are stacked to obtain comprehensive features. Then, average pooling is performed on the comprehensive features to achieve feature dimensionality reduction and noise suppression. Finally, multi-scale hydrological spatial features in the rainfall-spring flow process are extracted.
[0043] CNN models are used to extract local features and temporal correlations from rainfall data. The input data sequence X to the CNN model is represented as: X = { x 1 , x 2 ,..., x T}, x T ∈R d , In the formula: T The time step represents the time dimension of the data sequence; d is the input feature dimension, representing the number of features in the input data at each time step; R is the real number field; The comprehensive feature H obtained after the output features of the convolutional layer are stacked with multiple sets of convolutional kernels is represented as: , In the formula: Concat It is K Each column vector (or row vector) is placed side-by-side by column (or by channel) to form a matrix; K The number of convolution kernels, H The final output of the convolution kernel; Subsequently, the convolutional features can be reduced in dimensionality and suppressed for noise using average pooling, where the mathematical form of average pooling is: , In the formula: This represents the output value of the j-th feature channel after average pooling. h j The pooling window size corresponding to the j-th feature channel. poolThis is an average pooling operation; r This represents the size of the pooling window, and also the index range of positions within the pooling window; j This is the feature channel index, and also indicates that the current process is the first [number]th [channel] in that channel. j A pooled window, h jr Let be the feature value at position r within the j-th feature channel and the j-th pooling window.
[0044] Through the above convolution and pooling operations, the model can effectively extract multi-scale data features during the rainfall process, providing information-rich input feature representations for subsequent sequence modeling.
[0045] S300: Based on multi-scale hydrological spatial characteristics, a surrogate model is established using Gaussian processes (Mockus et al., 1978). The EI (Expected Improvement) function (Jones, 2001) is selected as the acquisition function. A Bayesian optimization (Bayes) framework is constructed through the surrogate model and the acquisition function. The Bayesian optimization framework is used to solve for the optimal value of the log-likelihood function corresponding to the parameters (such as learning rate, regularization coefficient, etc.) through gradient descent. Through multiple iterations, global optimization of all parameters is achieved to obtain the optimal parameter combination.
[0046] The mathematical expression of the EI function, which serves as the acquisition function in the Bayesian optimization algorithm, is as follows: , In the formula: α EI ( x ) is the EI function with respect to the independent variable x The function value at that point is the expected improvement value; Φ is the cumulative distribution function of the standard normal distribution. f Let be the probability density function of the standard normal distribution. f ( x + This represents the best predicted value obtained by the surrogate model so far. m ( x ) is the proxy model in the independent variable x The predicted mean at that location, s ( x ) is the proxy model in the independent variable x The standard deviation of the prediction at point Z; Z is the standardized variable.
[0047] in, ; The selection logic for the EI function is: when the independent variable... x Predicted mean at m (x (Better than the current best prediction) f ( x + When this happens, the EI function tends to select that independent variable. x The optimal predicted value is the best objective function value that the surrogate model believes to be so far; the optimal objective function value f ( x + It is usually defined as: , In the formula: x i For the first in the training set i One sample point; X train The training dataset contains all evaluated datasets. x i ; m ( x i )for x i The predicted mean at the location; Or when the independent variable x Predicted standard deviation at s ( x When the value is large, the EI function will assign a value to the independent variable. x Certain selection weights (the weights are determined adaptively by the data).
[0048] The specific process of using a Bayesian optimization framework and gradient descent to solve for the optimal value of the log-likelihood function corresponding to the parameters is as follows: S310: The probability density functions corresponding to each parameter are: , In the formula: x 1 represents the set of the first t observed data points, i.e. x 1 , x 2 , ..., x t ; t is the number of observed data points, η is the noise term, y t Let be the observed value corresponding to the t-th observation data point. m t For the surrogate model at the t-th observation data point x t The predicted mean at point T is the time step.
[0049] S320: Due to the set of observation data points x 1:t and corresponding observations f(x1:t ) Since all parameters are known, the parameter values that maximize the aforementioned probability density function are transformed into solving for the maximum value of the log-likelihood function, which is: .
[0050] S330: Maximizing the above log-likelihood function is equivalent to minimizing its negative function, i.e., minimizing: .
[0051] S340: Following the mechanism and process of solving the optimal value of probability density as described above, Bayesian optimization is performed on all other parameters in the BI-LSTM coupled prediction model and surrogate model one by one, and finally the optimal values of all parameters are obtained, forming the optimal parameter combination.
[0052] S400: Import the multi-scale hydrological spatial features extracted from S200 into the LSTM framework to construct a BI-LSTM coupled prediction model (e.g., ...). Figure 2 As shown in Table 1, the model parameters are then substituted into the BI-LSTM coupled prediction model. The flow rate of the karst spring is predicted through the synergistic effect of the forget gate, input gate, cell state update and output gate of the model.
[0053] Table 1 Model Parameter Settings The specific structures, mathematical expressions, and functions of the forgetting gate, input gate, cell state update, and output gate in the BI-LSTM coupled prediction model are as follows: Forget gate: Used to selectively filter historical information, deciding whether to discard or retain information from a previous moment; it uses an activation function. s Process the hidden state from the previous moment h t-1 and current input x t Generate the gate value; The mathematical expression is: , In the formula: f t σ is the activation value of the forget gate, with a value range of [0,1]; σ is the Sigmoid activation function, with an output range of [0,1]. W f This is the weight matrix; h t-1 The hidden state from the previous moment; x t The input vector at the current time step; b f For the bias term of the forget gate; Input gate: Used to determine whether the current input information is updated to the cell state, providing new information to be stored in the memory unit; includes input gate activation value. i t It is a gated unit based on the activation function σ, which determines h t-1 and x t Whether the information in the memory is updated; and selectively updating the cell state to provide new information to be stored in the memory unit; The mathematical expression for candidate memory cell states is: , The mathematical expression for the input gate activation value is: , In the formula: represents the candidate memory cell state; tanh is the hyperbolic tangent activation function with an output range of [-1, 1]; W C , w i These are the weight matrices for the input gate and the candidate memory cell states, respectively. h t-1 The hidden state from the previous moment; x t This is the input vector at the current moment; i t The input gate activation value, which takes the range [0,1]. bi , b C These are the bias terms for the input gate and the candidate memory cell state, respectively; Cell state update: The cell state is updated by updating the cell state at the previous time step. C t-1 Output with the Forgot Gate f t Multiplication enables selective forgetting; then the input gate output is superimposed. i t With candidate state The product of and generates the updated cell state.
[0054] The mathematical expression is: , In the formula: C t This represents the current state of the cell. Output gate: Used to control the output of cell state, and determines the final model output information through gating mechanism; The mathematical expression for the output gate activation value is: , The mathematical expression for the hidden state at the current moment is: , In the formula: h t This represents the hidden state at the current moment, i.e., the output value of the BI-LSTM coupled prediction model; O t The activation value of the output gate, with a value range of [0,1]; W o This is the weight matrix of the output gate; b o This is the bias term for the output gate.
[0055] The final BI-LSTM coupled prediction model uses Adam as the optimizer, with a maximum of 2400 iterations, 8 iterations per epoch, and 400 training epochs. The CNN convolutional kernel size is 3, and the number of kernels is 128. The first LSTM layer has 89 nodes, the second LSTM layer has 70 nodes, the learning rate is 0.001473, the learning rate decay factor is 0.5 (the closer to 1.0, the faster the learning rate decays), and the Dropout rate is 0.210. The L2 regularization coefficient is 0.00001245, and the training set to test set ratio is 7:3. The final performance of the BI-LSTM coupled prediction model is as follows: Figure 3 .
[0056] S500: To more comprehensively evaluate the predictive performance of the BI-LSTM coupled prediction model, RMSE, MAE, and R... 2 The evaluation is conducted using four indicators, including MAPE, combined with residual analysis plots, Taylor plots, and time series comparison plots.
[0057] The aforementioned BI-LSTM coupled prediction model was trained, validated, and tested with CNN, RNN, and LSTM models on the spring flow dataset. The observed flow sequences and predicted flow curves for each model were plotted. The standard deviation and variance of each model compared to the measured values were calculated, and residual analysis plots, Taylor plots, and comparison plots with measured spring flow were generated. The results are shown below. Figure 5 , 6 As shown in Figure 7.
[0058] Comprehensive evaluation results show that the BI-LSTM coupled prediction model's predictions of spring flow closely match the actual trend, with most predicted values closely matching the actual values. Deviations are only observed at inflection points and for some peak values, demonstrating its ability to predict flow changes effectively. The BI-LSTM coupled prediction model more accurately captures the details and local characteristics of actual sequence fluctuations, and its generalization of typical hydrological features such as flood peak response is more accurate. Through optimization of model parameters, cumulative errors are controlled to some extent, and the number of outliers is significantly reduced. Residual analysis plots, Taylor plots, and radar charts demonstrate the predictive performance of each model. Compared with other benchmark models, the BI-LSTM coupled prediction model exhibits the best performance (R0).2 =0.974, RMSE=12.3, MAPE=2.89%, indicating that the hybrid framework proposed in this invention has strong anti-interference ability and good generalization ability in long-term series prediction.
[0059] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.
Claims
1. A method for predicting karst spring flow based on Bayesian optimization algorithm and deep learning, characterized in that: The process includes data acquisition and preprocessing, CNN spatial feature extraction, Bayesian parameter optimization, and karst spring flow prediction. The specific details of each step are as follows: A. Data acquisition and preprocessing: A right-angled triangular weir was set up below the target spring in the study area and a water level gauge was placed. The water level in front of the weir was continuously observed over a long period of time. The observed water level in front of the weir was converted into the spring flow rate using the weir flow formula. Collect rainfall data around the spring, and integrate the spring flow data with the rainfall data to form the original dataset; The original dataset is preprocessed to remove outliers, unify the scale, and suppress noise, resulting in a standardized dataset. B. CNN Spatial Feature Extraction: The standardized dataset is input into the CNN model. The output features of the convolutional layer are obtained through the convolution operation of the CNN model. Then, multiple sets of convolutional kernels are stacked to obtain comprehensive features. Then, average pooling is performed on the comprehensive features to achieve feature dimensionality reduction and noise suppression. Finally, multi-scale hydrological spatial features in the rainfall-spring flow process are extracted. C. Bayesian parameter optimization: Based on multi-scale hydrological spatial characteristics, a surrogate model is established using a Gaussian process, and the EI function is selected as the acquisition function. A Bayesian optimization framework is constructed through the surrogate model and the acquisition function. The optimal value of the log-likelihood function corresponding to the parameter is solved by the gradient descent method using the Bayesian optimization framework. Through multiple iterations, the global optimization of all parameters is achieved to obtain the optimal parameter combination. D. Karst spring flow prediction: The multi-scale hydrological spatial features extracted in step B are imported into the LSTM framework to construct a BI-LSTM coupled prediction model. Then, the optimal parameter combination obtained in step C is substituted into the BI-LSTM coupled prediction model. Through the synergistic effect of the model's forget gate, input gate, cell state update and output gate, the karst spring flow is predicted.
2. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 1, characterized in that: In step A, the original dataset is preprocessed using MATLAB packages, including the following sub-steps: A10. Outlier Removal: The 3σ principle outlier detection method is used to identify and remove extreme outlier data points in the original dataset, retaining valid data. A20. Data Normalization: Perform IQR normalization on the dataset after removing extreme outliers; A30. Noise Suppression: For the dataset after IQR normalization, calculate the average value of all data within a preset window before and after each data point, and replace the original value of the corresponding data point with this average value to obtain the standardized dataset.
3. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 2, characterized in that: In the IQR normalization process, the IQR calculation formula is as follows: IQR = Q 3- Q 1, In the formula: Q 3 represents the data value at the 75th percentile after the original dataset is sorted in ascending order; Q 1 represents the data value at the 25th percentile after the original dataset is sorted in ascending order.
4. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 1, characterized in that: In step B, the specific mathematical expression and process of the convolution operation, stacking of multiple convolutional kernels, and average pooling operation of the CNN model are as follows: B10. The process data sequence X input to the CNN model is represented as: X = { x 1 , x 2 ,..., x T }, x T ∈R d , In the formula: T For time step, d Let R be the input feature dimension, and R be the real number field. B20. The comprehensive feature H obtained after stacking multiple sets of convolutional kernels from the output features of the convolutional layer is represented as follows: , In the formula: Concat It is K Each column vector (or row vector) is placed side-by-side by column (or by channel) to form a matrix; K The number of convolution kernels, H The final output of the convolution kernel; B30. The mathematical form of the average pooling operation is: , In the formula: This represents the output value of the j-th feature channel after average pooling. h j The pooling window size corresponding to the j-th feature channel. pool This is an average pooling operation; r This represents the size of the pooling window, and also the index range of positions within the pooling window; j This is the feature channel index, and also indicates that the current process is the first [number]th [channel] in that channel. j A pooled window, h jr Let be the feature value at position r within the j-th feature channel and the j-th pooling window.
5. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 1, characterized in that: In step C, the EI function, which serves as the acquisition function in the Bayesian optimization algorithm, has the following mathematical expression: , In the formula: α EI ( x ) is the EI function with respect to the independent variable x The function value at that point, i.e., the expected improvement value; Φ is the cumulative distribution function of the standard normal distribution. φ Let be the probability density function of the standard normal distribution. f ( x + This represents the best predicted value obtained by the surrogate model so far. μ ( x ) is the proxy model in the independent variable x The predicted mean at that location, σ ( x ) is the proxy model in the independent variable x The predicted standard deviation at point Z; Z is a standardized variable; in, ; The selection logic for the EI function is: when the independent variable... x Predicted mean at μ ( x (Better than the current best prediction) f ( x + When this happens, the EI function tends to select that independent variable. x The optimal predicted value is the best objective function value that the surrogate model believes to be so far. Wherein, the optimal objective function value f ( x + It is usually defined as: , In the formula: x i For the first in the training set i One sample point; X train The training dataset contains all evaluated datasets. x i ; μ ( x i )for x i The predicted mean at the location; Or when the independent variable x Predicted standard deviation at σ ( x When the value is large, the EI function will assign a value to the independent variable. x A certain selection weight.
6. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 5, characterized in that: In step C, the specific process of using a Bayesian optimization framework to solve for the optimal value of the log-likelihood function corresponding to the parameters through gradient descent is as follows: C10: The probability density functions corresponding to each parameter are: , In the formula: x 1 represents the set of the first t observed data points, i.e. x 1 , x 2 , ..., x t ; t is the number of observed data points, η is the noise term, y t Let be the observed value corresponding to the t-th observation data point. μ t For the surrogate model at the t-th observation data point x t The predicted mean at point T is the time step; C20: Due to the set of observation data points x 1:t and corresponding observations f(x 1:t ) Since all parameters are known, the parameter values that maximize the aforementioned probability density function are transformed into solving for the maximum value of the log-likelihood function, which is: , C30: Maximizing the above log-likelihood function is equivalent to minimizing its negative function, i.e., minimizing: , C40: Following the mechanism and process of solving the optimal value of probability density as described above, Bayesian optimization is performed on all other parameters in the BI-LSTM coupled prediction model and the surrogate model one by one, and finally the optimal values of all parameters are obtained, forming the optimal parameter combination.
7. The method for predicting karst spring flow based on Bayesian optimization algorithm and deep learning according to claim 1, characterized in that: In step D, the specific structure, mathematical expression, and function of the forget gate, input gate, cell state update, and output gate of the BI-LSTM coupled prediction model are as follows: Forget gate: Used to selectively filter historical information, deciding whether to discard or retain information from a previous moment; it uses an activation function. σ Process the hidden state from the previous moment h t-1 and current input x t Generate the gate value; The mathematical expression is: , In the formula: f t σ is the activation value of the forget gate, with a value range of [0,1]; σ is the Sigmoid activation function, with an output range of [0,1]. W f This is the weight matrix; h t-1 The hidden state from the previous moment; x t The input vector at the current time step; b f For the bias term of the forget gate; Input gate: Used to determine whether the current input information is updated to the cell state, providing new information to be stored in the memory unit; includes input gate activation value. i t It is a gated unit based on the activation function σ, which determines h t-1 and x t Has the information in the file been updated? And selectively update cell states to provide new information to be stored in memory units; The mathematical expression for candidate memory cell states is: , The mathematical expression for the input gate activation value is: , In the formula: represents the candidate memory cell state; tanh is the hyperbolic tangent activation function with an output range of [-1, 1]; W C , w i These are the weight matrices for the input gate and the candidate memory cell states, respectively. h t-1 The hidden state from the previous moment; x t This is the input vector at the current moment; i t The input gate activation value, which takes the range [0,1]. bi , b C These are the bias terms for the input gate and the candidate memory cell state, respectively; Cell state update: The cell state is updated by updating the cell state C from the previous time step. t-1 Output with the Forgot Gate f t Multiplication enables selective forgetting; then the input gate output is superimposed. i t With candidate state The product of these factors generates the updated cell state. The mathematical expression is: , In the formula: C t This represents the current state of the cell. Output Gates: Used to control the output of cell states, and the final model output information is determined through gating mechanisms; The mathematical expression for the output gate activation value is: , The mathematical expression for the hidden state at the current moment is: , In the formula: h t This represents the hidden state at the current moment, i.e., the output value of the BI-LSTM coupled prediction model; O t The activation value of the output gate, with a value range of [0,1]; W o This is the weight matrix of the output gate; b o This is the bias term for the output gate.
8. The method for predicting karst spring flow based on Bayesian optimization algorithm and deep learning according to any one of claims 1 to 7, characterized in that: It also includes a model performance evaluation step, which is performed after step D, using RMSE, MAE, and R... 2 The predictive performance of the BI-LSTM coupled prediction model is comprehensively evaluated using four evaluation metrics, including MAPE, residual analysis plots, Taylor plots, and time series comparison plots.
9. The karst spring flow prediction method based on Bayesian optimization algorithm and deep learning according to claim 8, characterized in that: The RMSE is the root mean square error, which measures the sample standard deviation of the difference between the predicted and observed values. The formula for calculating RMSE is: , The MAE is the Mean Absolute Error, which measures the average absolute error between the predicted and actual values. The formula for calculating MAE is: , The R 2 The coefficient of determination is used to measure how well the model fits the data. R 2 The calculation formula is: , The MAPE stands for Mean Absolute Percentage Error, which measures the average percentage difference between the predicted value and the actual observed value. The formula for calculating MAPE is: , In the formula: m For the number of samples, y i Let i be the true value of the i-th sample. Let be the predicted value for the i-th sample. This is the average of the true values of all samples.
10. A computer-readable storage medium having a computer program stored thereon, characterized in that: The computer program can be executed by one or more processors to implement the karst spring flow prediction method based on Bayesian optimization algorithm and deep learning as described in any one of claims 1 to 9.