Spring flow prediction method fusing physical mechanism and deep learning

Abstract

The present application relates to the technical field of hydrological data processing, and discloses a spring flow prediction method combining physical mechanism and deep learning, comprising the following steps: collecting input variables of a karst spring water system, and obtaining time delay input sequences after feature mode decomposition; inputting the time delay input sequences into a CNN-BiGRU-KAN hybrid neural model of a spring flow prediction model to obtain multi-scale enhanced features and a first spring flow prediction flow; inputting the time delay input sequences into a lightweight multi-layer perceptron structure to obtain a dissipation term; inputting the time delay input sequences and the dissipation term into a water balance equation of the spring flow prediction model to obtain a second spring flow prediction flow; and constructing a multi-objective joint loss function of the spring flow prediction model. The present application combines physical mechanism and a hybrid neural model based on deep neural network to predict spring flow, makes up for the problem of insufficient physical consistency of pure data models, and improves the accuracy of hydrological parameter inversion and flow prediction.

Description

Technical Field

[0001] This invention relates to the field of hydrological data processing technology, and in particular to a method for predicting spring flow that integrates physical mechanisms and deep learning. Background Technology

[0002] Karst springs, as an indispensable water resource, support the production and lives of approximately one-quarter of the world's population. They are not only the lifeline of regional water supply but also the cornerstone of maintaining unique spring ecosystems and socio-economic development. Currently, under the dual pressures of global climate change and intensive human activities, karst spring systems are facing unprecedented stress. On the one hand, climate change leads to frequent rainfall and extreme hydrological events, altering the system's original replenishment patterns and intensity, directly interfering with the seasonal changes of springs, and exacerbating the risk of springs drying up or experiencing sudden floods. On the other hand, with rapid socio-economic development, the impact of human activities is becoming increasingly prominent. The surge in water demand from industry, agriculture, and cities has led to excessive groundwater extraction, which is the most direct cause of the decline and even drying up of many karst springs. Drastic changes in land use patterns, such as urbanization, industrial and mining construction, and agricultural expansion, damage the vegetation and soil structure of the replenishment areas, weakening water conservation capacity and further exacerbating the depletion of spring resources.

[0003] Existing methods for predicting spring flow mainly fall into three categories, each with its own limitations: (1) Mathematical statistical regression (such as regression analysis, time series analysis, etc.), which predicts spring flow based on the mathematical statistical relationships of historical observation data. However, it ignores the physical mechanisms, does not delve into hydrogeological processes, and is difficult to reflect time lags, making it mostly suitable for short-term predictions. (2) Physical mechanism models, which take groundwater dynamics and water balance as the core, solve the flow control equations by numerical methods, and can realistically reproduce the physical processes of precipitation infiltration, groundwater movement and spring discharge, with high prediction accuracy. These include groundwater flow numerical simulation and multi-field coupling models. However, they have the following limitations: they rely on a large amount of observation data and detailed hydrogeological data. The generalization process of boundaries and hydrogeological conditions may artificially introduce uncertainty, which in turn affects the accuracy of the model and the reliability of the results. The modeling process is complex and computationally intensive, and it has poor applicability to areas lacking data. (3) Data-driven models: Data-driven shallow machine learning models (such as SVR, ANN, XGBoost, etc.) still have limitations when facing highly nonlinear complex data; Pure data-driven deep learning models have certain advantages in long-term series processing and nonlinear feature extraction, but due to the lack of clear physical mechanism constraints, their model structure and prediction results are usually insufficiently interpretable, and they are also prone to the problem of declining generalization ability under extreme climate or abnormal perturbation conditions. Summary of the Invention

[0004] The purpose of this invention is to solve the above-mentioned technical problems and provide a spring flow prediction method that integrates physical mechanisms and deep learning. By combining physical mechanisms with a hybrid neural model based on deep neural networks for spring flow prediction, it can, on the one hand, make up for the lack of physical consistency of pure data models to a certain extent and enhance the interpretability of the model; on the other hand, it can also give full play to the advantages of deep learning in complex nonlinear mapping and improve the accuracy of hydrological parameter inversion and flow prediction.

[0005] To achieve the above-mentioned objectives, the embodiments of the present invention provide the following technical solutions:

[0006] A spring flow prediction method integrating physical mechanisms and deep learning includes the following steps:

[0007] Step 1: Collect the input variables of the karst spring water system and obtain the time-delay input sequence through characteristic mode decomposition;

[0008] Step 2: Input the time-delayed input sequence into the CNN-BiGRU-KAN hybrid neural model of the spring water flow prediction model to obtain multi-scale enhanced features and the predicted flow of the first spring.

[0009] Step 3: After decomposing the feature modes, the time-delayed input sequence is input into a lightweight multilayer perceptron structure to obtain the dissipation term;

[0010] Step 4: Input the time-delay input sequence and the dissipation term into the water balance equation of the spring flow prediction model to obtain the predicted flow of the second spring.

[0011] Step 5: Construct a multi-objective joint loss function for the spring flow prediction model based on the predicted flow rates of the first and second springs.

[0012] In the above scheme, a spring flow prediction model that integrates the physical mechanism of water balance and deep learning is constructed. The core is a CNN-BiGRU-KAN hybrid neural model to achieve accurate fitting of spring flow and dissipation process. A physical constraint model is constructed based on the water balance equation. Through multi-objective joint loss function and joint training strategy, the data fitting accuracy and physical rationality are synergistically optimized. Finally, a spring flow prediction model that takes into account both prediction accuracy and physical rationality explanation is established.

[0013] Furthermore, in step 1, the input variables include precipitation P(t), daily average temperature T(t), groundwater level H(t), lateral runoff R(t), and factors related to the impact of human activities.

[0014] Furthermore, the step of obtaining the time-delayed input sequence after feature mode decomposition includes:

[0015] The precipitation P(t) is divided into K frequency bands, k=0,1,...,K; and K filters are generated simultaneously.

[0016] Maximizing the correlation kurtosis is used as the objective function for filter optimization:

[0017] ;

[0018] Among them, CM M (u k Let f be the objective function for the correlation kurtosis with shift order M; M is the shift order of the correlation kurtosis, and m is the periodic delay step, m=0,1,...,M; k (l) represents the coefficients of the k-th filter; u k (n) represents the original signal x(n) passing through the k-th filter f of length L. k The filtered modal signal; n represents the number of sampling points, n=1,2,...,N, where N is the total number of sampling points in the original signal; T s The fault cycle is expressed in units of sampling points; u k (n-mT s ) is the n-mT s The modal signal obtained by filtering the signal at each sampling point through the k-th filter of length L;

[0019] ;

[0020] Among them, CC(u p ,u q ) represents the correlation coefficient between the p-th mode signal and the q-th mode signal; u p (n) represents the modal signal output by the p-th filter; For N sampling points u p The mean of (n); u q (n) represents the modal signal output by the q-th filter; For N sampling points u q The mean of (n);

[0021] Discarding the mode with the smaller kurtosis among the two modes with the highest correlation coefficients and retaining the mode with the larger kurtosis, we obtain the precipitation decomposition result P`(t).

[0022] Furthermore, the CNN-BiGRU-KAN hybrid neural model in step 2 includes a deep neural network, a two-layer bidirectional gated recurrent unit, and a KAN network;

[0023] The time-delayed input sequence is fed into a deep neural network to obtain the local mutation feature F. CNN :

[0024] ;

[0025] in, For convolution operation; W1 and W2 are the weight matrices of the two convolutional layers, respectively; b1 and b2 are the bias terms of the two convolutional layers, respectively; ReLU is the ReLU activation function; BN is the batch normalization layer operation;

[0026] The time-delayed input sequence is input into a two-layer bidirectional gated recurrent unit to obtain the global time-series feature F. BiGRU :

[0027] ;

[0028] ;

[0029] ;

[0030] ;

[0031] Among them, z t To update the door; W is the Sigmoid activation function. z To update the gate weight matrix; h t-1 The hidden state of the previous moment; r t To reset the door; W r To reset the gate weight matrix; Let h be the candidate hidden state; W be the candidate hidden state weight matrix; h ... t The current hidden state; To perform element-wise multiplication, the hidden states of the forward-gated recurrent unit and the backward-gated recurrent unit at the same time are concatenated and output to obtain the global temporal feature F. BiGRU ;

[0032] Based on local mutation feature F CNN and global temporal features F BiGRU Obtain multi-scale enhanced features F fusion :

[0033] ;

[0034] Where ReLU is the ReLU activation function; Dropout is the Dropout regularization; W linear These are the weights for the linear mapping.

[0035] Furthermore, in step 2, the step of obtaining the predicted flow rate of the first spring includes: processing the multi-scale enhanced feature F... fusion Input a KAN network and output the predicted flow rate Q of the first spring. data (t).

[0036] Furthermore, in step 4, the water balance equation is:

[0037] ;

[0038] Among them, Q phys P(t) is the predicted flow rate of the second spring; P(t) is the original precipitation without characteristic mode decomposition; R(t) is the lateral runoff; E(t) is the actual evapotranspiration, and the quantitative relationship E(t) = g × T(t) is constructed using the temperature method, where T(t) is the daily average temperature, g is the evapotranspiration conversion coefficient, and D(t) is the dissipation term. ΔH is the change in water level per unit time, Δt is the change in water level per unit time; α is the rainfall infiltration recharge coefficient; C is the water storage coefficient.

[0039] Furthermore, the rainfall infiltration recharge coefficient α, water storage coefficient C, and evapotranspiration conversion coefficient g are specifically as follows:

[0040] ;

[0041] ;

[0042] ;

[0043] Among them, w α w C w g These are the trainable parameters.

[0044] Furthermore, in step 5, the multi-objective joint loss function Loss is:

[0045] ;

[0046] Among them, L MSE The data fitting loss is based on the predicted flow rate Q of the first spring. data The mean square error between (t) and the measured flow rate Q(t) is obtained; L phys The loss due to physical consistency is represented by the predicted flow rate Q from the second spring. phys The mean square error between the measured flow rate Q(t) and the actual flow rate Q(t) is obtained; This is the L2 regularization term, used to suppress model overfitting; This is the loss balance coefficient.

[0047] In the above scheme, combining physical mechanisms with deep learning models for spring flow prediction can, on the one hand, compensate for the lack of physical consistency of pure data models to a certain extent and enhance the interpretability of the models; on the other hand, it can also leverage the advantages of deep learning in complex nonlinear mappings to improve the accuracy of hydrological parameter inversion and flow prediction.

[0048] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0049] (1) Closed-loop coupled physical information fusion architecture: The dissipation term D(t) adaptively generated by a lightweight multilayer perceptron structure is explicitly embedded into the water balance equation of the closed watershed. This forms a closed-loop mechanism of "prediction-physical feedback-parameter correction", which makes up for the shortcomings of traditional physical models in quantifying unmonitored leakage and artificial extraction.

[0050] (2) Joint inversion of physical parameters and network weights: The process of independent parameter calibration of traditional hydrological models is abandoned. The hydrogeological parameters (α, C, g) and neural network weights are placed under the joint loss function for end-to-end synchronous training. The differentiable mapping function ensures that the parameters always have hydrogeological significance in the iteration.

[0051] (3) Fusion of KAN network and multi-scale feature extraction: In the field of hydrological prediction, KAN network is integrated and a learnable spline activation function is used to replace the traditional fully connected layer, which fundamentally reduces the model's dependence on the linear weight matrix. Combined with CNN and BiGRU, the ability to analyze extreme nonlinear features is greatly improved. Attached Figure Description

[0052] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0053] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0054] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0055] It should be noted that similar reference numerals and letters in the following figures denote similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, the terms "first," "second," etc., are used only for distinguishing descriptions and should not be construed as indicating or implying relative importance, or suggesting any such actual relationship or order between these entities or operations. Additionally, the terms "connected," "linked," etc., can refer to a direct connection between elements, components, modules, etc., or an indirect connection via other elements, components, modules, etc.

[0056] Example 1:

[0057] This invention is achieved through the following technical solution: a method for predicting spring flow that integrates physical mechanisms and deep learning, such as... Figure 1 As shown, it includes the following steps:

[0058] Step 1: Collect the input variables of the karst spring water system and obtain the time-delay input sequence through characteristic mode decomposition.

[0059] This scheme standardizes the problem of predicting the flow rate of karst spring systems into a multivariate hydrological time series prediction problem. Let the discrete time scale be t = 1, 2, ..., T, and the input variable vector obtained from observations at any time t be represented as follows: , Let m be the set of real numbers, m be the dimension, and x be the input variable vector. t The parameters include precipitation P(t), daily average temperature T(t), groundwater level H(t), and lateral runoff R(t), while also needing to observe factors related to human activities. The target variable for prediction is spring water flow. A time-delayed input sequence is constructed based on a sliding time window, expressed as:

[0060] ;

[0061] ;

[0062] Among them, X t X is a time-delayed input sequence; Len is the length of the sliding time window, used to characterize the lag effects of precipitation infiltration, groundwater recharge, and aquifer storage and release processes in karst spring systems. t As input to the spring water flow prediction model, it enables accurate prediction of spring water flow at current and future times.

[0063] To improve the ability of the subsequent spring flow prediction model to decouple features from precipitation P(t), this scheme uses the Eigenmode Decomposition (FMD) algorithm to analyze the time-delayed input sequence X. tPrecipitation P(t) is preprocessed. The Eigenmode Decomposition algorithm adopts a non-recursive decomposition method, selecting different models by initializing the adaptive finite impulse response (FIR) filter and updating the filter coefficients, while considering the impulsiveness and periodicity of the signal, thus achieving good results in precipitation feature extraction.

[0064] (a) Filter expression.

[0065] In signal decomposition tasks, adaptive finite impulse response (FIR) filters selectively enhance or suppress specific frequencies or time domains, thus "filtering" specific signals and aiding in the decomposition process. An FIR filter is a digital filter implemented in the time domain through convolution; its output depends only on the current and past input values, contains no feedback loop, and its expression is:

[0066] ;

[0067] Where y(n) is the output of the filter; h(l) is the filter coefficient; x(nl) is the sampled value of the input signal x[n] after passing through l sampling points; L is the filter length (corresponding to the filter order L-1), l=0,1,...,L-1.

[0068] (ii) Filter update and period estimation.

[0069] The Eigenmode Decomposition (EMD) algorithm improves upon traditional filters by adding coefficient adjustments. It uses the Hanning window (cutoff frequency) as the initialization method to initialize the adaptive finite impulse response (FIR) filter (hereinafter referred to as the filter) and iteratively optimizes the filter coefficients to make the filtered signal infinitely close to the objective function. During initialization, the original signal bandwidth is first uniformly divided into K frequency bands, and the upper and lower cutoff frequencies of these bands are assigned as follows:

[0070] ;

[0071] Where k is the frequency band number, k=0,1,...,K-1 (the range of the 0th frequency band is [0,f)). s / (2K)]); f s f is the sampling frequency of the original signal (i.e., precipitation P(t)); k,low f is the lower cutoff frequency corresponding to the k-th frequency band; k,high Let L be the upper cutoff frequency corresponding to the k-th frequency band. The above formula can generate K filters with different cutoff frequencies, each with a length of L. Using a Hanning window, K evenly distributed filters are generated, forming a filter bank covering the entire frequency band.

[0072] The Eigenmode Decomposition (EMD) algorithm uses maximizing the correlation kurtosis (CK) as the objective function for filter optimization, transforming EMD theory into a constraint problem to be solved. The constraint equations are:

[0073] ;

[0074] Among them, CM M (u k Let f be the objective function for the correlation kurtosis with shift order M; M is the shift order of the correlation kurtosis, and m is the number of periodic delay steps (loop variable), m=0,1,...,M; k (l) represents the coefficients of the k-th filter; u k (n) represents the original signal x(n) passing through the k-th filter f of length L. k The filtered modal signal; n represents the number of sampling points, n=1,2,...,N, where N is the total number of sampling points in the original signal; T s The fault cycle is expressed in units of sampling points; u k (n-mT s ) is the n-mT s The modal signal obtained by filtering the signal at each sampling point through the k-th filter of length L. This is achieved by finding a set of filter coefficients {f}. k (l)}, making the objective function CM M (u k (To reach maximum, molecules) By u k (n) at different time points n-mT s The sum of the squares of the products reflects the signal u k (n) Correlation across multiple delay periods; denominator Used to normalize molecules, making comparisons between signals at different sampling points more reasonable.

[0075] ;

[0076] Where x(n-l+1) is the signal at the (n-l+1)th sampling point. This formula obtains the modal signal u output by the kth filter by summing the products of the coefficients of different l with the corresponding sampling points. k (n).

[0077] Since the signal period is difficult to estimate accurately in practice, the Improved Maximum Correlated Kurtosis Deconvolution (IMCKD) method is used to estimate the signal period from the measured signal based on autocorrelation theory. Specifically, the autocorrelation spectrum will produce a local maximum at the periodic location, and the point where the autocorrelation spectrum reaches a local maximum after crossing zero is selected as the estimated period. As the filter is updated, the estimated period will be continuously updated and become more accurate.

[0078] (iii) Modal selection.

[0079] To avoid generating a large number of redundant modes by updating all initialization filters, the Eigenmode Decomposition (EMD) algorithm proposes a mode selection strategy based on the correlation coefficient (CC). Eigenmode Decomposition selects the mode with the largest kurtosis (CK). To avoid multiple modes containing the same component, the two modes with the largest correlation coefficients are locked, defined as follows:

[0080] ;

[0081] Among them, CC(u p ,u q Let be the correlation coefficient between the p-th mode signal and the q-th mode signal, where q = 0, 1, ..., K-1, and p = 0, 1, ..., K-1; u p (n) represents the modal signal output by the p-th filter; For N sampling points u p The mean of (n); u q (n) represents the modal signal output by the q-th filter; For N sampling points u q The mean of (n) is used. The resulting correlation coefficient ranges from -1 to 1. The closer to 1, the stronger the correlation between the two modes, meaning they contain more identical components. The mode with the lower kurtosis (CK) among the two modes with the highest correlation coefficient is discarded, and the mode with the higher kurtosis (CK) is retained. This eliminates mode aliasing and redundant modes, retaining the mode with more information.

[0082] The Eigenmode Decomposition (EMD) algorithm considers both the impulsiveness and periodicity of the signal, exhibiting stronger robustness against interference and noise. It extracts decomposed modes using an adaptive finite impulse response (FIR) filter, which is not limited by filter shape, bandwidth, or center frequency, resulting in a more thorough decomposition. In this scheme, Eigenmode Decomposition is used to perform mode decomposition on precipitation P(t), yielding the decomposed result P'(t). The original time-delayed input sequence X is then replaced with the decomposed result P'(t). t The precipitation P(t) in the middle is used to obtain a new time-delay input sequence. This serves as the input for subsequent spring flow prediction models. Alternatively, the decomposition result P`(t) can be directly added to the original time-delay input sequence X. t In this process, a new time-delay input sequence is obtained. .

[0083] The spring water flow prediction model includes a first branch and a second branch; the first branch is a CNN-BiGRU-KAN hybrid neural model that takes a time-delayed input sequence. After input, output the predicted flow rate Q of the first spring. data (t); the second branch is the water balance equation, which inputs the time-delay sequence. After input, output the predicted flow rate Q of the second spring. phys (t); Finally, the spring water flow prediction model uses a multi-objective joint loss function, based on the Q values ​​output by the two branches respectively. data (t) and Q phys (t) The optimization objectives are to unify data fitting accuracy, physical constraint consistency and model generalization ability.

[0084] Step 2: Input the time-delayed input sequence into the CNN-BiGRU-KAN hybrid neural model of the spring water flow prediction model to obtain multi-scale enhanced features and the predicted flow of the first spring.

[0085] The CNN-BiGRU-KAN hybrid neural model includes a deep neural network (CNN), a two-layer bi-directional gated recurrent unit (BiGRU), a KAN network, and a time-delayed input sequence. The local mutation feature F output by the deep neural network is obtained by inputting the deep neural network and the two-layer bidirectional gated recurrent unit respectively. CNN The global timing characteristics F of the output of the two-layer bidirectional gated cyclic unit BiGRU After fusion, multi-scale enhancement features F are formed. fusion Finally, multi-scale enhancement features F fusion Input a KAN network, and the KAN network outputs the predicted flow rate Q of the first spring. data (t).

[0086] This scheme employs a deep neural network (CNN) with two stacked one-dimensional convolutional layers (Conv1D) to extract local abrupt change features, accurately capturing the rapid impact of short-term factors such as sudden precipitation changes, abrupt temperature fluctuations, and human disturbances on spring water flow. The first convolutional layer uses 64 kernels of size 3, and the second convolutional layer uses 128 kernels of size 3. Both convolutional layers use the same padding method to ensure the time-delayed input sequence is properly padded. The length remains consistent before and after the convolution operation, preserving the complete temporal location information of the hydrological sequence. After the convolution operation, a batch normalization (BN) layer, a ReLU activation function, and a dropout regularization layer are sequentially applied. The batch normalization layer alleviates internal covariate bias during CNN training, improving training stability; the ReLU activation function introduces non-linear mapping capabilities into the CNN; and the regularization layer suppresses overfitting by randomly deactivating some neurons. The output of a deep neural network (CNN) can be represented as:

[0087] ;

[0088] Among them, F CNN These are local mutation features output by a deep neural network. For convolution operations; W1 and W2 are the weight matrices of the two convolutional layers, respectively; b1 and b2 are the bias terms of the two convolutional layers, respectively; ReLU is the ReLU activation function; BN is the batch normalization layer operation. After extraction by the two convolutional layers of the CNN, the output is a 128-dimensional local mutation feature F. CNN This provides high-quality shallow features for subsequent long-term time-dependent encoding.

[0089] To efficiently capture the time-delayed input sequence of karst spring systems To address the long-short-term dependence and bidirectional lag effect, this invention employs a two-layer bidirectional gated recurrent unit (BiGRU) for timing encoding. The BiGRU includes both a forward GRU and a backward GRU, enabling it to encode time-delayed input sequences. BiGRU uses parallel encoding in both the beginning and end directions to fully explore the correlation between historical and future states. Internally, BiGRU controls the retention and forgetting of information through update and reset gates. The calculation expression is:

[0090] ;

[0091] ;

[0092] ;

[0093] ;

[0094] Among them, z t This is used to update the gate and control the proportion of information retained from the previous hidden state. W is the Sigmoid activation function. z To update the gate weight matrix; h t-1 The hidden state of the previous moment; r t To reset the gate, used to control the degree of forgetting of historical information; W r To reset the gate weight matrix; Let h be the candidate hidden state; W be the candidate hidden state weight matrix; h ... t The current hidden state; Element-wise multiplication is performed. The hidden states of the forward GRU and the backward GRU at the same time step are concatenated to obtain the 100-dimensional global temporal feature F. BiGRU It fully characterizes the long-term evolution of karst spring systems.

[0095] To mitigate the gradient degradation problem in deep networks and preserve effective features from shallow layers, this invention constructs a cross-layer residual connection structure from CNN to BiGRU. The 128-dimensional local mutation features F output by the CNN are processed through a linear fully connected layer. CNN Mapping to a 100-dimensional feature space yields residual features, which are then mapped to the global temporal features F output by the BiGRU. BiGRU Dimensionality is consistent. The mapped residual features are compared with the global temporal features F. BiGRU The elements are added together and fused. The fused result is then processed by Dropout regularization and ReLU activation function to obtain the multi-scale enhanced feature F that integrates local mutation features and global temporal features. fusion The calculation expression is:

[0096] ;

[0097] Where ReLU is the ReLU activation function; Dropout is the Dropout regularization; W linear These are the weights for the linear mapping.

[0098] This invention abandons the traditional fully connected layer with a fixed activation function and enhances feature F at multiple scales. fusion As input, a KAN (Kolmogorov-Arnold) network is used as the first branch. Traditional linear weights are replaced with learnable one-dimensional B-spline basis functions. High-precision approximation of strongly nonlinear relationships is achieved by dynamically adjusting the spline grid parameters of the B-spline basis functions. Multi-scale enhancement features F... fusion Input a KAN network and output the predicted flow rate Q of the first spring. data (t). Compared to traditional fully connected layers, KAN networks can more accurately fit the complex nonlinear relationship of precipitation-replenishment-discharge in karst spring systems, while also possessing stronger interpretability, which is the core guarantee for the model's prediction accuracy.

[0099] Step 3: After decomposing the feature modes, the time-delayed input sequence is input into the lightweight multilayer perceptron structure to obtain the dissipation term.

[0100] To quantify hidden water loss processes that are difficult to parameterize using physical formulas, such as those caused by artificial extraction, unmonitored seepage, and water consumption within a watershed, this invention independently sets up a water loss term fitting branch. It uses precipitation P'(t) after eigenmode decomposition, daily average temperature T(t), groundwater level H(t), lateral runoff R(t), and relevant factors related to human activities to fit the loss term. The lightweight multilayer perceptron (MLP) structure includes a linear transformation layer, a ReLU activation function, a Dropout regularization layer, and a linear output layer connected in sequence. The output feature is then forcibly rectified using the ReLU activation function at the end of the linear output layer to ensure the output result is non-negative, resulting in a physically meaningful water loss term D(t). This water loss term D(t) participates in the subsequent physical calculation of the water balance equation and is simultaneously optimized during backpropagation, serving as a core bridging variable connecting CNN, BiGRU, and the water balance equation.

[0101] Step 4: Input the time-delay input sequence and the dissipation term into the water balance equation of the spring flow prediction model to obtain the predicted flow of the second spring.

[0102] This invention introduces the dissipation term D(t) obtained from the time-delayed input sequence after eigenmode decomposition into the water balance equation. This allows the CNN-BiGRU-KAN hybrid neural model to be embedded into the water balance equation, constructing core physical constraints and ensuring that the output of the spring flow prediction model strictly satisfies hydrophysical laws. The water balance equation is as follows:

[0103] ;

[0104] Among them, Q phys P(t) is the predicted flow rate of the second spring; P(t) is the original precipitation without characteristic mode decomposition; R(t) is the lateral runoff; E(t) is the actual evapotranspiration, and the quantitative relationship E(t) = g × T(t) is constructed using the temperature method, where T(t) is the daily average temperature, g is the evapotranspiration conversion coefficient, and D(t) is the dissipation term. ΔH is the change in water level per unit time, Δt is the change in water level per unit time; α is the rainfall infiltration recharge coefficient; C is the water storage coefficient.

[0105] To ensure that the physical parameters always satisfy the hydrophysical meaning during the training of the spring flow prediction model, and to avoid problems such as gradient discontinuity and training oscillation caused by hard truncation, this invention adopts a parameter constraint mechanism of "initial trainable parameters + differentiable smooth mapping" to obtain α, C, and g. First, the trainable parameters w are randomly initialized. α w C w g Using it as the network weight input optimizer, it can directly receive the gradient from backpropagation and synchronously update it to participate in the computation iteration, adaptively determine the convergence condition, and after iteration, transfer w through a mapping relationship. αw C w g The optimal values ​​are converted into physical parameters α, C, and g for output. Specifically, the trainable parameters w are converted into g using a differentiable function. α w C w g The mappings are to physical parameters α, C, and g, which have physical meaning, and the mapping relationships are as follows:

[0106] ;

[0107] ;

[0108] ;

[0109] Specifically, the rainfall infiltration recharge coefficient α is mapped to the (0,1) interval through a Sigmoid activation function to ensure that the effective infiltration amount does not exceed the original rainfall P(t); the water storage coefficient C and the evapotranspiration conversion coefficient g are mapped to the (0,+∞) interval through an exponential function to ensure that the parameters are always positive. Since the mapping relationship is differentiable throughout, the gradient can be smoothly transferred to the trainable parameter w during training. α w C w g Ultimately, the physical parameters α, C, and g are inverted without the need for additional independent calibration. Substituting the optimal physical parameters α, C, and g into the water balance equation yields the predicted flow rate Q of the second spring. phys (t).

[0110] Step 5: Construct a multi-objective joint loss function for the spring flow prediction model based on the predicted flow rates of the first and second springs.

[0111] This invention constructs a multi-objective joint loss function for a spring water flow prediction model, unifying the optimization objectives of data fitting accuracy, physical constraint consistency, and model generalization ability. The complete optimization process is as follows:

[0112] (1) Calculate the data fitting loss. The predicted flow rate Q of the first spring, output by the KAN network, is used as the basis. data The mean square error (MSE) between the measured flow rate Q(t) and the actual flow rate Q(t) is used as the data fitting loss L. MSE The spring flow prediction model is guided to learn the nonlinear mapping relationship between hydrological and meteorological factors and spring flow:

[0113] ;

[0114] Among them, Q i Let be the measured flow rate of the i-th sample; Predict the flow rate of the first spring for the i-th sample; is the mean squared error loss function; N is the number of samples.

[0115] (2) Calculate the physical consistency loss. Calculate the predicted flow rate Q of the second spring. phys The mean square error (MSE) between the measured flow rate Q(t) and the actual flow rate Q(t) is calculated using a similar formula to L. MSE ) as the physical consistency loss for residual construction L phys Simultaneously, the ReLU activation function constraint method is used to achieve the Q... phys (t) Non-negative control: The computational flow output by the physical module is mapped through the ReLU activation function, automatically filtering all negative results, so that the output flow is always greater than or equal to 0, and the non-negative constraint of the physical computational flow is efficiently realized.

[0116] (3) Construct a multi-objective joint loss function. The data fitting loss L... MSE Physical consistency loss L phys The weighted summation with the L2 regularization term yields the multi-objective joint loss function Loss:

[0117] ;

[0118] in, This is the L2 regularization term, used to suppress model overfitting; This is the loss balancing coefficient, used to adjust the relative weights of the various loss terms. .

[0119] (4) End-to-end backpropagation optimization. The Adam optimizer is used to perform iterative optimization within a unified computation graph. The gradient of the multi-objective joint loss function Loss is simultaneously backpropagated to all network weights and trainable parameters w. α w C w g Simultaneously update feature extraction capabilities and physical parameter values. During training, data fitting loss drives the model to fit the observed data, while physical consistency loss guides the model to conform to hydrological laws. Together, they promote the model to converge towards "high accuracy + high physical rationality".

[0120] To improve the generalization performance and training stability of the model, this invention employs an improved Artificial Goat Population Optimization (AGTO) algorithm to perform global optimization on key hyperparameters such as sliding time window length, initial learning rate, number of training iterations, and batch size, using the minimum set mean squared error (MSE) as the fitness function. Building upon the traditional GTO algorithm, the AGTO algorithm introduces an Elite Backward Learning (EOBL) strategy during the initialization phase to improve the quality and diversity of the initial population. During the position update phase, it introduces the Cauchy Inverse Accumulation (CICD) operator and the Tangential Flight (TFO) operator to improve the position update mechanisms of "following silverbacks" and "competing for adult females," balancing global exploration and local exploitation capabilities. This effectively avoids premature convergence and local optima, and the optimal hyperparameter combination is obtained through global optimization for joint model training.

[0121] Hydrological verification experiments were conducted using the aforementioned hybrid model and optimal parameters. The model demonstrated exceptional fitting and generalization performance in predicting spring flow rates. After training and optimization, the final predicted flow rate Q of the first spring was determined. data (t) is used as the final prediction result. Specific model evaluation metrics are: Nash efficiency coefficient (NSE) as high as 0.971, mean absolute percentage error (MAPE) as low as 5.251%, and root mean square error (RMSE) as low as 0.239. Experimental results fully demonstrate that this invention achieves high-precision fitting of complex nonlinear mapping relationships in karst hydrological systems while effectively suppressing data noise interference, solving the problem of performance degradation of single deep learning models under data scarcity or extreme hydrological events.

[0122] The above description is merely a specific 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 technical scope 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 spring flow prediction method integrating physical mechanisms and deep learning, characterized in that, Includes the following steps: Step 1: Collect the input variables of the karst spring water system and obtain the time-delay input sequence through characteristic mode decomposition; Step 2: Input the time-delayed input sequence into the CNN-BiGRU-KAN hybrid neural model of the spring water flow prediction model to obtain multi-scale enhanced features and the predicted flow of the first spring. Step 3: After decomposing the feature modes, the time-delayed input sequence is input into a lightweight multilayer perceptron structure to obtain the dissipation term; Step 4: Input the time-delay input sequence and the dissipation term into the water balance equation of the spring flow prediction model to obtain the predicted flow of the second spring. Step 5: Construct a multi-objective joint loss function for the spring flow prediction model based on the predicted flow rates of the first and second springs.

2. The spring water flow prediction method integrating physical mechanisms and deep learning according to claim 1, characterized in that, In step 1, the input variables include precipitation P(t), daily average temperature T(t), groundwater level H(t), lateral runoff R(t), and factors related to the impact of human activities.

3. The spring flow prediction method integrating physical mechanisms and deep learning according to claim 1, characterized in that, The step of obtaining the time-delayed input sequence after feature mode decomposition includes: The precipitation P(t) is divided into K frequency bands, k=0,1,...,K; and K filters are generated simultaneously. Maximizing the correlation kurtosis is used as the objective function for filter optimization: ； Among them, CM M (u k Let f be the objective function for the correlation kurtosis with shift order M; M is the shift order of the correlation kurtosis, and m is the periodic delay step, m=0,1,...,M; k (l) represents the coefficients of the k-th filter; u k (n) represents the original signal x(n) passing through the k-th filter f of length L. k The filtered modal signal; n represents the number of sampling points, n=1,2,...,N, where N is the total number of sampling points in the original signal; T s The fault cycle is expressed in units of sampling points; u k (n-mT s ) is the n-mT s The modal signal obtained by filtering the signal at each sampling point through the k-th filter of length L; ； Among them, CC(u p ,u q ) represents the correlation coefficient between the p-th mode signal and the q-th mode signal; u p (n) represents the modal signal output by the p-th filter; For N sampling points u p The mean of (n); u q (n) represents the modal signal output by the q-th filter; For N sampling points u q The mean of (n); Discarding the mode with the smaller kurtosis among the two modes with the highest correlation coefficients and retaining the mode with the larger kurtosis, we obtain the precipitation decomposition result P`(t).

4. The spring water flow prediction method integrating physical mechanisms and deep learning according to claim 1, characterized in that, The CNN-BiGRU-KAN hybrid neural model in step 2 includes a deep neural network, a two-layer bidirectional gated recurrent unit, and a KAN network. The time-delayed input sequence is fed into a deep neural network to obtain the local mutation feature F. CNN : ； in, For convolution operation; W1 and W2 are the weight matrices of the two convolutional layers, respectively; b1 and b2 are the bias terms of the two convolutional layers, respectively; ReLU is the ReLU activation function; BN is the batch normalization layer operation; The time-delayed input sequence is input into a two-layer bidirectional gated recurrent unit to obtain the global time-series feature F. BiGRU : ； ； ； ； Among them, z t To update the door; W is the Sigmoid activation function. z To update the gate weight matrix; h t-1 The hidden state of the previous moment; r t To reset the door; W r To reset the gate weight matrix; Let h be the candidate hidden state; W be the candidate hidden state weight matrix; h ... t The current hidden state; To perform element-wise multiplication, the hidden states of the forward-gated recurrent unit and the backward-gated recurrent unit at the same time are concatenated and output to obtain the global temporal feature F. BiGRU ; Based on local mutation feature F CNN and global temporal features F BiGRU Obtain multi-scale enhanced features F fusion : ； Where ReLU is the ReLU activation function; Dropout is the Dropout regularization; W linear These are the weights for the linear mapping.

5. The spring flow prediction method integrating physical mechanisms and deep learning according to claim 4, characterized in that, Step 2, the step of obtaining the predicted flow rate of the first spring, includes: processing the multi-scale enhanced features F... fusion Input a KAN network and output the predicted flow rate Q of the first spring. data (t).

6. The spring water flow prediction method integrating physical mechanisms and deep learning according to claim 1, characterized in that, In step 4, the water balance equation is: ； Among them, Q phys P(t) is the predicted flow rate of the second spring; P(t) is the original precipitation without characteristic mode decomposition; R(t) is the lateral runoff; E(t) is the actual evapotranspiration, and the quantitative relationship E(t) = g × T(t) is constructed using the temperature method, where T(t) is the daily average temperature, g is the evapotranspiration conversion coefficient, and D(t) is the dissipation term. ΔH is the change in water level per unit time, Δt is the change in water level per unit time; α is the rainfall infiltration recharge coefficient; C is the water storage coefficient.

7. The spring water flow prediction method integrating physical mechanisms and deep learning according to claim 6, characterized in that, The rainfall infiltration recharge coefficient α, water storage coefficient C, and evapotranspiration conversion coefficient g are specifically as follows: ； ； ； Among them, w α w C w g These are the trainable parameters.

8. The spring flow prediction method integrating physical mechanisms and deep learning according to claim 1, characterized in that, In step 5, the multi-objective joint loss function Loss is: ； Among them, L MSE The data fitting loss is based on the predicted flow rate Q of the first spring. data The mean square error between (t) and the measured flow rate Q(t) is obtained; L phys The loss due to physical consistency is represented by the predicted flow rate Q from the second spring. phys The mean square error between the measured flow rate Q(t) and the actual flow rate Q(t) is obtained; This is the L2 regularization term, used to suppress model overfitting; This is the loss balance coefficient.

Images