Hourly runoff forecasting method based on CEEMDAN-VMD-TFLS

Abstract

The application discloses a kind of CEEMDAN-VMD-TFLS-based hour runoff prediction method, comprising: the collected hour scale historical runoff sequence and rainfall sequence are handled, obtain standardized measured hour runoff data and rainfall data;Standardized measured hour runoff data is decomposed using CEEMDAN, and IMF component and k-order residual component are obtained, the fuzzy entropy of each IMF component is calculated, and based on the fuzzy entropy characteristics of each IMF component, it is reconstructed into high, medium and low frequency sequence using spectral clustering method, and the high frequency sequence is secondarily decomposed using VMD;Key input factors are screened using partial autocorrelation function combined with feature importance analysis and aligned with time series to form an input feature set, a TFLS hybrid prediction model is constructed, and a runoff prediction result is obtained by inputting the predicted features;The runoff prediction result is superimposed and reconstructed to obtain the hour runoff prediction value.The application can effectively weaken the influence of non-stationary and multi-scale noise of runoff sequence on modeling, improve the accuracy and stability of hour runoff prediction.

Description

Technical Field

[0001] This invention relates to the field of runoff forecasting technology, and in particular to an hourly runoff forecasting method based on CEEMDAN-VMD-TFLS. Background Technology

[0002] The uncertainty of runoff processes directly affects the scientific validity and safety of major decisions such as flood control and disaster reduction, optimal allocation of water resources, and ecological water use security. Therefore, developing high-precision, stable, and reliable runoff forecasts has become a core research direction in the field of hydrology and water resources. In recent years, climate change and intensified human activities have significantly increased the randomness and abruptness of watershed hydrological processes. Runoff sequences exhibit nonlinear, highly uncertain, and complex non-stationary characteristics, placing increasingly higher demands on the accuracy and scale of runoff forecasts. Conventional time series methods are prone to accumulating forecast biases when processing such signals, making it difficult to meet practical application needs. Furthermore, existing forecasting methods suffer from low accuracy and short lead times in hourly-scale forecasts, failing to meet the needs of practical engineering projects such as reservoir flood control scheduling. To address this, researchers have gradually introduced signal processing and intelligent modeling techniques, extracting core variation components at different frequencies of runoff sequences through multi-scale decomposition and reconstruction, providing technical support for improving forecast accuracy.

[0003] Existing hydrological forecasting methods can be broadly classified into process-driven and data-driven methods. However, with the frequent occurrence of extreme weather and rapid changes in the underlying surface, the non-stationarity and multi-scale characteristics of runoff sequences have become more prominent. Traditional single models or simple neural networks are difficult to adapt to such highly nonlinear signals, especially in multi-step predictions at the hourly scale, where the accumulation of errors is significant and the limitations are obvious.

[0004] In hourly runoff forecasting, the sequence is superimposed with various high-frequency disturbances such as sudden weather changes, human intervention, and local confluence, resulting in a low signal-to-noise ratio and strong volatility. If hourly runoff data is not properly decomposed and frequency-banded, and is directly input into a single deep learning model, the model may struggle to simultaneously account for both high-frequency short-term fluctuations and low-frequency long-term trends. This leads to biases in peak value characterization and baseflow simulation, impacting the reliability of flood warnings and real-time management. Therefore, it is crucial to explore the evolution patterns and intrinsic relationships of runoff sequences on an hourly timescale through the fusion of modal decomposition and deep learning techniques. Developing an hourly runoff forecasting method and system based on the fusion of runoff signal decomposition and a hybrid deep learning model is of significant practical importance for improving the timeliness of basin flood forecasts, optimizing the scientific nature of reservoir management decisions, and enhancing the level of refined water resource management. Summary of the Invention

[0005] Purpose of the invention: The purpose of this invention is to provide an hourly runoff forecasting method based on CEEMDAN-VMD-TFLS, which solves the problems in hourly runoff forecasting. Due to the significant non-stationarity and multi-scale characteristics of runoff sequences caused by multiple sources of high-frequency disturbances such as weather changes, artificial scheduling, and local confluence, single deep learning models cannot simultaneously take into account high-frequency short-term fluctuations and low-frequency long-term trends. The accumulation of multi-step prediction errors is serious, and the accuracy of peak characterization and baseflow simulation is insufficient.

[0006] Technical solution: The hourly runoff forecasting method based on CEEMD-TFLS described in this invention includes:

[0007] S1. Collect hourly measured flow data from the watershed forecast section and upstream hydrological stations, as well as hourly precipitation data from rain gauge stations within the watershed, to construct a complete hourly historical runoff and precipitation sequence; perform missing value imputation, outlier handling, and normalization on the hourly historical runoff and precipitation sequences to obtain standardized measured hourly runoff and precipitation data.

[0008] S2. The standardized measured hourly runoff data are decomposed using complete ensemble empirical mode decomposition to obtain IMF components and k-th order residual components; the fuzzy entropy of each IMF component is calculated, and based on the fuzzy entropy characteristics of each IMF component, it is reconstructed into high, medium and low frequency sequences using spectral clustering method, and the high frequency sequence is subjected to VMD secondary decomposition.

[0009] S3. Use partial autocorrelation function combined with feature importance analysis to screen key input factors and align them with time series to form an input feature set; construct a TFLS hybrid forecasting model, input the features to be forecasted to obtain runoff prediction results;

[0010] S4. Determine the basic weights based on the historical energy proportions of each modal component, and adjust the weights by combining the correlation strength between rainfall and runoff. Then, overlay and reconstruct the runoff prediction results to obtain the hourly runoff forecast values.

[0011] Furthermore, the standardized measured hourly runoff data are decomposed using complete ensemble empirical mode decomposition, including:

[0012] (21) Gaussian white noise was added to the standardized measured hourly runoff data to form the overall decomposed signal;

[0013] (22) The overall average value is calculated by repeating the EMD decomposition process I times, and it is defined as the first IMF component of the original signal, where I is the number of times the EMD decomposition is repeated;

[0014] (23) Calculate the first-order residual based on the standardized measured hourly runoff data with added Gaussian white noise and the first IMF component of the original signal;

[0015] (24) After adding specific noise to the first-order residual, perform EMD decomposition again and calculate the average value to obtain the second IMF component;

[0016] (25) Calculate the k-th order residual;

[0017] (26) After adding specific noise to the k-th residual and performing EMD decomposition again, the average value is obtained to obtain the (k+1)th IMF component;

[0018] (27) Repeat steps (25)-(26) until the residuals cannot be decomposed, and obtain the final residuals;

[0019] (28) The standardized measured hourly runoff data after decomposition are obtained based on the (k+1)th IMF component and the final residual.

[0020] Furthermore, calculating the fuzzy entropy of each IMF component includes:

[0021] (31) Take the phase space dimension m, and construct N-m+1 m-dimensional vectors X. m，i N represents the length of the IMF component time series of the fuzzy entropy to be calculated;

[0022] (32) Define vector X m，i X m，j absolute distance It is the maximum value of the difference between its corresponding elements;

[0023] (33) Introduce the fuzzy function μ(d ij m Define X(n, r) m，i X m，j Similarity The expression is:

[0024]

[0025] In the formula, r is the tolerance width and n is the tolerance gradient;

[0026] (34) Define the function :

[0027]

[0028] (35) Change the phase space dimension to m+1, repeat steps (31)-(34), and obtain the function. :

[0029]

[0030] (36) When the sample has N data points, the fuzzy entropy is defined as:

[0031]

[0032] Furthermore, the VMD secondary decomposition of the high-frequency sequence includes:

[0033] (41) By using the Intrinsic Mode Function (IMF) k Converted into a bandwidth-limited amplitude-frequency modulation signal μ k The minimum sum of bandwidths estimated by intrinsic mode functions is used as the constraint condition for decomposition:

[0034]

[0035] And it satisfies:

[0036]

[0037] In the formula, To standardize the measured hourly runoff data, k is the number of modal components obtained from the decomposition; ω is the k-th modal component; k Let be the center angular frequency of the k-th modal component; K is the total number of modal components; It is the Dirac function; This represents the partial derivative operation with respect to the time variable t; It is a complex exponential modulation term;

[0038] (42) Introduce Lagrange multipliers and quadratic penalty terms to transform constrained optimization into unconstrained optimization;

[0039] (43) The optimal values ​​are obtained by iteratively updating the modal function, center frequency and Lagrange multipliers using the alternating multiplier direction method;

[0040] In the formula, For the remaining components Wiener filtering results; The center frequency of the power spectrum of the modal function; Lagrange multipliers The frequency domain form of (t); Fourier transform of the input sequence f(t); For the k-th modal component u k Fourier transform of (t); n is the number of iterations; α is the bandwidth constraint penalty factor; τ is the update step size of the Lagrange multiplier; ω is the frequency variable; ω k The center frequency of the k-th modal component;

[0041] (44) Repeat steps (41)-(43) and stop iterating when the convergence condition is met.

[0042] Furthermore, constructing the TFLS hybrid forecast model includes:

[0043] The original time series data is synchronously input into the LSTM branch and the Transformer branch. The Transformer branch needs to be mapped through the embedding layer before receiving the input and then overlaid with position encoding to give the sequence position vector information.

[0044] The LSTM branch utilizes its internal forget gate, input gate, and output gate to perform time-step processing on the input sequence, extracting the local dynamic fluctuation features and short-term time dependencies within the sequence, and outputting a hidden state vector containing temporal information.

[0045] Global modeling is based on multi-head attention. The Transformer branch receives a feature vector with positional encoding, and then performs N... x Processed by N stacked Transformer units, x This represents the number of encoder units stacked in the Transformer branch; each Transformer unit uses a multi-head attention mechanism to calculate the correlation weights across the entire sequence.

[0046] Deeply fuse and output the model features, and perform parallel feature concatenation or weighted fusion of the hidden state vector output by LSTM and the global feature vector output by the last layer of the Transformer branch;

[0047] The fused composite feature vector is then subjected to nonlinear transformation through a feedforward network, followed by dimensionality mapping through a linear layer, and finally the prediction result is obtained through a Softmax regression layer.

[0048] Furthermore, in step S3, the input feature set is used as the input variable of the TFLS hybrid prediction model. After normalization, it is divided into training set, validation set and test set according to time order. During the training phase, the input feature set is used to learn the mapping relationship between each component sequence and the runoff target sequence. During the prediction phase, the input feature set corresponding to the time to be predicted is input into the trained TFLS hybrid prediction model.

[0049] Further, step S4 includes:

[0050] For each IMF component, its historical energy is calculated to determine the basic weights;

[0051] The corrected weights are obtained based on the correlation coefficient between each component and the concurrent rainfall sequence. The final fusion weights are obtained by combining the basic weights and the corrected weights.

[0052] The final hourly runoff forecast value is obtained based on the final fusion weight.

[0053] Furthermore, the formula for calculating historical energy is as follows:

[0054]

[0055] Where T is the length of the historical sequence, IMF i (t) represents the i-th EMD decomposition result.

[0056] Furthermore, the weighting formula is revised as follows:

[0057]

[0058] In the formula, r i Let be the correlation coefficient between the i-th IMF component and the rainfall sequence. Let be the correlation coefficient between the j-th IMF component and the rainfall sequence, and n be the total number of components.

[0059] Furthermore, the final fusion weights are calculated using a weighted combination, as shown in the following formula:

[0060]

[0061] In the formula, α is the fusion weight balance coefficient. Based on the weights, To adjust the weights.

[0062] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: (1) It specifically solves the technical bottlenecks of low accuracy and poor stability of existing hourly runoff forecasting methods when dealing with nonlinear and non-stationary runoff sequences; (2) It breaks through the limitations of traditional single decomposition methods and adopts a multi-step decomposition strategy of "CEEMDAN decomposition + high-frequency component VMD secondary decomposition", combined with fuzzy entropy complexity analysis and spectral clustering to accurately remove high-frequency interference and weaken the non-stationarity of the sequence; Before modeling, the PACF partial autocorrelation function is introduced to screen features, remove redundant information and strengthen the contribution of key factors, thereby improving the modeling efficiency. (3) The construction of the TFLS hybrid forecast model integrates the complementary advantages of Transformer and LSTM and realizes bidirectional information interaction, adapting to the differences in multi-frequency component features; for low-frequency trend components, LSTM accurately captures local dependencies and short-term features through gating mechanism, alleviating the gradient vanishing problem; for high-frequency fluctuation components, Transformer extracts global correlation features through multi-head attention mechanism, strengthening the ability to capture long-period dependencies; the two achieve dynamic fusion of local-global features through residual connection and linear mapping, improving the accuracy of peak characterization and baseflow simulation. Attached Figure Description

[0063] Figure 1 This is a flowchart of the method of the present invention;

[0064] Figure 2This is a structural diagram of the TFLS model of the present invention;

[0065] Figure 3 This is a simulation result diagram of the CEEMDAN-VMD-TFLS model during the testing period of this invention;

[0066] Figure 4 This is a simulation result diagram of the VMD-TFLS model during the testing period of this invention;

[0067] Figure 5 This is a simulation result diagram of the TFLS model during the testing period of this invention;

[0068] Figure 6 This is a simulation result diagram of the Transformer model during the testing phase;

[0069] Figure 7 This is a simulation result of the LSTM model during the testing period. Detailed Implementation

[0070] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0071] like Figure 1 As shown, the hourly runoff forecasting method based on CEEMD-TFLS of this invention includes:

[0072] S1. Collect hourly measured flow data in the lead-up to the watershed forecast section and hourly precipitation data from rain gauges within the watershed to construct complete hourly historical runoff and precipitation sequences. Perform missing value imputation, outlier handling, and normalization on the hourly historical runoff and precipitation sequences to obtain standardized measured hourly runoff and precipitation data.

[0073] S2. The standardized measured hourly runoff data are decomposed using the CEEMDAN method to obtain the IMF components and the k-th order residual components, and the fuzzy entropy value of each IMF component is calculated; the fuzzy entropy of each IMF component is calculated, and the average similarity is solved. To reduce the complexity, the fuzzy entropy characteristics of each IMF component are used to reconstruct it into high, medium and low frequency sequences using spectral clustering. The high frequency sequence is then subjected to VMD secondary decomposition to further stabilize the high frequency signal.

[0074] The adaptive decomposition of standardized measured hourly runoff data using the Complete Ensemble Empirical Mode Decomposition (CEEMDAN) method includes:

[0075] (1) Gaussian white noise was added to the standardized measured hourly runoff data to form the overall decomposed signal:

[0076]

[0077] In the formula, The noise figure represents the control over the added noise and the signal-to-noise ratio of the original signal. (i=1,2,…,I) represents a set of Gaussian white noise with a mean of 0 and a variance of 1, where I represents the ensemble order, i.e., achieving I-fold EMD decomposition. To standardize measured hourly runoff data, The sequence after adding Gaussian white noise to the standardized measured hourly runoff data.

[0078] (2) The overall average value is calculated by repeating the EMD decomposition process I times, and it is defined as the first IMF component of the original signal. ,Right now:

[0079]

[0080] In the formula, (i=1,2,…,I) represents the i-th EMD decomposition result, and I represents the number of times the EMD decomposition is repeated.

[0081] (3) Calculate the first-order residual :

[0082]

[0083] In the formula, This is the first IMF component.

[0084] (4) After adding specific noise to the first-order residual, perform EMD decomposition again and calculate the average value to obtain the second IMF component. ,Right now:

[0085]

[0086] In the formula, The noise figure represents the noise figure used to control the added noise and the signal-to-noise ratio of the first-order residual. (i=1,2,⋯,I) represents a set of Gaussian white noise with a mean of 0 and a variance of 1, where I represents the ensemble order, i.e., achieving I-fold EMD decomposition, and E k (·) represents the operators that have been defined during the process of obtaining the components.

[0087] (5) Calculate the k-th order residual :

[0088]

[0089] In the formula, For the (k-1)th order residual, IMF k This is the k-th IMF component.

[0090] (6) After adding specific noise to the k-th residual and performing EMD decomposition again, the average value is calculated to obtain the (k+1)-th IMF component. （k+1） (t), that is:

[0091]

[0092] In the formula, E represents the noise figure, which controls the added noise and the k-th order residual noise ratio. k (·) represents the operators that have been defined during the process of obtaining the components.

[0093] (7) Repeat steps (5)-(6) until the residuals cannot be decomposed, thus obtaining the final residuals. :

[0094]

[0095] (8) Standardized measured hourly runoff data The expression is:

[0096]

[0097] In the formula, K represents the total number of IMF components decomposed from the original sequence at a given precision.

[0098] Specifically, the fuzzy entropy of each IMF component is calculated, and the average similarity is solved. The complexity of quantification includes:

[0099] (1) Take the phase space dimension m, and construct N-m+1 m-dimensional vectors X. m，i ,Right now:

[0100]

[0101]

[0102] In the formula, N represents the length of the IMF component time series for which the fuzzy entropy is to be calculated, that is, the total number of data points contained in the IMF component; For vector [x i x i+1 , ..., x i+1-m The mean of ].

[0103] (2) Define vector X m，i X m，j absolute distance d ij m The maximum value of the difference between its corresponding elements:

[0104]

[0105] In the formula, i, j = 1, 2, ..., Nm, i ≠ j.

[0106] (3) Introduce the fuzzy function μ(d ij m Define X(n, r) m，i X m，j Similarity The expression is:

[0107]

[0108] In the formula, r is the tolerance width and n is the tolerance gradient.

[0109] (4) Define the function :

[0110]

[0111] This function represents the average fuzzy similarity between all pairs of embedded vectors when the phase space dimension is m, and is used to characterize the overall similarity or regularity of the IMF component in the m-dimensional phase space; where r is the tolerance width and n is the boundary gradient parameter of the fuzzy function.

[0112] (5) Change the phase space dimension to m+1, and repeat steps (1)-(4) to obtain the function. :

[0113]

[0114] This function represents the average fuzzy similarity between all pairs of embedded vectors when the phase space dimension increases to m+1.

[0115] (6) When the sample has N data points, the fuzzy entropy is defined as:

[0116]

[0117] VMD-based decomposition of high-frequency sequences includes:

[0118] (1) By using the Intrinsic Mode Function (IMF) k Converted into a bandwidth-limited amplitude-frequency modulation signal μ k The minimum sum of bandwidths estimated by intrinsic mode functions is used as the constraint condition for decomposition:

[0119]

[0120] And it satisfies:

[0121]

[0122] In the formula, k is the number of modal components obtained by decomposition; ω is the k-th modal component; k Let be the center angular frequency of the k-th modal component; K is the total number of modal components; It is the Dirac function; This represents the partial derivative operation with respect to the time variable t; This is a complex exponential modulation term.

[0123] (2) Introducing Lagrange multipliers The quadratic penalty term α transforms constrained optimization into unconstrained optimization.

[0124] (3) The alternating multiplier direction method is used to analyze the mode function μ. k Center frequency and Lagrange multipliers Iterative updates are performed to obtain the optimal value.

[0125]

[0126]

[0127]

[0128] In the formula, For the remaining components Wiener filtering results; The center frequency of the power spectrum of the modal function; Lagrange multipliers The frequency domain form of (t); Fourier transform of the input sequence f(t); For the k-th modal component u k Fourier transform of (t); n is the number of iterations; α is the bandwidth constraint penalty factor; τ is the update step size of the Lagrange multiplier; ω is the frequency variable; ω k Let be the center frequency of the k-th modal component.

[0129] (4) Repeat steps (1)-(3) until the condition is met. If convergence is considered, the iteration stops, where f is the threshold.

[0130] S3. Use the partial autocorrelation function (PACF) combined with feature importance analysis to screen key input factors, and align them according to the relationship between forecast time and lag time to form an input feature set; construct a TFLS hybrid forecast model, input the features to be forecasted to obtain runoff prediction results for each component or time scale;

[0131] The high-frequency subsequence, mid-frequency sequence, and low-frequency sequence are respectively input into the TFLS model, which is composed of LSTM and Transformer coupling. See the TFLS model diagram. Figure 2The high-frequency sub-model TFLS was obtained through training. H TFLS (Intermediate Frequency Submodel) M and low-frequency sub-model TFLS L This enables segmented forecasting of runoff components in different frequency bands.

[0132] The input feature set, as the input variable of the TFLS hybrid prediction model, is normalized and then divided into training set, validation set and test set in a 3:1:1 ratio according to time order.

[0133] The key input factors refer to the candidate feature set composed of standardized measured hourly runoff and precipitation data preprocessed in step S1 and the modal components at various scales generated in step S2. This set includes the historical dependence of runoff itself, precipitation driving factors, and modal decomposition components. Using the runoff and precipitation sequences from step S1 and the modal components obtained from step S2 as candidate factors, PACF calculations are used to eliminate intermediate lag interference, quantify the direct linear correlation between each candidate factor and the target runoff sequence, and screen out key factors highly correlated with runoff forecasting.

[0134] The TFLS model is trained using the training set data, and key hyperparameters such as regularization parameters, number of iterations, number of attention heads, and LSTM hidden layer dimensions are optimized and adjusted to dynamically optimize model performance. An early stopping mechanism is introduced during training to avoid overfitting and ensure the model's generalization ability on the data.

[0135] During the training phase, the mapping relationship between each component sequence and the runoff target sequence is learned using the input feature set;

[0136] In the prediction phase, the input feature set corresponding to the time to be predicted is input into the trained TFLS hybrid prediction model to obtain the runoff prediction results for each component or time scale.

[0137] The PACF is calculated using the following formula:

[0138]

[0139] In the formula, k is the lag order (k = 1, 2, ..., K, where K is the maximum lag order); x t x represents the value of the target runoff sequence at time t; t-k Let be the value of the candidate factor sequence at time tk; For x based on the first k-1 lags t The linear predicted value; For x based on the first k-1 lags t-k The linear predicted value.

[0140] The prediction model TFLS integrates LSTM and Transformer deep learning models through parallel heterogeneous feature fusion, specifically including:

[0141] (1) The original time series data is synchronously input into the LSTM branch and the Transformer branch. The Transformer branch needs to be mapped through the embedding layer before receiving the input and then overlaid with position encoding to give the sequence position vector information.

[0142] (2) The LSTM branch uses its internal forget gate, input gate and output gate to process the input sequence step by step, extract the local dynamic fluctuation features and short-term time dependencies within the sequence, and output the hidden state vector h containing time series information. t ;

[0143] (3) Global modeling is performed based on multi-head attention. The Transformer branch receives a feature vector with positional encoding and then performs N... x Processed by N stacked Transformer units (N x This represents the number of stacked encoder units in the Transformer branch, used to characterize the network depth of the Transformer branch. Each unit utilizes a multi-head attention mechanism to calculate association weights across the entire sequence, achieving global feature extraction and long-distance dependency modeling, and employs residual connections and layer normalization to ensure the stability of feature propagation.

[0144] (4) Perform deep fusion and output of model features, and output the hidden state vector h output by LSTM. t Parallel feature concatenation or weighted fusion is performed with the global feature vector output from the last layer of the Transformer branch.

[0145] Assume the feature vector output by LSTM is The feature vector output by Transformer is Its coupling and output process can be represented as:

[0146]

[0147] in, This indicates a feature splicing operation.

[0148] (5) The fused model feature vectors are sequentially transformed nonlinearly through a feedforward network, then dimension-mapped through a linear layer, and finally passed through a softmax regression layer to calculate the final prediction result.

[0149]

[0150] Wherein, FFN stands for feedforward neural network mapping. This is the output layer weight matrix. This is the output layer bias vector. ŷ represents the feature vector resulting from the fusion of the LSTM and Transformer branches, and ŷ represents the final runoff prediction value.

[0151] Standardize and dimensionally integrate the features selected in step S3 to construct a multi-dimensional input matrix; divide the training / validation / test sets into a 3:1:1 ratio. Train the TFLS model parameters and optimize the hyperparameters; introduce an early stopping mechanism to prevent overfitting. Input the feature variables to be predicted to obtain the hourly-scale prediction results for each component;

[0152] The standardized measured hourly runoff is decomposed using CEEMD and VMD, and combined with the TFLS model to obtain the CEEMD-VMD-TFLS model for hourly runoff prediction.

[0153] Through the synergistic effect of secondary decomposition, precise screening and hybrid modeling, the error of the final reconstructed forecast result can be effectively reduced, and it is more adaptable to extreme scenarios, meeting the high accuracy and high stability requirements of hourly-scale forecasts for practical engineering projects such as reservoir scheduling.

[0154] S4. The basic weights are determined based on the historical energy proportion of each modal component, and the weights are adjusted by combining the rainfall-runoff correlation strength. The runoff prediction results of each sub-model are then weighted and superimposed to obtain the final hourly runoff forecast value. The accuracy and stability under different hydrological scenarios such as rainstorms and dry seasons are quantitatively evaluated using RMSE and NSE indices.

[0155] Among them, the historical energy proportion of each modal component objectively reflects its contribution to the original runoff. The higher the energy of the component, the greater its impact on the runoff process. Determining the basic weight in this way is in line with physical meaning. Then, the weight is corrected by the rainfall-runoff correlation, which can strengthen the constraint of hydrophysical mechanism on the component weight, so that the weight takes into account both historical energy contribution and runoff generation and confluence patterns, thereby improving the physical rationality and forecast accuracy of the model.

[0156] The specific method for determining the weights is as follows:

[0157] (1) Calculate the energy percentage of each component:

[0158] For each IMF component Calculate its historical energy :

[0159]

[0160] Basic weights for:

[0161]

[0162] Where T is the length of the historical sequence, n is the total number of components, and the weight sum is 1. Let the historical energy of the j-th IMF component be represented;

[0163] (2) Correction based on rainfall-runoff correlation:

[0164] Calculate IMF for each component i The correlation coefficient r with the concurrent rainfall series P(t) i Correcting weight w i,corr for:

[0165]

[0166] Final fusion weight w i Weighted combinations can be used:

[0167]

[0168] in, Let α be the correlation coefficient between the j-th IMF component and the rainfall sequence P(t), and let α be the fusion weight balancing coefficient used to adjust the base weights. Corrected weight w for precipitation-runoff correlation i,corr The proportion of α in the final weight. The candidate value range of α can be set to [0.5, 0.7]. This range allows the historical energy contribution to play a dominant role in the final weight, while retaining the ability of rainfall-runoff correlation to correct the weight. The specific value is determined by optimization on the validation set, that is, α that minimizes the prediction error of the validation set is selected as the final balance coefficient within the candidate range.

[0169] (3) Weighted superposition reconstruction:

[0170] Final hourly runoff forecast for:

[0171]

[0172] Example:

[0173] This embodiment uses the Tingzikou Station in the Jialing River Basin as an example to illustrate how the method of this invention is used to forecast the hourly flow rate at the Tingzikou Station.

[0174] The Tingzikou Hydropower Project, a key control project on the main stream of the Jialing River, has a total reservoir capacity of 4.067 billion m³ and a maximum flood control capacity of 1.44 billion m³, playing a crucial role in the basin's flood control system.

[0175] Currently, hourly runoff data and rainfall data from multiple stations in the upstream basin of the Tingzikou Hydropower Project from 2014 to 2024 are being collected.

[0176] First, using hourly runoff from Tingzikou station from 2014 to 2024 and hourly rainfall from multiple upstream stations as data samples, missing values ​​were filled, outliers were handled, and normalization was performed on the runoff and rainfall sequences.

[0177] The preprocessed standardized hourly runoff data is decomposed using the Complete Ensemble Empirical Mode Decomposition (CEEMDAN) method, which decomposes the runoff sequence at the forecast section (Tingzikou Station) into a series of intrinsic mode function (IMF) components and one residual component.

[0178] For each IMF component obtained from the decomposition, its fuzzy entropy (FE) value is calculated. Based on the characteristics of each IMF component and its fuzzy entropy, spectral clustering (SC) is used to reconstruct all IMF components into three classes of sequences: high-frequency, medium-frequency, and low-frequency. Variational mode decomposition (VMD) is used to perform a secondary decomposition on the high-frequency sequences to further reduce the instability of the IMF signal. Partial autocorrelation coefficient (PACF) combined with feature importance analysis is used to accurately screen potential input factors, and a Transformer-LSTM fusion (TFLS) hybrid forecasting model is constructed. The screened feature factors are input into the TFLS model, and the TFLS model is trained and optimized. The hourly runoff prediction results during the test period are as follows:

[0179] To verify the reliability and practicality of the CEEMD-TFLS model in this invention, the runoff and precipitation series from 2014 to 2024 were decomposed into three parts: the six-year runoff and precipitation series from 2014 to 2020 were used as training data, the runoff and precipitation series from 2021 to 2022 were used as validation data, and the runoff and precipitation series from 2023 to 2024 were used as test data. CEEMD and VMD were used to perform a secondary decomposition of the runoff series, and the hourly runoff was predicted using the hybrid deep learning model TFLS. The prediction results are as follows: Figure 3 Simultaneously, the TFLS model prediction using only VMD decomposition of the sequence is calculated, and the prediction results are as follows: Figure 4 The results also include hourly runoff predictions using single-model TFLS, Transformer, and LSTM models, as shown in the figures below. Figure 5 , Figure 6 and Figure 7 The simulation performance of the model was evaluated using the commonly used Nash coefficient (NSE) and minimum root mean square error (RMSE) indices. The results are shown in Table 1.

[0180] Table 1 Experimental Results

[0181]

[0182] Table 1 shows that the CEEMD-VMD-TFLS model demonstrates significant superiority over other models in both 24-hour forecasts and the average forecast for the previous 24 hours. Furthermore, the VMD-TFLS model, after VMD decomposition, also outperforms the single model, indicating that the secondary decomposition technique combining CEEMD and VMD performs better in hourly runoff sequence simulation and prediction compared to conventional single decomposition methods. In addition, the TFLS model based on a hybrid deep learning framework also outperforms individual prediction models such as LSTM and Transformer, demonstrating that the complementary advantages of Transformer and LSTM can effectively improve the accuracy of the prediction model.

[0183] This invention is applicable not only to hourly runoff forecasting but also to daily, monthly, and other multi-timescale runoff forecasting. This method provides a general technical reference for processing nonlinear and non-stationary time series, offering reliable technical support for basin flood early warning and refined water resource management. It also provides valuable insights for innovation in time series modeling in related fields, possessing significant theoretical importance, practical application value, and broad prospects for promotion. Furthermore, it provides new technical means and methodological support for reservoir scheduling and basin flood control and disaster reduction.

Claims

1. A method for hourly runoff forecasting based on CEEMD-TFLS, characterized in that, include: S1. Collect hourly measured flow data from the watershed forecast section and upstream hydrological stations, as well as hourly precipitation data from rain gauge stations within the watershed, to construct a complete hourly historical runoff and precipitation sequence; perform missing value imputation, outlier handling, and normalization on the hourly historical runoff and precipitation sequences to obtain standardized measured hourly runoff and precipitation data. S2. The standardized measured hourly runoff data are decomposed using complete ensemble empirical mode decomposition to obtain IMF components and k-th order residual components; the fuzzy entropy of each IMF component is calculated, and based on the fuzzy entropy characteristics of each IMF component, it is reconstructed into high, medium and low frequency sequences using spectral clustering method, and the high frequency sequence is subjected to VMD secondary decomposition. S3. Use partial autocorrelation function combined with feature importance analysis to screen key input factors and align them with time series to form an input feature set; construct a TFLS hybrid forecasting model, input the features to be forecasted to obtain runoff prediction results; S4. Determine the basic weights based on the historical energy proportions of each modal component, and adjust the weights by combining the correlation strength between rainfall and runoff. Then, overlay and reconstruct the runoff prediction results to obtain the hourly runoff forecast values.

2. The hourly runoff forecasting method according to claim 1, characterized in that, The standardized measured hourly runoff data were decomposed using complete ensemble empirical mode decomposition, including: (21) Gaussian white noise was added to the standardized measured hourly runoff data to form the overall decomposed signal; (22) The overall average value is calculated by repeating the EMD decomposition process I times, and it is defined as the first IMF component of the original signal, where I is the number of times the EMD decomposition is repeated; (23) Calculate the first-order residual based on the standardized measured hourly runoff data with added Gaussian white noise and the first IMF component of the original signal; (24) After adding specific noise to the first-order residual, perform EMD decomposition again and calculate the average value to obtain the second IMF component; (25) Calculate the k-th order residual; (26) After adding specific noise to the k-th residual and performing EMD decomposition again, the average value is obtained to obtain the (k+1)th IMF component; (27) Repeat steps (25)-(26) until the residuals cannot be decomposed, and obtain the final residuals; (28) The standardized measured hourly runoff data after decomposition are obtained based on the (k+1)th IMF component and the final residual.

3. The hourly runoff forecasting method according to claim 1, characterized in that, The calculation of the fuzzy entropy of each IMF component includes: (31) Take the phase space dimension m, and construct N-m+1 m-dimensional vectors X. m，i N represents the length of the IMF component time series of the fuzzy entropy to be calculated; (32) Define vector X m，i X m，j absolute distance It is the maximum value of the difference between its corresponding elements; (33) Introduce the fuzzy function μ(d ij m Define X(n, r) m，i X m，j Similarity The expression is: In the formula, r is the tolerance width and n is the tolerance gradient; (34) Define the function : (35) Change the phase space dimension to m+1, repeat steps (31)-(34), and obtain the function. : (36) When the sample has N data points, the fuzzy entropy is defined as:

4. The hourly runoff forecasting method according to claim 1, characterized in that, The VMD secondary decomposition of high-frequency sequences includes: (41) By using the intrinsic mode function (IMF) k Converted into a bandwidth-limited amplitude-frequency modulation signal μ k The minimum sum of bandwidths estimated by intrinsic mode functions is used as the constraint condition for decomposition: And it satisfies: In the formula, To standardize the measured hourly runoff data, k is the number of modal components obtained from the decomposition; ω is the k-th modal component; k Let be the center angular frequency of the k-th modal component; K is the total number of modal components; It is the Dirac function; This represents the partial derivative operation with respect to the time variable t; It is a complex exponential modulation term; (42) Introduce Lagrange multipliers and quadratic penalty terms to transform constrained optimization into unconstrained optimization; (43) The optimal values ​​are obtained by iteratively updating the modal function, center frequency and Lagrange multipliers using the alternating multiplier direction method; In the formula, For the remaining components Wiener filtering results; The center frequency of the power spectrum of the modal function; Lagrange multipliers The frequency domain form of (t); Fourier transform of the input sequence f(t); For the k-th modal component u k Fourier transform of (t); n is the number of iterations; α is the bandwidth constraint penalty factor; τ is the update step size of the Lagrange multiplier; ω is the frequency variable; ω k The center frequency of the k-th modal component; (44) Repeat steps (41)-(43) and stop iterating when the convergence condition is met.

5. The hourly runoff forecasting method according to claim 1, characterized in that, Constructing a TFLS hybrid forecasting model includes: The original time series data is synchronously input into the LSTM branch and the Transformer branch. The Transformer branch needs to be mapped through the embedding layer before receiving the input and then overlaid with position encoding to give the sequence position vector information. The LSTM branch utilizes its internal forget gate, input gate, and output gate to perform time-step processing on the input sequence, extracting the local dynamic fluctuation features and short-term time dependencies within the sequence, and outputting a hidden state vector containing temporal information. Global modeling is based on multi-head attention. The Transformer branch receives a feature vector with positional encoding, and then performs N... x Processed by N stacked Transformer units, x This represents the number of encoder units stacked in the Transformer branch; each Transformer unit uses a multi-head attention mechanism to calculate the correlation weights across the entire sequence. Deeply fuse and output the model features, and perform parallel feature concatenation or weighted fusion of the hidden state vector output by LSTM and the global feature vector output by the last layer of the Transformer branch; The fused composite feature vector is then subjected to nonlinear transformation through a feedforward network, followed by dimensionality mapping through a linear layer, and finally the prediction result is obtained through a Softmax regression layer.

6. The hourly runoff forecasting method according to claim 1, characterized in that, In step S3, the input feature set is used as the input variable of the TFLS hybrid prediction model. After normalization, it is divided into training set, validation set and test set according to time order. During the training phase, the input feature set is used to learn the mapping relationship between each component sequence and the runoff target sequence. During the prediction phase, the input feature set corresponding to the time to be predicted is input into the trained TFLS hybrid prediction model.

7. The hourly runoff forecasting method according to claim 1, characterized in that, Step S4 includes: For each IMF component, its historical energy is calculated to determine the basic weights; The corrected weights are obtained based on the correlation coefficient between each component and the concurrent rainfall sequence. The final fusion weights are obtained by combining the basic weights and the corrected weights. The final hourly runoff forecast value is obtained based on the final fusion weight.

8. The hourly runoff forecasting method according to claim 7, characterized in that, The formula for calculating historical energy is as follows: Where T is the length of the historical sequence, IMF i (t) represents the i-th EMD decomposition result.

9. The hourly runoff forecasting method according to claim 7, characterized in that, The corrected weight formula is as follows: In the formula, r i Let be the correlation coefficient between the i-th IMF component and the rainfall sequence. Let be the correlation coefficient between the j-th IMF component and the rainfall sequence, and n be the total number of components.

10. The hourly runoff forecasting method according to claim 7, characterized in that, The final fusion weights are calculated using a weighted combination, as shown in the following formula: In the formula, α is the fusion weight balance coefficient. Based on the weights, To adjust the weights.

Images