Monthly runoff prediction method based on improved snow ablation optimization algorithm

Abstract

The invention belongs to the technical field of runoff forecasting, and particularly relates to a monthly runoff forecasting method based on an improved snow ablation optimization algorithm, which solves the technical problem of low forecasting precision caused by still existence of high-frequency non-stationary components after monthly runoff primary decomposition, and comprises the following steps: constructing an ISAO; iSAO is utilized to optimize parameters of the CEEMDAN, and the CEEMDAN is utilized to decompose a monthly runoff sequence into a plurality of IMF components and a residual sequence; classifying and reconstructing each IMF component by using a K-means clustering algorithm to obtain a high-frequency sequence, an intermediate-frequency sequence and a low-frequency sequence; iSAO is utilized to optimize parameters of SVMD, and the SVMD is utilized to secondarily decompose the high-frequency sequence into a plurality of IMF components; optimizing the parameters of the BiLSTM by using ISAO, and performing prediction on the residual sequence, the intermediate frequency sequence, the low frequency sequence and the IMF component after secondary decomposition by using the BiLSTM; and integrating the prediction results to obtain a final monthly runoff prediction result. According to the method, the complexity of the original monthly runoff sequence is remarkably reduced, and the monthly runoff prediction precision based on the deep learning model is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of hydrological forecasting, and in particular to the technical field of runoff forecasting, and specifically is a monthly runoff prediction method based on an improved snow melting optimization algorithm. Background Art

[0002] As one of the important research areas of hydrology, medium- and long-term runoff prediction can play an important role in solving the mismatch between natural and human water use and guiding the rational development and management of water resources in river basins. However, under the influence of climate change and human activities, runoff sequences often show nonlinear and non-stationary characteristics, making medium- and long-term runoff prediction extremely challenging. Building an accurate prediction model is of great significance for practical engineering applications.

[0003] At present, the models widely used in runoff prediction mainly include hydrological models based on physical processes and hydrological models based on data-driven. Hydrological models based on physical processes usually combine the influence of meteorological factors and use traditional runoff theory to predict future runoff. However, the shortcomings of this type of model are also obvious, mainly manifested in: meteorological factor data are difficult to obtain, and many empirical parameters need to be determined in traditional hydrological theory. In contrast, data-driven hydrological models do not need to clarify the hydrological physical process, and can be predicted only based on data such as precipitation, evapotranspiration and runoff. This type of model can be further divided into three categories, namely statistical models based on mathematical statistics, machine learning models and deep learning models. Traditional models and machine learning have fast prediction speeds, but ignore the temporal relationship of samples. Therefore, in recent years, deep learning models that consider temporal relationships have been widely used. Long short-term memory neural network LSTM is currently the most popular deep learning network in the field of time series prediction. It inherits the hidden layer structure of RNN sequential connection and replaces simple hidden layer neurons with long short-term memory units, thereby solving the long-term dependency problem caused by gradient explosion and gradient disappearance. However, LSTM can only use information from historical moments, and does not consider information from future moments, so it cannot learn the relationship between the load data in the two time directions. As an improvement of LSTM, the bidirectional long short-term memory network (BiLSTM) can simultaneously consider the forward and backward information in the sequence, providing a more comprehensive time dynamic analysis capability. When dealing with complex sequence prediction problems, it not only has higher prediction accuracy, but is also more effective in understanding the internal structure of the sequence. Therefore, this type of method has broad application prospects.

[0004] In order to further improve the runoff prediction capability based on deep learning models, current related research mainly starts from the following two aspects: On the one hand, starting from reducing the complexity of runoff sequences, the "decomposition-prediction-reconstruction" strategy in the field of time series prediction is introduced. The "decomposition-prediction-reconstruction" strategy can further explore the data characteristics of the runoff series, so it can effectively improve the runoff prediction accuracy based on deep learning. However, related studies have shown that although signal processing can weaken the redundant content in the original sequence and eliminate noise interference, the decomposition characteristics obtained by decomposing the signal by a single decomposition model will still generate high-frequency information containing a large amount of rainfall and runoff sequences. This phenomenon is more obvious in the first few modal components. The components containing high-frequency information will significantly interfere with the prediction and affect the prediction accuracy. To address this problem, the method of secondary decomposition was proposed and applied. This method improves the accuracy and stability of the prediction model and has achieved good results to a certain extent. Unfortunately, in the field of runoff forecasting, secondary decomposition is mainly used to decompose the first high-frequency component after the first decomposition, without considering other components containing high-frequency information. On the other hand, parameter optimization based on optimization algorithm is studied. For a certain runoff sequence decomposition method and deep learning prediction, it is particularly important to choose the right parameters, which is directly related to the runoff sequence decomposition effect and prediction accuracy. Commonly used parameter optimization algorithms include particle swarm optimization PSO, genetic algorithm GA, whale optimization algorithm WOA, snow melt optimization algorithm SAO, etc. Among them, SAO is a novel meta-heuristic algorithm inspired by nature, which shows excellent balance, search efficiency and adaptability when dealing with complex optimization problems, especially multi-peak and high-dimensional problems. However, SAO also has some defects, such as it is easy to fall into local optimality during population initialization and position update. Summary of the invention

[0005] In order to overcome the technical defects in the existing runoff forecasting field that high-frequency non-stationary components still exist after the monthly runoff is decomposed once, and the decomposition and deep learning model parameter optimization are difficult, resulting in low forecasting accuracy, the present invention provides a monthly runoff prediction method based on an improved snow melting optimization algorithm.

[0006] The present invention provides a monthly runoff prediction method based on an improved snow ablation optimization algorithm. S1, when generating the initial population of the snow ablation optimization algorithm SAO, a Sobol sequence is introduced, and an adaptive t-distribution variation strategy is added to the SAO to construct an improved snow ablation optimization algorithm ISAO; S2, the parameters of a fully adaptive noise ensemble empirical mode decomposition algorithm CEEMDAN are optimized by using ISAO, and the monthly runoff sequence is first decomposed into a plurality of intrinsic mode function components IMF and a residual sequence by using CEEMDAN;

[0007] S3, using the K-means clustering algorithm to classify and reconstruct the IMF components obtained in step S2 to obtain a high-frequency sequence, a medium-frequency sequence, and a low-frequency sequence;

[0008] S4. The parameters of the successive variational mode decomposition algorithm SVMD are optimized using the improved snow ablation optimization algorithm ISAO, and the high-frequency sequence is decomposed into multiple IMF components using SVMD.

[0009] S5. Optimize the parameters of the bidirectional long short-term memory network model BiLSTM using the improved snow melting optimization algorithm ISAO, and use BiLSTM to predict the residual sequence, intermediate frequency sequence, low frequency sequence and the IMF component after the secondary decomposition in step S4 respectively;

[0010] S6. Integrate each prediction result to obtain the final monthly runoff prediction result.

[0011] Preferably, the sub-steps of S1 are:

[0012] S11, use Sobol sequence to generate uniform distribution in the unit hypercube [0,1] n The Sobol sequence point set in the improved snow melting optimization algorithm ISAO is used as the initial population, which improves the uniformity of the initial population generated by the traditional SAO using a random algorithm, thereby ensuring that the coverage of the search space is more uniform.

[0013] The generation process of Sobol sequence points depends on the binary representation of integers and a set of direction vectors, where the direction vectors are generated by a specific polynomial. Assume that X i represents the i-th Sobol sequence point, whose j-th dimension coordinate x ij It is expressed as:

[0014]

[0015] In formula (1), v jk is the kth component of the direction vector associated with the jth dimension; a ik is the kth bit in the binary representation of integer i; m is the length of the direction vector;

[0016] S12. An adaptive t-distribution mutation strategy is added to the snow melting optimization algorithm SAO to construct an improved snow melting optimization algorithm ISAO, introducing appropriate randomness and adaptive adjustment to enhance the robustness and global search capability of ISAO; Adaptive t-distribution mutation is a mutation operation method used in optimization algorithms. The core idea of ​​this strategy is to generate the variation amplitude by using the probability density function of the t-distribution, thereby introducing more diverse mutations to avoid the algorithm from falling into the local optimal solution, and adaptively adjust the mutation parameters to meet the needs of the search process; Among them, the adaptive t-distribution mutation strategy generates the variation amplitude by using the probability density function of the t-distribution, and the probability density function of the adaptive t-distribution mutation strategy containing the degree of freedom parameter a is:

[0017]

[0018] In formula (2), Γ is the gamma function; a is the degree of freedom; b is the standardized random variable; and π is 180°.

[0019] Preferably, the sub-steps of S2 are:

[0020] S21, initializing the population size, maximum number of iterations and optimization dimension of the improved snow ablation optimization algorithm ISAO, setting the range of the noise amplitude weight Nstd and the number of noise additions NR in CEEMDAN, and selecting fuzzy entropy as the optimization objective function;

[0021] S22. Add K times of Gaussian white noise with mean 0 to the original sequence x(t) to construct a new sequence x with K trials i (t), where i = 1, 2, ..., K:

[0022] x i (t) = x(t) + εδ i (t) (3),

[0023] In formula (3), ε is the weight coefficient of Gaussian white noise; δ i (t) is the Gaussian white noise generated during the i-th processing;

[0024] S23, for the new sequence x i (t) Perform EMD decomposition to obtain the first IMF component, and take its mean as the first IMF component obtained by CEEMDAN decomposition, that is:

[0025]

[0026] r1(t)=x(t)-IMF1(t)(5), where IMF1(t) is the first IMF component obtained by CEEMDAN decomposition; r1(t) is the undecomposed residual sequence obtained by subtracting the first IMF component from the original sequence, i.e., the residual sequence after the first decomposition of the new sequence;

[0027] S24, after adding specific noise to the residual sequence of the jth stage obtained after decomposition, continue to perform EMD decomposition:

[0028]

[0029] r j (t) = r j-1 (t)-IMF j (t) (7),

[0030] In formula (6) and (7), IMF j (t) is the jth IMF component decomposed by CEEMDAN; E j-1 (δ i (t)) is the j-1th IMF component after EMD decomposition of the new sequence; ε j-1 is the weight coefficient of adding noise to the residual sequence of the j-1th stage by CEEMDAN; r j (t) is the residual sequence of the jth stage;

[0031] S25, repeat step S24 to perform multiple EMD decompositions to obtain the IMF components of the next stage until the termination condition is met, where the termination condition is the residual sequence r of the nth decomposition n (t) is a monotone sequence, the iteration stops and the original sequence x(t) is decomposed into:

[0032]

[0033] S26. Calculate the fitness values ​​of all subsequences decomposed by CEEMDAN according to the fuzzy entropy, and compare them with the fitness values ​​calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameters Nstd and NR of CEEMDAN. Otherwise, the parameters will not be updated, and the number of ISAO iterations will be increased by 1. The optimization iteration will continue until the maximum number of ISAO iterations is reached. ISAO continuously adjusts the parameters of CEEMDAN in an iterative manner, so that the fitness of the decomposition result is gradually optimized, and finally outputs the optimal CEEMDAN parameter combination, and returns the IMF component and residual sequence obtained by decomposition.

[0034] Preferably, the sub-steps of S3 are:

[0035] S31, determine the cluster number K of the K-means clustering algorithm to be 3, corresponding to high-frequency sequence, medium-frequency sequence and low-frequency sequence respectively; set K points as initial cluster center points, and use the IMF components obtained in step S2 as sample points, calculate the Euclidean distance from the sample points to each initial cluster center point, and classify the sample points into the category closest to a certain initial cluster center point to form K initial classifications, where the Euclidean distance calculation formula is:

[0036]

[0037] In formula (9), x i is the i-th variable in sample X, y i is the i-th variable in sample Y;

[0038] S32, determining a criterion function for evaluating clustering results, wherein sample entropy is used as the criterion function to characterize the complexity of the sequence. The closer the sample entropy value is to 0, the closer the sample point is to its initial cluster center point, and the better the clustering effect is;

[0039] S33, recalculate the cluster center μ for each initial classification after division k , through iteration, we continuously search for the appropriate cluster midpoint μ k , so that the sample entropy is minimized, and finally the IMF components decomposed by CEEMDAN are clustered to form high-frequency sequences, medium-frequency sequences and low-frequency sequences.

[0040] Preferably, the sub-steps of S4 are:

[0041] S41, initializing the population size, maximum number of iterations and optimization dimension of the improved snow ablation optimization algorithm ISAO, setting the range of the penalty coefficient α in SVMD, and selecting fuzzy entropy as the optimization objective function;

[0042] S42. Assume that the high frequency sequence f(t) is decomposed into two parts: L-order IMF component and residual sequence:

[0043] f(t)=u L (t)+f r (t) (10),

[0044] In formula (10), the high-frequency sequence includes the unprocessed part f r (t) and the sum of the IMF components obtained previously

[0045] S43, L-order IMF components need to achieve minimization constraints, the constraints are as follows:

[0046]

[0047] In formula (11), J1 is the objective function of the variational problem; u L is the modal function; ω L represents the center frequency of the L-order IMF component; is the partial derivative with respect to time t; δ(t) is the Dirac function; * represents the convolution symbol; is the Hilbert transform; π is 180°; j is the imaginary unit; Indicates that the spectrum is shifted to the baseband;

[0048] S44, in u L (t) Based on the effective component, the residual sequence f r (t) should be minimized. In order to ensure the effective implementation of this constraint, a filter should be selected, where the frequency response of the filter is:

[0049]

[0050] Residual sequence f r The constraint established between (t) and the frequency response of the filter is:

[0051]

[0052] S45. According to the two constraints in formula (11) and formula (13), it is impossible to distinguish the two modes of L-order IMF component and L-1IMF component. Therefore, the idea of ​​J2 constraint is used to select a suitable filter again. The frequency response of the filter is:

[0053]

[0054] The constraints thus established are:

[0055]

[0056] The final constraint ensures that the sequence can be completely reconstructed, and is expressed as follows:

[0057]

[0058] In formula (16), f u (t) is the residual sequence;

[0059] S46. Calculate the fitness values ​​of all subsequences obtained by SVMD decomposition according to the permutation entropy, and compare them with the fitness value calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameter α of SVMD; otherwise, the parameter will not be updated; and the number of ISAO iterations will be increased by 1, and the optimization iteration will continue until the maximum number of ISAO iterations is reached; ISAO continuously adjusts the parameters of SVMD through iterations, so that the fitness of the decomposition result is gradually optimized, and finally, the optimal parameters of SVMD and the IMF components obtained by decomposition are output.

[0060] Preferably, the sub-steps of S5 are:

[0061] S51. Initialize the population size, maximum number of iterations and optimization dimension of the improved snow melting optimization algorithm ISAO, construct the initial matrix of each optimization parameter in BiLSTM, and determine the Dropout layer and the fully connected layer;

[0062] S52, randomly generate an initial population X, and determine the initial velocity and initial position of the particles;

[0063] S53, assigning values ​​to the BiLSTM parameters, training the BiLSTM network under different parameters, and recording each training process;

[0064] S54. BiLSTM is optimized through ISAO. ISAO gradually adjusts the parameters of BiLSTM through multiple iterations, so that BiLSTM moves towards the goal of minimizing the mean square error MSE in each round of iteration. When the maximum number of iterations preset by ISAO is reached, ISAO selects the parameter combination with the smallest MSE as the optimal solution and assigns it to BiLSTM. At this time, the optimized BiLSTM is used to predict the residual sequence, intermediate frequency sequence, low frequency sequence and IMF component after secondary decomposition, respectively, to obtain the final prediction result.

[0065] Compared with the prior art, the technical solution provided by the present invention has the following technical effects: in view of the problem of low forecast accuracy caused by spectral aliasing and high-frequency non-stationary components that still exist after the primary decomposition of monthly runoff, the present invention first proposes an improved snow melting optimization algorithm, and based on the improved optimization algorithm, implements parameter optimization of the primary decomposition CEEMDAN model, the secondary decomposition SVMD model and the deep learning model BiLSTM to ensure that each model exerts the best effect; secondly, in order to fully decompose all high-frequency non-stationary components existing in the original runoff sequence, the K-means clustering algorithm is used to cluster and reconstruct all IMF components generated after the primary decomposition, cluster the high-frequency sequences in each component after the primary decomposition, and then decompose the high-frequency sequences obtained by clustering into multiple IMF components again; the method described in the present invention significantly reduces the complexity of the original monthly runoff sequence, improves the monthly runoff prediction accuracy based on the deep learning model, and is an effective method for monthly runoff prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.

[0067] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0068] Figure 1 The present invention is a structural block diagram of a monthly runoff prediction method based on an improved snow melting optimization algorithm described in an embodiment of the present invention. DETAILED DESCRIPTION

[0069] In order to more clearly understand the above-mentioned objectives, features and advantages of the present invention, the scheme of the present invention will be further described below. It should be noted that the embodiments of the present invention and the features in the embodiments can be combined with each other without conflict.

[0070] In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein; it is obvious that the embodiments in the specification are only part of the embodiments of the present invention, rather than all of the embodiments.

[0071] The specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0072] In one embodiment, Figure 1 As shown, a monthly runoff prediction method based on an improved snow melting optimization algorithm is disclosed, and the steps are:

[0073] S1. Introduce the Sobol sequence when generating the initial population of the snow melting optimization algorithm SAO, and add the adaptive t distribution mutation strategy to the snow melting optimization algorithm SAO to construct the improved snow melting optimization algorithm ISAO;

[0074] S2. The parameters of the fully adaptive noise ensemble empirical mode decomposition algorithm CEEMDAN are optimized using the improved snow ablation optimization algorithm ISAO, and the monthly runoff series is first decomposed into multiple intrinsic mode function components IMF and residual series using CEEMDAN;

[0075] S3, using the K-means clustering algorithm to classify and reconstruct the IMF components obtained in step S2 to obtain a high-frequency sequence, a medium-frequency sequence, and a low-frequency sequence;

[0076] S4. The parameters of the successive variational mode decomposition algorithm SVMD are optimized using the improved snow ablation optimization algorithm ISAO, and the high-frequency sequence is decomposed into multiple IMF components using SVMD.

[0077] S5. Optimize the parameters of the bidirectional long short-term memory network model BiLSTM using the improved snow melting optimization algorithm ISAO, and use BiLSTM to predict the residual sequence, intermediate frequency sequence, low frequency sequence and the IMF component after the secondary decomposition in step S4 respectively;

[0078] S6. Integrate each prediction result to obtain the final monthly runoff prediction result.

[0079] Based on the above embodiment, in a preferred embodiment, the sub-steps of S1 are:

[0080] S11, use Sobol sequence to generate uniform distribution in the unit hypercube [0,1] n The Sobol sequence point set in the improved snow melting optimization algorithm ISAO is used as the initial population, which improves the uniformity of the initial population generated by the traditional SAO using a random algorithm, thereby ensuring that the coverage of the search space is more uniform.

[0081] The generation process of Sobol sequence points depends on the binary representation of integers and a set of direction vectors, where the direction vectors are generated by a specific polynomial. Assume that X i represents the i-th Sobol sequence point, whose j-th dimension coordinate x ij It is expressed as:

[0082]

[0083] In formula (1), v jk is the kth component of the direction vector associated with the jth dimension; a ikis the kth bit in the binary representation of integer i; m is the length of the direction vector;

[0084] S12. An adaptive t-distribution mutation strategy is added to the snow melting optimization algorithm SAO to construct an improved snow melting optimization algorithm ISAO, introducing appropriate randomness and adaptive adjustment to enhance the robustness and global search capability of ISAO; Adaptive t-distribution mutation is a mutation operation method used in optimization algorithms. The core idea of ​​this strategy is to generate the variation amplitude by using the probability density function of the t-distribution, thereby introducing more diverse mutations to avoid the algorithm from falling into the local optimal solution, and adaptively adjust the mutation parameters to meet the needs of the search process; Among them, the adaptive t-distribution mutation strategy generates the variation amplitude by using the probability density function of the t-distribution, and the probability density function of the adaptive t-distribution mutation strategy containing the degree of freedom parameter a is:

[0085]

[0086] In formula (2), Γ is the gamma function; a is the degree of freedom; b is the standardized random variable; and π is 180°.

[0087] Based on the above embodiment, in a preferred embodiment, the sub-steps of S2 are:

[0088] S21, initializing the population size, maximum number of iterations and optimization dimension of the improved snow ablation optimization algorithm ISAO, setting the range of the noise amplitude weight Nstd and the number of noise additions NR in CEEMDAN, and selecting fuzzy entropy as the optimization objective function; in this embodiment, the population size, maximum number of iterations and optimization dimension of ISAO are 5, 100 and 2 respectively, the range of Nstd is [0.2, 0.5], and the range of NR is [50, 500];

[0089] S22. Add K times of Gaussian white noise with mean 0 to the original sequence x(t) to construct a new sequence x with K trials i (t), where i = 1, 2, ..., K:

[0090] x i (t) = x(t) + εδ i (t) (3)

[0091] In formula (3), ε is the weight coefficient of Gaussian white noise; δ i (t) is the Gaussian white noise generated during the i-th processing;

[0092] S23, for the new sequence x i (t) Perform EMD decomposition to obtain the first IMF component, and take its mean as the first IMF component obtained by CEEMDAN decomposition, that is:

[0093]

[0094] r1(t)=x(t)-IMF1(t)(5), where IMF1(t) is the first IMF component obtained by CEEMDAN decomposition; r1(t) is the undecomposed residual sequence obtained by subtracting the first IMF component from the original sequence, i.e., the residual sequence after the first decomposition of the new sequence;

[0095] S24, after adding specific noise to the residual sequence of the jth stage obtained after decomposition, continue to perform EMD decomposition:

[0096]

[0097] r j (t) = r j-1 (t)-IMF j (t) (7),

[0098] In formula (6) and (7), IMF j (t) is the jth IMF component decomposed by CEEMDAN; E j-1 (δ i (t)) is the j-1th IMF component after EMD decomposition of the new sequence; ε j-1 is the weight coefficient of adding noise to the residual sequence of the j-1th stage by CEEMDAN; r j (t) is the residual sequence of the jth stage;

[0099] S25, repeat step S24 to perform multiple EMD decompositions to obtain the IMF components of the next stage until the termination condition is met, where the termination condition is the residual sequence r of the nth decomposition n (t) is a monotone sequence, the iteration stops and the original sequence x(t) is decomposed into:

[0100]

[0101] S26. Calculate the fitness values ​​of all subsequences decomposed by CEEMDAN according to the fuzzy entropy, and compare them with the fitness values ​​calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameters Nstd and NR of CEEMDAN. Otherwise, the parameters will not be updated, and the number of ISAO iterations will be increased by 1. The optimization iteration will continue until the maximum number of ISAO iterations is reached. ISAO continuously adjusts the parameters of CEEMDAN in an iterative manner, so that the fitness of the decomposition result is gradually optimized, and finally outputs the optimal CEEMDAN parameter combination, and returns the IMF component and residual sequence obtained by decomposition.

[0102] Based on the above embodiment, in a preferred embodiment, the sub-step of S3 is:

[0103] S31, determine the cluster number K of the K-means clustering algorithm to be 3, corresponding to high-frequency sequence, medium-frequency sequence and low-frequency sequence respectively; set K points as initial cluster center points, and use the IMF components obtained in step S2 as sample points, calculate the Euclidean distance from the sample points to each initial cluster center point, and classify the sample points into the category closest to a certain initial cluster center point to form K initial classifications, where the Euclidean distance calculation formula is:

[0104]

[0105] In formula (9), x i is the i-th variable in sample X, y i is the i-th variable in sample Y;

[0106] S32, determining a criterion function for evaluating clustering results, wherein sample entropy is used as the criterion function to characterize the complexity of the sequence. The closer the sample entropy value is to 0, the closer the sample point is to its initial cluster center point, and the better the clustering effect is;

[0107] S33, recalculate the cluster center μ for each initial classification after division k , through iteration, we continuously search for the appropriate cluster midpoint μ k , so that the sample entropy is minimized, and finally the IMF components decomposed by CEEMDAN are clustered to form high-frequency sequences, medium-frequency sequences and low-frequency sequences.

[0108] Based on the above embodiment, in a preferred embodiment, the sub-step of S4 is:

[0109] S41, initializing the population size, maximum number of iterations and optimization dimension of the improved snow ablation optimization algorithm ISAO, setting the range of the penalty coefficient α in SVMD, and selecting fuzzy entropy as the optimization objective function; in this embodiment, the population size, maximum number of iterations and optimization dimension of ISAO are 5, 100 and 1 respectively, and the range of α is [10, 20000];

[0110] S42. Assume that the high frequency sequence f(t) is decomposed into two parts: L-order IMF component and residual sequence:

[0111] f(t)=u L (t)+f r (t) (10),

[0112] In formula (10), the high-frequency sequence includes the unprocessed part f r (t) and the sum of the IMF components obtained previously

[0113] S43, L-order IMF components need to achieve minimization constraints, the constraints are as follows:

[0114]

[0115] In formula (11), J1 is the objective function of the variational problem; u L is the modal function; ω L represents the center frequency of the L-order IMF component; is the partial derivative with respect to time t; δ(t) is the Dirac function; * represents the convolution symbol; is the Hilbert transform; л is 180°; j is the imaginary unit; Indicates that the spectrum is shifted to the baseband;

[0116] S44, in u L (t) Based on the effective component, the residual sequence f r (t) should be minimized. In order to ensure the effective implementation of this constraint, a filter should be selected, where the frequency response of the filter is:

[0117]

[0118] Residual sequence f r The constraint established between (t) and the frequency response of the filter is:

[0119]

[0120] S45. According to the two constraints in formula (11) and formula (13), it is impossible to distinguish the two modes of L-order IMF component and L-1IMF component. Therefore, the idea of ​​J2 constraint is used to select a suitable filter again. The frequency response of the filter is:

[0121]

[0122] The constraints thus established are:

[0123]

[0124] The final constraint ensures that the sequence can be completely reconstructed, and is expressed as follows:

[0125]

[0126] In formula (16), f u (t) is the residual sequence;

[0127] S46. Calculate the fitness values ​​of all subsequences obtained by SVMD decomposition according to the permutation entropy, and compare them with the fitness value calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameter α of SVMD; otherwise, the parameter will not be updated; and the number of ISAO iterations will be increased by 1, and the optimization iteration will continue until the maximum number of ISAO iterations is reached; ISAO continuously adjusts the parameters of SVMD through iterations, so that the fitness of the decomposition result is gradually optimized, and finally, the optimal parameters of SVMD and the IMF components obtained by decomposition are output.

[0128] Based on the above embodiment, in a preferred embodiment, the sub-step of S5 is:

[0129] S51. Initialize the population size, maximum number of iterations and optimization dimension of the improved snow melting optimization algorithm ISAO, construct the initial matrix of each optimization parameter in BiLSTM, and determine the Dropout layer and the fully connected layer;

[0130] S52, randomly generate an initial population X, and determine the initial velocity and initial position of the particles;

[0131] S53, assigning values ​​to the BiLSTM parameters, training the BiLSTM network under different parameters, and recording each training process;

[0132] S54. BiLSTM is optimized through ISAO. ISAO gradually adjusts the parameters of BiLSTM through multiple iterations, so that BiLSTM moves towards the goal of minimizing the mean square error MSE in each round of iteration. When the maximum number of iterations preset by ISAO is reached, ISAO selects the parameter combination with the smallest MSE as the optimal solution and assigns it to BiLSTM. At this time, the optimized BiLSTM is used to predict the residual sequence, intermediate frequency sequence, low frequency sequence and IMF component after secondary decomposition, respectively, to obtain the final prediction result.

[0133] The above is only a specific implementation of the present invention, which enables those skilled in the art to understand or implement the present invention. Although detailed descriptions are given with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the aforementioned embodiments, or replace some or all of the technical features therein by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments, and they should all be covered by the protection scope of the claims.

Claims

1. A monthly runoff prediction method based on an improved snow melt optimization algorithm, characterized in that: S1. Introduce the Sobol sequence when generating the initial population of the snow melting optimization algorithm SAO, and add the adaptive t-distribution mutation strategy to the snow melting optimization algorithm SAO to construct the improved snow melting optimization algorithm ISAO; S2, using ISAO to optimize the parameters of the fully adaptive noise ensemble empirical mode decomposition algorithm CEEMDAN, and using CEEMDAN to decompose the monthly runoff series into multiple intrinsic mode function components IMF and residual series for the first time; S3, using the K-means clustering algorithm to classify and reconstruct the IMF components obtained in step S2 to obtain a high-frequency sequence, a medium-frequency sequence, and a low-frequency sequence; S4. The parameters of the successive variational mode decomposition algorithm SVMD are optimized using the improved snow ablation optimization algorithm ISAO, and the high-frequency sequence is decomposed into multiple IMF components using SVMD. S5. Optimize the parameters of the bidirectional long short-term memory network model BiLSTM using the improved snow melting optimization algorithm ISAO, and use BiLSTM to predict the residual sequence, intermediate frequency sequence, low frequency sequence and the IMF component after the secondary decomposition in step S4 respectively; S6. Integrate each prediction result to obtain the final monthly runoff prediction result.

2. The monthly runoff prediction method based on the improved snow melting optimization algorithm according to claim 1 is characterized in that: The sub-steps of S1 are: S11, use Sobol sequence to generate uniform distribution in the unit hypercube [0,1] n The Sobol sequence point set in is used as the initial population of ISAO; The generation process of Sobol sequence points depends on the binary representation of integers and a set of direction vectors, where the direction vectors are generated by a specific polynomial. Assume that X i represents the i-th Sobol sequence point, whose j-th dimension coordinate x ij It is expressed as: In formula (1), v jk is the kth component of the direction vector associated with the jth dimension; a ik is the kth bit in the binary representation of integer i; m is the length of the direction vector; S12. An adaptive t-distribution mutation strategy is added to SAO to construct ISAO. The adaptive t-distribution mutation strategy generates the variation amplitude by using the probability density function of the t-distribution. The probability density function of the adaptive t-distribution mutation strategy containing the degree of freedom parameter a is: In formula (2), Γ is the gamma function; a is the degree of freedom; b is the standardized random variable; π is 180°.

3. The monthly runoff prediction method based on the improved snow melting optimization algorithm according to claim 2 is characterized in that: The sub-steps of S2 are: S21, initializing the population size, maximum number of iterations and optimization dimension of ISAO, setting the range of the noise amplitude weight Nstd and the number of noise additions NR in CEEMDAN, and selecting fuzzy entropy as the optimization objective function; S22. Add K times of Gaussian white noise with mean 0 to the original sequence x(t) to construct a new sequence x with K trials i (t), where i = 1, 2, ..., K: x i (t)=x(t)+εδ i (t)(3), In formula (3), ε is the weight coefficient of Gaussian white noise; δ i (t) is the Gaussian white noise generated during the i-th processing; S23, for the new sequence x i (t) Perform EMD decomposition to obtain the first IMF component, and take its mean as the first IMF component obtained by CEEMDAN decomposition, that is: r1(t)=x(t)-IMF1(t) (5), In formulas (4) and (5), IMF1(t) is the first IMF component obtained by CEEMDAN decomposition; r1(t) is the residual sequence without decomposition obtained by subtracting the first IMF component from the original sequence, that is, the residual sequence after the first decomposition of the new sequence; S24, after adding specific noise to the residual sequence of the jth stage obtained after decomposition, continue to perform EMD decomposition: r j (t)=r j-1 (t)-IMF j (t) (7), In formula (6) and (7), IMF j (t) is the jth IMF component decomposed by CEEMDAN; E j-1 (δ i (t)) is the j-1th IMF component after EMD decomposition of the new sequence; ε j-1 is the weight coefficient of adding noise to the residual sequence of the j-1th stage by CEEMDAN; r j (t) is the residual sequence of the jth stage; S25, repeat step S24 to perform multiple EMD decompositions to obtain the IMF components of the next stage until the termination condition is met, where the termination condition is the residual sequence r of the nth decomposition n (t) is a monotone sequence, the iteration stops and the original sequence x(t) is decomposed into: S26. Calculate the fitness values ​​of all subsequences decomposed by CEEMDAN according to the fuzzy entropy, and compare them with the fitness values ​​calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameters Nstd and NR of CEEMDAN. Otherwise, the parameters will not be updated, and the number of ISAO iterations will be increased by 1. The optimization iteration will continue until the maximum number of ISAO iterations is reached. ISAO continuously adjusts the parameters of CEEMDAN in an iterative manner, so that the fitness of the decomposition result is gradually optimized, and finally outputs the optimal CEEMDAN parameter combination, and returns the IMF component and residual sequence obtained by decomposition.

4. The monthly runoff prediction method based on the improved snow melting optimization algorithm according to claim 1 is characterized in that: The sub-steps of S3 are: S3 1. Determine the number of clusters K of the K-means clustering algorithm to be 3, corresponding to high-frequency sequences, medium-frequency sequences, and low-frequency sequences respectively; set K points as initial cluster centers, and use the IMF components obtained in step S2 as sample points. Calculate the Euclidean distance from the sample points to each initial cluster center, and classify the sample points into the category closest to a certain initial cluster center to form K initial classifications. The Euclidean distance calculation formula is: In formula (9), x i is the i-th variable in sample X, y i is the i-th variable in sample Y; S32, determining a criterion function for evaluating clustering results, wherein sample entropy is used as the criterion function to characterize the complexity of the sequence. The closer the sample entropy value is to 0, the closer the sample point is to its initial cluster center point, and the better the clustering effect is; S33, recalculate the cluster center μ for each initial classification after division k , through iteration, we continuously search for the appropriate cluster midpoint μ k , so that the sample entropy is minimized, and finally the IMF components and clusters obtained by CEEMDAN decomposition are formed into high-frequency sequences, medium-frequency sequences and low-frequency sequences.

5. The monthly runoff prediction method based on the improved snow melting optimization algorithm according to claim 1 is characterized in that: The sub-steps of S4 are: S41, initializing the population size, maximum number of iterations and optimization dimension of the improved snow ablation optimization algorithm ISAO, setting the range of the penalty coefficient α in SVMD, and selecting fuzzy entropy as the optimization objective function; S42. Assume that the high frequency sequence f(t) is decomposed into two parts: L-order IMF component and residual sequence: f(t)=u L (t)+f r (t) (10), In formula (10), the high-frequency sequence includes the unprocessed part f r (t) and the sum of the IMF components obtained previously S43, L-order IMF components need to achieve minimization constraints, the constraints are as follows: In formula (11), J1 is the objective function of the variational problem; u L is the modal function; ω L represents the center frequency of the L-order IMF component; is the partial derivative with respect to time t; δ(t) is the Dirac function; * represents the convolution symbol; is the Hilbert transform; π is 180°; j is the imaginary unit; Indicates that the spectrum is shifted to the baseband; S44, in u L (t) Based on the effective component, the residual sequence f r (t) should be minimized. In order to ensure the effective implementation of this constraint, a filter should be selected, where the frequency response of the filter is: Residual sequence f r The constraint established between (t) and the frequency response of the filter is: S45. According to the two constraints in formula (11) and formula (13), it is impossible to distinguish the two modes of L-order IMF component and L-1IMF component. Therefore, the idea of ​​J2 constraint is used to select a suitable filter again. The frequency response of the filter is: The constraints thus established are: The final constraint ensures that the sequence can be completely reconstructed, and is expressed as follows: In formula (16), f u (t) is the residual sequence; S46. Calculate the fitness values ​​of all subsequences obtained by SVMD decomposition according to the permutation entropy, and compare them with the fitness value calculated by the previous iteration of ISAO. If the fitness value is less than the result of the previous iteration, ISAO will update the parameter α of SVMD; otherwise, the parameter will not be updated; and the number of ISAO iterations will be increased by 1, and the optimization iteration will continue until the maximum number of ISAO iterations is reached; ISAO continuously adjusts the parameters of SVMD through iterations, so that the fitness of the decomposition result is gradually optimized, and finally, the optimal parameters of SVMD and the IMF components obtained by decomposition are output.

6. The monthly runoff prediction method based on the improved snow melting optimization algorithm according to claim 1 is characterized in that: The sub-steps of S5 are: S51. Initialize the population size, maximum number of iterations and optimization dimension of the improved snow melting optimization algorithm ISAO, construct the initial matrix of each optimization parameter in BiLSTM, and determine the Dropout layer and the fully connected layer; S52, randomly generate an initial population X, and determine the initial velocity and initial position of the particles; S53, assigning values ​​to the BiLSTM parameters, training the BiLSTM network under different parameters, and recording each training process; S54. BiLSTM is optimized through ISAO. ISAO gradually adjusts the parameters of BiLSTM through multiple iterations, so that BiLSTM moves towards the goal of minimizing the mean square error MSE in each round of iteration. When the maximum number of iterations preset by ISAO is reached, ISAO selects the parameter combination with the smallest MSE as the optimal solution and assigns it to BiLSTM. At this time, the optimized BiLSTM is used to predict the residual sequence, intermediate frequency sequence, low frequency sequence and IMF component after secondary decomposition, respectively, to obtain the final prediction result.

Images