Multi-mode daily runoff prediction method and system based on VMD-CEEDMAN

Abstract

The present invention relates to the technical field of water runoff prediction, and specifically discloses a multi-mode daily runoff prediction method and system based on VMD-CEEDMAN. Aiming at the periodicity and randomness of daily runoff, VMD and CEEMDAN are firstly used to decompose historical runoff and water level data respectively, extract information of different frequency bands, reduce the complexity of original data sequence, give play to the advantages of different decomposition algorithms, and more completely highlight the fluctuation trend; then the decomposed daily runoff data are divided into stationary sequence and non-stationary sequence, and the stationary sequence includes medium-low frequency sequence and high frequency sequence; different prediction models are adopted for different sequences; finally, the prediction results of different subsequences are reconstructed to obtain the daily runoff prediction result.

Description

Technical Field

[0001] The invention relates to the technical field of water runoff prediction, and in particular to a multi-mode daily runoff prediction method and system based on VMD-CEEDMAN. Background Art

[0002] Runoff is an important indicator reflecting the hydrological process of a river basin. The simulation and prediction of the changing trend of runoff in a river basin is an important topic in the field of hydrological research. Accurate runoff forecast is also an important prerequisite for ensuring water security in the river basin, preventing water and drought disasters, and improving the comprehensive development and utilization of water resources. Currently, the methods for runoff prediction are mainly divided into mechanism-based and non-mechanism-based.

[0003] Mechanism-based prediction methods are often based on some distributed hydrological models, such as the SWAT model and the DHSVM model. Although these models can simulate some basin hydrological processes very well, the model construction has high requirements for basin underlying surface and hydrological data and the parameter setting process is complex. There are problems such as parameter uncertainty and poor model universality, and it is difficult to apply to areas where basic data are lacking. The process of runoff formation is not a simple physical process, but also includes chemical and biological processes, and is disturbed by human activities. Mechanism-based physical models are difficult to fully consider these factors and often cannot obtain good prediction results. In particular, in recent years, many water conservancy projects have been built in various places, resulting in significant changes in the evolution of river confluence; superimposed on the significant changes in the meteorological and hydrological characteristics of the basin caused by global climate change, the continuity and adaptability of the traditional basin runoff generation and confluence model have been destroyed, greatly increasing the uncertainty of the prediction results.

[0004] Non-mechanistic methods generalize complex physical processes into black box models to improve the model migration ability. However, due to the nonlinearity, mutation, and noise characteristics of hydrological systems, the prediction effects of models based on statistics and traditional machine learning are often unsatisfactory. With the development of artificial neural networks and the improvement of deep learning technology, methods such as convolutional neural networks (CNNs) and recurrent neural networks (RNNs) have shown good nonlinear mapping capabilities, powerful feature extraction technology, and self-learning generalization capabilities. They can adapt well to complex time series data such as runoff, and are therefore widely used in the research of hydrological forecasting and environmental monitoring. However, although the CNN model can perform convolution operations on each time series and extract local features from runoff information well, CNN is not sensitive to the time order of runoff information and cannot complete the prediction task alone; and ordinary RNNs are prone to gradient expansion and gradient disappearance problems.

[0005] Long and Short-term Memory (LSTM) is a special structure recurrent neural network designed to alleviate the long-term dependency problem and gradient vanishing problem of RNN network. It has obvious advantages in processing data with long time steps and nonlinear sequences, but due to its recursive nature, it may be lacking in processing local features. In addition, a single LSTM network is difficult to capture the various complex patterns of data and is more sensitive to noise, especially when the input data quality is poor, which may affect its prediction performance. Research results have shown that a hybrid network model composed of multiple neural networks can help improve prediction accuracy.

[0006] In summary, the existing runoff hydrological forecasting methods have the following problems: First, affected by factors such as the dynamics of the basin's underlying surface, climate change, and human interference, the non-stationarity of the hydrological time series has increased, and traditional mechanism-based prediction methods are difficult to solve problems such as large hydrological forecast deviations and poor model universality. Second, compared with traditional methods, artificial intelligence algorithms have many advantages, but different algorithms have different characteristics and scopes of application, and a single algorithm is difficult to meet the accuracy requirements of hydrological forecasts. Summary of the invention

[0007] The purpose of the present invention is to provide a multi-mode daily runoff prediction method and system based on VMD-CEEDMAN to solve the above technical problems.

[0008] The purpose of the present invention can be achieved through the following technical solutions:

[0009] The multi-mode daily runoff prediction method based on VMD-CEEDMAN includes the following steps:

[0010] Step S1: collecting historical data, the historical data including daily runoff data and water level data; preprocessing and time sorting the historical data to obtain complete historical time series data;

[0011] Step S2: performing VMD decomposition on the daily runoff data to obtain a number of submodal components of different frequencies; performing stationarity determination on the submodal components according to a preset significance level threshold, dividing the submodal components into stationary submodal components and non-stationary submodal components, and dividing the stationary submodal components into medium-low frequency components and high frequency components;

[0012] Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components;

[0013] Step S3: Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component and obtain the average importance score Imp of each feature ave , the average importance score is the average value of the impurity reduction; set the impurity threshold I, select Imp ave Features <I are recorded as key features;

[0014] Step S4: obtaining a non-stationary input sequence according to the key features and their corresponding historical data; obtaining a medium-low frequency input sequence according to the daily runoff data corresponding to the medium-low frequency components, and obtaining a high frequency input sequence according to the daily runoff data corresponding to the high frequency components;

[0015] Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix;

[0016] Step S5: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction;

[0017] The prediction results of multi-feature prediction and single-feature prediction are obtained and recorded as sub-modal component values; the sub-modal component values ​​are accumulated and reconstructed to obtain the daily runoff prediction result.

[0018] As a further solution of the present invention: in step S1, the process of preprocessing the historical data includes:

[0019] Abnormal values ​​in the historical data are obtained, and the abnormal values ​​are blanked to obtain missing values, and the missing values ​​are filled by using a linear interpolation method.

[0020] As a further solution of the present invention: in step S2, the process of determining the stability of the submodal component includes:

[0021] The significance level threshold Ht=0.05;

[0022] The stability determination is based on the ADF test method, assuming that the submodal component has a unit root, and the statistics Sta and P value of the submodal component are obtained through the ADF test;

[0023] If the P value of the submodal component is ≤Ht, then the submodal component does not have a unit root under the preset significance level threshold, that is, it is a stationary submodal component; if the P value of the submodal component is greater than Ht, then the submodal component has a unit root under the preset significance level threshold, that is, it is a non-stationary submodal component.

[0024] As a further solution of the present invention: in step S2, the process of dividing the stationary submodal components includes:

[0025] The division of the stationary submodal component is based on the zero-crossing method, and the stationary submodal component is recorded as a detection signal; the frequency of the detection signal when it crosses the zero point is obtained, which is recorded as the zero-crossing rate Zcr, and the zero point is the point where the signal value changes from positive to negative or from negative to positive;

[0026] If the zero-crossing rate Zcr of the detection signal is greater than 0.2, the stationary submodal component corresponding to the detection signal is a high-frequency component; if the zero-crossing rate Zcr is less than or equal to 0.2, the stationary submodal component corresponding to the detection signal is a medium-low frequency component.

[0027] As a further solution of the present invention: in step S2, the process of obtaining the trend component and the detail component includes:

[0028] The modal component is recorded as a new detection signal, the zero-crossing rate Zdr of the new detection signal is obtained, and the zero-crossing rate Zdr of the new detection signal is determined; if the new detection signal is a high-frequency component, it is recorded as a detail component; if the new detection signal is a medium-to-low-frequency component, it is recorded as a trend component.

[0029] As a further solution of the present invention: in step S4, the process of the normalization operation includes:

[0030] Set the normalization range threshold [-1, 1]; record the non-stationary input sequence, the medium-low frequency input sequence and the high frequency input sequence as the matrix to be normalized array; then the process of normalizing the matrix to be normalized is:

[0031]

[0032] Where array′ is the value after normalization, min(array) is the minimum value in the matrix to be normalized, and max(array) is the maximum value in the matrix to be normalized.

[0033] As a further solution of the present invention: a multi-mode daily runoff prediction system based on VMD-CEEDMAN, comprising:

[0034] Data collection module: collects historical data, including daily runoff data and water level data; pre-processes and time-sorts the historical data to obtain complete historical time series data;

[0035] Training data set module: Perform VMD decomposition on the daily runoff data to obtain several submodal components of different frequencies; perform stationarity determination on the submodal components according to a preset significance level threshold, divide the submodal components into stationary submodal components and non-stationary submodal components, and divide the stationary submodal components into medium-low frequency components and high frequency components;

[0036] Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components;

[0037] Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component and obtain the average importance score Imp of each feature ave , the average importance score is the average value of the impurity reduction; set the impurity threshold I, select Imp ave Features <I are recorded as key features;

[0038] According to the key features and their corresponding historical data, a non-stationary input sequence is obtained; according to the daily runoff data corresponding to the medium and low frequency components, a medium and low frequency input sequence is obtained, and according to the daily runoff data corresponding to the high frequency components, a high frequency input sequence is obtained;

[0039] Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix;

[0040] Model building module: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction;

[0041] Prediction module: obtain the prediction results of multi-feature prediction and single-feature prediction, recorded as sub-modal component values; accumulate and reconstruct the values ​​of each sub-modal component to obtain the daily runoff prediction result.

[0042] Beneficial effects of the present invention:

[0043] The existing runoff hydrological forecasting methods have the following problems: First, affected by factors such as the dynamics of the basin's underlying surface, climate change and human interference, the non-stationarity of the hydrological time series has increased, and the traditional mechanism-based prediction methods are difficult to solve the problems of large hydrological forecast deviations and poor model universality. Second, compared with traditional methods, artificial intelligence algorithms have many advantages, but different algorithms have different characteristics and scopes of application, and a single algorithm is difficult to meet the accuracy requirements of hydrological forecasts; in view of the above situation, in order to overcome the defects of the prior art, the present invention provides a method for runoff prediction by combining CEEMDAN and VMD signal decomposition methods for feature extraction, and then coupling the combined network model of Attention+BiGRU multi-feature prediction, BiLSTM single-feature prediction and TCN+BiLSTM single-feature prediction algorithms, which has a higher prediction accuracy and can track the strong trend characteristics and strong periodic characteristics of runoff data.

[0044] Compared with the prior art, the present invention aims at the periodicity and randomness of daily runoff. Firstly, VMD and CEEMDAN are used to decompose the historical runoff and water level data respectively, extract information of different frequency bands, reduce the complexity of the original data sequence, give play to the advantages of different decomposition algorithms, and highlight the fluctuation trend more completely; then the decomposed daily runoff data is divided into stationary sequence and non-stationary sequence, and different prediction models are adopted for different sequences; finally, the prediction results of different subsequences are reconstructed to obtain the daily runoff prediction results.

[0045] In the simulation analysis, it was found that LSTM is prone to gradient explosion and local optimal problems due to the influence of deep structure. The present invention first uses variational mode decomposition (PSO-VMD) optimized by particle swarm algorithm to decompose the runoff, and the obtained sub-signal tends to be more stable, which is more conducive to feature extraction of the model; then, according to the stability, frequency range and other characteristics of the sub-modal components, a bidirectional long short-term memory network (BiLSTM) and a time convolution network-long short-term memory neural network (TCN-LSTM) are used to predict different sub-modal components, which can overcome the shortcomings of a single LSTM. In terms of computational efficiency, since different decomposition methods are used for different signals and different models are used for prediction of different sub-modes, the extraction of key features is effectively accelerated, and the computational efficiency of the model is further improved. The present invention has obvious advantages in the accuracy of predicted values, the efficiency of the model and the computational efficiency, and can provide more timely and accurate hydrological information for comprehensive water resources management and flood and drought disaster prevention. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] The present invention will be further described below in conjunction with the accompanying drawings.

[0047] Figure 1It is a flow chart of the multi-mode daily runoff prediction method based on VMD-CEEDMAN of the present invention;

[0048] Figure 2 It is a schematic diagram of the LSTM network structure in the multi-mode daily runoff prediction method and system based on VMD-CEEDMAN of the present invention;

[0049] Figure 3 It is a schematic diagram of the BiLSTM network structure in the multi-mode daily runoff prediction method based on VMD-CEEDMAN of the present invention and the system;

[0050] Figure 4 It is a schematic diagram of the division results in the multi-mode daily runoff prediction method and system based on VMD-CEEDMAN of the present invention;

[0051] Figure 5 It is a schematic diagram comparing runoff prediction results between different models in the multi-mode daily runoff prediction method based on VMD-CEEDMAN of the present invention and the system. DETAILED DESCRIPTION

[0052] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0053] See also Figure 1 As shown, the present invention is a multi-mode daily runoff prediction method based on VMD-CEEDMAN, comprising the following steps:

[0054] Step S1: collecting historical data, the historical data including daily runoff data and water level data; preprocessing and time sorting the historical data to obtain complete historical time series data;

[0055] As a preferred embodiment of the present invention, the process of preprocessing the historical data includes:

[0056] Obtaining abnormal values ​​in the historical data, and blanking the abnormal values ​​to obtain missing values, and filling the missing values ​​by using a linear interpolation method;

[0057] Step S2: performing VMD decomposition on the daily runoff data to obtain a number of submodal components of different frequencies; performing stationarity determination on the submodal components according to a preset significance level threshold, dividing the submodal components into stationary submodal components and non-stationary submodal components, and dividing the stationary submodal components into medium-low frequency components and high frequency components;

[0058] Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components;

[0059] It should be noted that the process of decomposing the daily runoff data is VMD decomposition, and the specific process includes:

[0060] The particle swarm optimization algorithm is used to optimize the parameters of VMD decomposition to obtain several submodal components with different frequencies; the ADF test method is used to determine the stability of the obtained submodal components, and each submodal component is divided into a stable submodal component and a non-stationary submodal component; and the stable submodal component is divided into a medium-low frequency component and a high frequency component;

[0061] It should be noted that the VMD decomposition into Variational Modal Decomposition (VMD) is a new time-frequency analysis method that decomposes a multi-component signal into several single-component amplitude-modulated frequency-modulated signals at one time; the ADF test (Augmented Dickey-Fuller test) is a method for testing whether time series data has a unit root (Unit Root), which is often used to detect the stationarity of data;

[0062] The process of performing VMD decomposition on the daily runoff data includes:

[0063] The VMD decomposition is transformed into solving the following constraints:

[0064]

[0065] Wherein, f(t) is the original signal generated according to the daily runoff data, {u K} represents the decomposed K modal functions, u K is the modal function, {ω K} represents the center frequency of each mode, ω K is the modal center frequency, * represents the convolution operation, δ(t) represents the Dirac function, t is the time variable, is the symbol of partial derivative, j is the imaginary unit, |||| is the L2 norm;

[0066] And the above constraints are transformed into an unconstrained optimization problem by augmenting the Lagrangian function, and the following equation is obtained:

[0067]

[0068] Among them, α is the quadratic penalty factor, λ and λ(t) are both Lagrange multipliers, and <·> represents the inner product operation;

[0069] Solve the equation using the alternating direction multiplication method to obtain u K ,ω K The update formula of and λ is:

[0070]

[0071] Where ω is the frequency, are f(ω), u(ω), λ(ω), u respectively K n+1 (ω),u K The Fourier transform of (ω), ω K n+1 represents the central frequency of mode K at the center of the n+1th iteration, τ represents the penalty factor, that is, the update step size, which is used to control the update speed of the Lagrangian function;

[0072] Set the discrimination accuracy ε>0, and stop the iteration if the following equation is satisfied:

[0073]

[0074] It is worth noting that the VMD decomposition method can adaptively decompose the original runoff time series signal into a specified number of modal components, and each modal component has corresponding frequency and amplitude characteristics, which helps to extract the potential periodic information in the runoff data, reduce the complexity of the data signal, and improve the prediction accuracy of the model;

[0075] It should be noted that the process of performing CEEMDAN decomposition on the relevant water level data includes:

[0076] The input related water level data is recorded as the original signal s(t), and s(t) is decomposed into the intrinsic mode function and the residual signal, where:

[0077]

[0078] Where r(t) is the residual signal, IMF(t) is the intrinsic mode function, and IMF e (t) is the e-th eigenmode function, m is the total number of eigenmode functions;

[0079] The generation process of the intrinsic mode function includes:

[0080] Add noise nos(t) to the original signal s(t), the amplitude of the noise is Where std(s(t)) is the standard deviation of the original signal s(t), and J is the order;

[0081] Perform EMD decomposition on s(t) after adding noise to obtain the first-stage intrinsic mode component; repeat the above steps N times, randomly add different noise each time, and obtain the average value of IMF obtained after all N EMD decompositions Denoted as the first intrinsic mode function; where IMF1 k (t) represents the kth decomposition of the first intrinsic mode function;

[0082] At this time, the process of obtaining the first residual r1(t) is r1(t) = s(t) - IMF1(t); repeat the above steps until it cannot be decomposed;

[0083] It should be noted that the noise added to the original signal is Gaussian white noise;

[0084] It is worth noting that CEEMDAN decomposition (Complete Ensemble Empirical Mode Decomposition with Adaptive Noise) is an advanced time series analysis method, mainly used for signal processing and data decomposition. It was proposed by Torres et al. in 2011 as an improved version of the EMD algorithm, aiming to solve the mode aliasing problem in EMD and improve the stability and accuracy of the decomposition results. The EMD algorithm is empirical mode decomposition, which decomposes signals based on the time scale characteristics of the data itself without presetting any basis functions.

[0085] As a preferred embodiment of the present invention, the process of determining the stationarity of the submodal component includes:

[0086] The significance level threshold Ht=0.05;

[0087] The stability determination is based on the ADF test method, assuming that the submodal component has a unit root, and the statistics Sta and P value of the submodal component are obtained through the ADF test;

[0088] If the P value of the submodal component is ≤Ht, then the submodal component does not have a unit root under the preset significance level threshold, that is, it is a stable submodal component; if the P value of the submodal component is greater than Ht, then the submodal component has a unit root under the preset significance level threshold, that is, it is a non-stationary submodal component;

[0089] It should be noted that the basic idea of ​​the ADF test is to test whether a time series has a unit root by constructing an autoregressive model. The basic model is as follows:

[0090]

[0091] Among them, Δy t =y t -y t-1 Represents the time series y t The first-order difference of , α is the intercept term, β is the coefficient of the time trend term, and γ is the independent variable y lagged by one period. t-1 The coefficient of i is the first-order difference Δy at lag i t-i The coefficient of , p is the lag order, ε t is a random error term.

[0092] The null hypothesis of the ADF test is the existence of a unit root. The p-value of the submodal component is obtained by calculation and compared with the preset significance level threshold to determine whether to reject the null hypothesis. If the p-value is less than or equal to the preset significance level threshold Ht, the null hypothesis can be rejected at the preset significance level, and it is believed that the time series does not have a unit root, that is, it is a stationary series; conversely, if the p-value is greater than the significance level threshold, the null hypothesis cannot be rejected, indicating that the series may be non-stationary. The significance level set in the present invention is 0.05.

[0093] The ADF test method is a method specifically used for testing the stationarity of time series data. Another conventional method is the DF test, which is used to test whether a time series is a pure random sequence; while the ADF test focuses more on testing the stability of the time series, especially when dealing with non-stationary time series data and the presence of high-order lag correlations. The present invention takes into account the situation where the runoff data has high-order lag correlations, that is, the impact of the previous runoff on the current runoff has a lag, so the ADF test method is used.

[0094] As a preferred embodiment of the present invention, the process of dividing the stationary submodal components includes:

[0095] The division of the stationary submodal component is based on the zero-crossing method, and the stationary submodal component is recorded as a detection signal; the frequency of the detection signal when it crosses the zero point is obtained, which is recorded as the zero-crossing rate Zcr, and the zero point is the point where the signal value changes from positive to negative or from negative to positive;

[0096] If the zero-crossing rate Zcr of the detection signal is greater than 0.2, the stationary submodal component corresponding to the detection signal is a high-frequency component; if the zero-crossing rate Zcr of the detection signal is less than or equal to 0.2, the stationary submodal component corresponding to the detection signal is a medium-low frequency component;

[0097] The process of obtaining the trend component and the detail component includes:

[0098] The modal component is recorded as a new detection signal, the zero-crossing rate Zdr of the new detection signal is obtained, and the zero-crossing rate Zdr of the new detection signal is determined; if the new detection signal is a high-frequency component, it is recorded as a detail component; if the new detection signal is a medium-low frequency component, it is recorded as a trend component;

[0099] Step S3: Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component and obtain the average importance score Imp of each feature ave , the average importance score is the average value of the impurity reduction; set the impurity threshold I, select Imp ave Features <I are recorded as key features;

[0100] It can be understood that the trend component is usually the lowest frequency component among all modal components, which represents the long-term trend or background trend of the signal; the detail component is the other modal components except the trend component, which usually contains higher frequency information;

[0101] It should be noted that the present invention uses a random forest algorithm for feature selection. Random forest is an integrated learning method based on decision trees. First, a training data set is used to train a random forest model. During the training process, each tree is built on a subset of the data, and only a part of the features are considered at each node. In the process of building each tree, for each feature, the random forest will calculate the sum of the impurities reduced by the feature at all segmentation points. The greater the average impurity reduction of the feature, the more important the feature. Then, the features are sorted according to their importance. The higher the feature importance score, the greater the contribution of the feature to the model's prediction ability.

[0102] Step S4: obtaining a non-stationary input sequence according to the key features and their corresponding historical data; obtaining a medium-low frequency input sequence according to the daily runoff data corresponding to the medium-low frequency components, and obtaining a high frequency input sequence according to the daily runoff data corresponding to the high frequency components;

[0103] Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix;

[0104] As a preferred embodiment of the present invention, the process of the normalization operation includes:

[0105] Set the normalization range threshold [-1, 1]; record the non-stationary input sequence, the medium-low frequency input sequence and the high frequency input sequence as the matrix to be normalized array; then the process of normalizing the matrix to be normalized is:

[0106]

[0107] Where array′ is the value after normalization operation, min(array) is the minimum value in the matrix to be normalized, and max(array) is the maximum value in the matrix to be normalized;

[0108] Step S5: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction;

[0109] Obtain the prediction results of multi-feature prediction and single-feature prediction, recorded as sub-modal component values; accumulate and reconstruct the sub-modal component values ​​to obtain the daily runoff prediction result;

[0110] It should be noted that the Long Short-Term Memory (LSTM) network is a recurrent neural network with a special structure designed to alleviate the long-term dependency problem and gradient vanishing problem of the RNN network. The core concept of the LSTM network is its unique cell state, also known as the memory cell c. It can capture and store important pieces of information in one time step, and maintain the continuity and accessibility of the information in the following time series, spanning a longer time span than conventional short-term memory. The LSTM model is managed by a gating mechanism consisting of an input gate, a forget gate, and an output gate. The LSTM network structure is as follows: Figure 2 shown.

[0111] The input gate determines which input information should be written into the cell state at the current moment. Its calculation formula is as follows:

[0112] i T =σ(W i *[h T-1 , x T ]+b i );

[0113] Among them, W i represents the weight matrix, b i represents the bias term, [h T-1 , x T ] represents the concatenation vector of the hidden state of the previous time step and the input of the current time step, iT is the input gate activation value of T, σ is the sigmoid activation function, T is the index of the time step, indicating the Tth time point in the sequence;

[0114] The forget gate determines the old information in the cell state that should be discarded. It uses the sigmoid activation function to generate a value between 0 and 1 to proportionally choose to retain or forget past information. The closer the value is to 1, the more it is forgotten; the closer the value is to 0, the more it is retained. Its calculation formula is as follows:

[0115] f T =σ(W f *[h T-1 , x T ]+b f );

[0116] Where W f represents the weight matrix, b f represents the bias term, f T Represents the activation value of the forget gate at time T;

[0117] After that, the cell state is updated according to the output results of the input gate and the forget gate. First, a new candidate cell state C is calculated. T ′, and then combine this candidate state with the result of multiplying the forget gate by the cell state at the previous moment to achieve selective retention and update of information. The calculation formula is as follows:

[0118] C T ′=tanh(W C *[h T-1 , x T ]+b C )

[0119] C T =f T *C T-1 +i T *C T '

[0120] where h T-1 is the hidden state of the previous time step, which is the output of the LSTM unit at the previous moment and is used to capture previous information. C is the weight matrix;

[0121] The output gate determines the current cell state C. T How much information is used to calculate the hidden state h at the current moment T , and used as the input for the next time step. The calculation formula is as follows:

[0122] y T =σ(W y*[h T-1 , x T ]+b y )

[0123] h T =y T *tanh(C T )

[0124] where Y T Represents the output gate value, h T represents the hidden state at the current moment, W y is the weight matrix, y T is the activation value of the output gate, which is actually calculated by the sigmoid function, and y represents the output;

[0125] The LSTM network has a unidirectional structure, and the model can only learn past information. Unlike LSTM, BILSTM consists of LSTM in both the forward and backward directions and has a bidirectional structure. BiLSTM can not only learn past information, but also future information. Through bidirectional historical information data, it can better mine the context information of time series data. The BiLSTM network structure is as follows: Figure 3 shown.

[0126] It should be noted that the gated recurrent unit is a variant of LSTM. It simplifies the network structure while maintaining the same functional effects of LSTM. Compared with LSTM, GRU consists of update gates and reset gates, which can reduce the number of parameters during model training and improve the speed and accuracy of model training. GRU can only learn historical information, while BiGRU is similar to BiLSTM. It can learn past and future information at the same time through forward and backward GRU.

[0127] It should be noted that TCN is a network structure model suitable for time series tasks built on the basis of a one-dimensional convolutional network. To deal with time series problems, TCN modified the convolutional layer structure and introduced causal convolution and dilated convolution. Compared with the traditional CNN network, causal convolution only focuses on the information of the past and current moments, without considering future information, that is, the output value y at time t′ t′ Only the input value R before and at time t′ n ={x0, x1, ...x t′-1 , x t′ The traditional CNN network increases the receptive field by adding pooling layers, but this will cause some information to be lost. The TCN network introduces dilated convolution to increase the receptive field, so that it can obtain a wider range of dependent information without losing information. The expression of the dilated convolution operation is as follows:

[0128]

[0129] Where t′ is the subscript of the sequence element; d is the expansion coefficient; x t′ ∈R n is a sequence element; F(t) is a sequence element x t′ dilated convolution; k′ represents the size of the convolution kernel, f(i′) is the weight value in the convolution kernel, corresponding to the weight of the filter at position i′;

[0130] As the number of network training layers increases, it is easy to encounter the following problems: (1) gradient vanishing or gradient exploding; (2) network degradation. To solve the above problems, TCN adds a residual block. The residual block contains two layers of convolution and nonlinear mapping, and each layer also contains weight regularization and Dropout regularization to avoid gradient vanishing or gradient exploding problems.

[0131] The multi-mode daily runoff prediction system based on VMD-CEEDMAN includes:

[0132] Data collection module: collects historical data, including daily runoff data and water level data; pre-processes and time-sorts the historical data to obtain complete historical time series data;

[0133] Training data set module: Perform VMD decomposition on the daily runoff data to obtain several submodal components of different frequencies; perform stationarity determination on the submodal components according to a preset significance level threshold, divide the submodal components into stationary submodal components and non-stationary submodal components, and divide the stationary submodal components into medium-low frequency components and high frequency components;

[0134] Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components;

[0135] Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component and obtain the average importance score Imp of each feature ave , the average importance score is the average value of the impurity reduction; set the impurity threshold I, select Imp ave Features <I are recorded as key features;

[0136] According to the key features and their corresponding historical data, a non-stationary input sequence is obtained; according to the daily runoff data corresponding to the medium and low frequency components, a medium and low frequency input sequence is obtained, and according to the daily runoff data corresponding to the high frequency components, a high frequency input sequence is obtained;

[0137] Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix;

[0138] Model building module: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction;

[0139] Prediction module: obtain the prediction results of multi-feature prediction and single-feature prediction, recorded as sub-modal component values; accumulate and reconstruct the values ​​of each sub-modal component to obtain the daily runoff prediction result.

[0140] The example verification is as follows:

[0141] 1. Data description:

[0142] This study uses the runoff data of a certain lake mouth area at 08:00 every day from January 1, 2016 to December 31, 2023 as the research object. The data collection frequency is 1 day (a total of 2922 samples). The data includes two indicators: runoff and water level. The missing values ​​are filled by linear interpolation. In order to prevent the model training process from being disturbed by data impurities, this paper divides the data set into training set, validation set and test set in a ratio of 7:2:1. According to the division results, the overall numerical range of runoff Q spans a large range, and the data shows obvious instability, which complicates the prediction of the model. Figure 4 shown.

[0143] 2. Analysis of prediction results:

[0144] In order to verify the effectiveness of the VMD-CEEMDAN-hybrid model combination proposed in this paper, the model is compared with the single BiGRU model, the single BILSTM model, the TCN-BILSTM combination model and the prediction method of this paper. Figure 5 The performance evaluation indicators of each model are shown in the table:

[0145]

[0146] To verify the effect of the proposed method on improving the prediction accuracy of the model, we first take TCN-BILSTM as an example. From the model error indicators, we can see that the hybrid model has lower RMSE and MAE values ​​and a higher determination coefficient R2 than the single model, indicating that the hybrid model has better prediction accuracy and goodness of fit. However, its MAPE value is higher than that of the single model, indicating that the hybrid model is not as good as the single model in predicting extreme values. The RMSE and MAE of the proposed method are reduced by 44.26% and 39.73% respectively compared with the TCN-BiLSTM hybrid model, while MAPE is reduced by nearly 10 percentage points compared with BiGRU, and the determination coefficient R 2 The improvement was 10.03%, indicating that the prediction method is superior to other models in prediction accuracy, goodness of fit, and degree of fit to extreme values.

[0147] The above is a detailed description of an embodiment of the present invention, but the content is only a preferred embodiment of the present invention and cannot be considered to limit the scope of implementation of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the scope of the patent coverage of the present invention.

Claims

1. The multi-mode daily runoff prediction method based on VMD-CEEDMAN is characterized by: The following steps are involved: Step S1: collecting historical data, the historical data including daily runoff data and water level data; preprocessing and time sorting the historical data to obtain complete historical time series data; Step S2: performing VMD decomposition on the daily runoff data to obtain a number of submodal components of different frequencies; performing stationarity determination on the submodal components according to a preset significance level threshold, dividing the submodal components into stationary submodal components and non-stationary submodal components, and dividing the stationary submodal components into medium-low frequency components and high frequency components; Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components; Step S3: Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component and obtain the average importance score Imp of each feature ave , the average importance score is the average value of the impurity reduction; Set the impurity threshold I and select Imp ave Features <I are recorded as key features; Step S4: obtaining a non-stationary input sequence according to the key features and their corresponding historical data; obtaining a medium-low frequency input sequence according to the daily runoff data corresponding to the medium-low frequency components, and obtaining a high frequency input sequence according to the daily runoff data corresponding to the high frequency components; Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix; Step S5: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction; Obtain the prediction results of multi-feature prediction and single-feature prediction, recorded as sub-modal component values; accumulate and reconstruct the sub-modal component values ​​to obtain the daily runoff prediction result; In step S2, the process of determining the stability of the submodal component includes: The significance level threshold Ht=0.05; The stability determination is based on the ADF test method, assuming that the submodal component has a unit root, and the statistics Sta and P value of the submodal component are obtained through the ADF test; If the P value of the submodal component is ≤Ht, then the submodal component does not have a unit root under the preset significance level threshold, that is, it is a stable submodal component; if the P value of the submodal component is greater than Ht, then the submodal component has a unit root under the preset significance level threshold, that is, it is a non-stationary submodal component; In step S4, the normalization process includes: Set the normalization range threshold [-1, 1]; record the non-stationary input sequence, medium-low frequency input sequence and high frequency input sequence as matrices to be normalized array ; Then the process of normalizing the matrix to be normalized is: ； in array ´ is the value after normalization, min( array ) is the minimum value in the matrix to be normalized, max( array ) is the maximum value in the matrix to be normalized.

2. The multi-mode daily runoff prediction method based on VMD-CEEDMAN according to claim 1 is characterized in that: In step S1, the process of preprocessing the historical data includes: Abnormal values ​​in the historical data are obtained, and the abnormal values ​​are blanked to obtain missing values, and the missing values ​​are filled by using a linear interpolation method.

3. The multi-mode daily runoff prediction method based on VMD-CEEDMAN according to claim 1 is characterized in that: In step S2, the process of dividing the stationary submodal components includes: The division of the stationary submodal component is based on the zero-crossing method, and the stationary submodal component is recorded as a detection signal; the frequency of the detection signal when it crosses the zero point is obtained, which is recorded as the zero-crossing rate Zcr, and the zero point is the point where the signal value changes from positive to negative or from negative to positive; If the zero-crossing rate Zcr of the detection signal is greater than 0.2, the stationary submodal component corresponding to the detection signal is a high-frequency component; if the zero-crossing rate Zcr is less than or equal to 0.2, the stationary submodal component corresponding to the detection signal is a medium-low frequency component.

4. The multi-mode daily runoff prediction method based on VMD-CEEDMAN according to claim 3 is characterized in that: In step S2, the process of obtaining the trend component and the detail component includes: The modal component is recorded as a new detection signal, the zero-crossing rate Zdr of the new detection signal is obtained, and the zero-crossing rate Zdr of the new detection signal is determined; if the new detection signal is a high-frequency component, it is recorded as a detail component; if the new detection signal is a medium-to-low-frequency component, it is recorded as a trend component.

5. The multi-mode daily runoff prediction system based on VMD-CEEDMAN is characterized by: include: Data collection module: collects historical data, including daily runoff data and water level data; pre-processes and time-sorts the historical data to obtain complete historical time series data; Training data set module: Perform VMD decomposition on the daily runoff data to obtain several submodal components of different frequencies; perform stationarity determination on the submodal components according to a preset significance level threshold, divide the submodal components into stationary submodal components and non-stationary submodal components, and divide the stationary submodal components into medium-low frequency components and high frequency components; Selecting water level data corresponding to the daily runoff data from the water level data and recording them as relevant water level data; performing CEEDMAN decomposition on the relevant water level data to obtain a number of modal components, and dividing the modal components into trend components and detail components; Obtain the trend component, detail component and non-stationary submodal component, all of which are recorded as non-stationary components; use the random forest algorithm to segment all the features of the non-stationary component, and obtain the average impurity reduction Imp of each feature ave ; Set the impurity threshold I and select Imp ave Features <I are recorded as key features; According to the key features and their corresponding historical data, a non-stationary input sequence is obtained; according to the daily runoff data corresponding to the medium and low frequency components, a medium and low frequency input sequence is obtained, and according to the daily runoff data corresponding to the high frequency components, a high frequency input sequence is obtained; Normalizing the non-stationary input sequence, the medium-low frequency input sequence, and the high frequency input sequence respectively to obtain a non-stationary input matrix, a medium-low frequency input matrix, and a high frequency input matrix; Model building module: Establish a BiGRU neural network with a self-attention mechanism, input the non-stationary input matrix into the BiGRU neural network for multi-feature prediction; input the medium and low frequency input sequence into the BiLSTM neural network for single feature prediction, and input the high frequency input sequence into the TCN+BiLSTM hybrid network for single feature prediction; Prediction module: obtains the prediction results of multi-feature prediction and single-feature prediction, recorded as sub-modal component values; The values ​​of each sub-modal component are accumulated and reconstructed to obtain the daily runoff prediction result; The process of determining the stability of the submodal component includes: The significance level threshold Ht=0.05; The stability determination is based on the ADF test method, assuming that the submodal component has a unit root, and the statistics Sta and P value of the submodal component are obtained through the ADF test; If the P value of the submodal component is ≤Ht, then the submodal component does not have a unit root under the preset significance level threshold, that is, it is a stable submodal component; if the P value of the submodal component is greater than Ht, then the submodal component has a unit root under the preset significance level threshold, that is, it is a non-stationary submodal component; The normalization process includes: Set the normalization range threshold [-1, 1]; record the non-stationary input sequence, medium-low frequency input sequence and high frequency input sequence as matrices to be normalized array ; Then the process of normalizing the matrix to be normalized is: ； in array ´ is the value after normalization, min( array ) is the minimum value in the matrix to be normalized, max( array ) is the maximum value in the matrix to be normalized.

Images