Runoff forecasting method of EMD-BiLSTM model based on snow melting multi-source factors

Abstract

The invention discloses a runoff forecasting method of an EMD-BiLSTM model based on a snow melting multi-source factor, and the method comprises the steps: firstly building a prediction data set based on multi-source accumulated snow remote sensing data in combination with a hydrological station control basin and a snow melting confluence production process in an alpine region, and carrying out the preprocessing of the data set, performing time-frequency analysis on the runoff sequence by adopting an empirical mode decomposition method, decomposing a complex runoff signal into a plurality of basic mode components, and averaging the mode components to obtain an interpretable runoff change signal; a BiLSTM model is utilized to train a data set, meanwhile, an attention mechanism is introduced, a Bayesian optimization algorithm and a gradient cutting technology are adopted in the training process to determine an optimal hyper-parameter, and the problem of gradient explosion in the back propagation process is prevented. According to the method, high-precision prediction of the daily runoff of the alpine hydrological control basin is achieved through the model, and the method has high stability and generalization ability and is suitable for runoff prediction and water resource management in alpine regions.

Description

Technical Field

[0001] The present invention relates to the field of watershed runoff prediction in alpine regions, and in particular to a runoff forecasting method based on an EMD-BiLSTM model with multi-source factors of snowmelt. Background Art

[0002] Runoff in alpine regions is primarily influenced by a complex combination of meteorological, geographical, and hydrological factors. Snowmelt, in particular, plays a significant role in the formation and variation of runoff. However, due to the nonlinear and non-stationary nature of the runoff process and its dynamic relationship with meteorological conditions, accurate runoff prediction presents significant challenges. Traditional hydrological models, such as physical or statistical models, can simulate runoff changes to a certain extent, but they have limitations when dealing with complex nonlinear and multidimensional data, and are particularly difficult to cope with the complex dynamic changes brought about by multi-source data factors.

[0003] In recent years, with the development of deep learning technology, machine learning models based on time series data have been widely used in the field of hydrological forecasting. Long short-term memory (LSTM) networks have been widely used in hydrological time series forecasting due to their significant advantages in capturing long-term dependencies in time series data. However, unidirectional LSTMs have certain limitations when processing symmetric time series patterns, making it difficult to fully utilize the sequential features of the input sequence. The bidirectional LSTM (BiLSTM) model, by introducing a bidirectional structure, can simultaneously consider the sequential relationships of the input sequence, significantly improving the accuracy of time series forecasting. In addition, empirical mode decomposition (EMD), as an effective signal decomposition method, can effectively handle non-stationary time series, decomposing complex signals into more physically meaningful modal components, thereby improving the interpretability and prediction accuracy of runoff series. Therefore, it is necessary to design a turbine PID parameter optimization method based on a gravitational search algorithm to address the above issues. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a runoff forecasting method based on the EMD-BiLSTM model of multi-source factors of snowmelt, aiming to solve the nonlinear and non-stationary problems of complex snowmelt runoff processes in alpine areas, improve the accuracy and stability of runoff prediction, and achieve efficient prediction of daily runoff in river basins in alpine areas by integrating multi-source factors such as meteorology, snowmelt, and hydrology, combined with deep learning and signal processing technology.

[0005] In order to achieve the above technical effects, the technical solution adopted by the present invention is: A runoff forecasting method based on the EMD-BiLSTM model with multiple snowmelt source factors includes: S1, snowmelt factor extraction: Based on the digital elevation model (DEM) data of the hydrological station, the low-lying areas are eliminated through the depression filling algorithm, and the hydrological basin analysis method is used to determine the scope of the hydrological station's control basin; The long-term cloudless snow remote sensing dataset was used to extract the snowmelt factors of the basin, including snow cover area, snow depth, and total snow water equivalent within the basin. S2, Construction and preprocessing of multi-source factor dataset: Analyze multi-source data, collect and construct prediction data sets for key factors; Fill missing values, process outliers, and normalize data to ensure data consistency and integrity; S3, non-stationary treatment of runoff series: Aiming at the non-stationary nature of runoff time series, the empirical mode decomposition (EMD) method is used to decompose the runoff signal into multiple intrinsic mode components (IMFs). The instantaneous frequency and amplitude information are extracted through the Hilbert-Huang transform to generate a runoff change signal with strong interpretability. S4, construction of BiLSTM model: Based on the preprocessed multi-source snowmelt factor data and the EMD-decomposed runoff series data, a bidirectional long short-term memory network BiLSTM model was constructed. Add an attention mechanism layer to the model to dynamically assign information weights at different time steps, enhancing the model's ability to extract temporal information; S5, model training and optimization: During model training, Bayesian optimization techniques are used to optimize hyperparameters, combined with gradient clipping techniques to prevent gradient explosion and ensure the stability and robustness of model training. S6, Result Analysis and Application: The constructed model is used to predict runoff for the data in the test set, and the predicted results are evaluated against the measured results.

[0006] Preferably, in step S1, the depression filling algorithm is to correct the elevation of low-lying areas in the digital elevation model (DEM) data to ensure the connectivity of the water flow path in the basin.

[0007] Furthermore, in step S1, based on the DEM data of the observation hydrological station, combined with depression filling and hydrological basin analysis methods, the control basin range of the hydrological station is obtained; using the constructed long-term cloudless snow remote sensing dataset, the snow cover, snow depth, snow water equivalent and other parameters of the hydrological control basin of the prediction station are obtained; among them, the long-term cloudless snow remote sensing dataset is derived from the public remote sensing snow datasets such as MODIS, FY3 and ERA5 on the domestic and foreign networks.

[0008] Preferably, the prediction data set of key factors includes previous runoff, basin snow cover area, basin snow depth, total snow water equivalent, precipitation, sunshine duration, basin average temperature, maximum temperature and minimum temperature.

[0009] Preferably, in step S2, the normalization process of the multi-source factor dataset adopts the Z-score normalization method, and the normalization formula is: ; in, is the normalized data, ranging from 0 to 1. is the sequence data value, is the mean of the sequence data, is the standard deviation of the sequence data.

[0010] Preferably, in step S3, the empirical mode decomposition (EMD) process includes: S301, determining the extreme points of the runoff sequence and adding Gaussian white noise to eliminate signal discontinuities; S302, constructing upper and lower envelopes of the runoff sequence using a cubic spline interpolation method, and calculating a mean envelope; The calculation formula of the mean envelope is as follows: ; S303, subtract the mean envelope from the runoff series data to obtain the intrinsic mode function IMF: ; if If the standard deviation of the cycles in the sequence is less than the preset threshold, it is considered that the first intrinsic mode function has been separated Otherwise, continue iterating until the conditions are met; repeat the above process, adding white noise sequences with different normal distributions each time, until all IMF components that meet the conditions are obtained; S304: Perform Hilbert-Huang transform on the residual and each IMF component to extract their instantaneous frequency and amplitude information to reflect the time-frequency characteristics of the runoff series. The residual formula is as follows, and the calculation formulas for the Hilbert-Huang transform, instantaneous frequency, and amplitude information are as follows: ; ; ; ; in, is the residual, is the instantaneous frequency Hilbert transform of the analytical signal, represents the Hilbert transform, is the instantaneous amplitude, is the instantaneous frequency.

[0011] Preferably, in step S4, a bidirectional long short-term memory (BiLSTM) model is constructed based on the preprocessed dataset. The BiLSTM model can simultaneously utilize both forward and backward information in the input sequence, improving the accuracy of time series prediction. Furthermore, to enhance the model's ability to extract important information, an attention mechanism layer is introduced. This mechanism further improves the model's predictive capabilities by assigning weights to information at different time steps. The BiLSTM model includes two LSTM layers, one for processing the forward and one for processing the backward time series.

[0012] Furthermore, the BiLSTM model structure is as follows: Input layer: The input layer receives sequence data, usually represented as a matrix of shape (T, D), where T is the number of time steps and D is the feature dimension of each time step.

[0013] BiLSTM layer: 1. Forward LSTM processes the input sequence from forward to backward to generate the forward hidden state , the forward LSTM output is 2. Backward LSTM processes the input sequence from back to front to generate backward hidden state , the backward LSTM output is , merge the hidden states at each time step to generate a comprehensive hidden state ,in .

[0014] Preferably, the method of introducing the attention mechanism is: Calculate each time step t and calculate the attention weight , the attention weight is calculated by an attention score function , and then normalized by the softmax function: ; ; Among them, c is the context vector, which can be the hidden state of the previous time step or other relevant information; finally, the context vector c is passed to the output layer to generate the final output y; the output layer is a fully connected layer used for regression to obtain the predicted value.

[0015] Preferably, in step S5, the Bayesian optimization technique is used to determine the optimal hyperparameter combination of the BiLSTM model, including the learning rate, the number of LSTM units, and the number of training rounds, and the optimization goal is to minimize the validation set loss.

[0016] Furthermore, Bayesian optimization techniques were used during training to optimize hyperparameters and further enhance model performance. Bayesian optimization intelligently searches for optimal hyperparameter combinations, improving the model's generalization and performance. Furthermore, gradient clipping effectively prevents gradient explosion and ensures a stable training process. By combining Bayesian optimization with gradient clipping, the model demonstrates enhanced robustness and stability in complex data.

[0017] Furthermore, the Bayesian optimization algorithm is used to optimize the hyperparameters in the LSTM model to determine the optimal model structure and parameters. The main steps are as follows: a) Select a set of initial hyperparameters and evaluate their performance on the objective function using accuracy and loss values; b) Build a surrogate model and use Gaussian Process Regression (GPR) to fit known hyperparameter-performance pairs; c) Select the next evaluation point and use the acquisition function to select the next evaluation point. The acquisition function calculation formula is as follows: ; in, is the probability improvement function, that is, the acquisition function. For the best performance currently, A small positive number.

[0018] Preferably, in step S5, the gradient clipping is performed by calculating the L2 norm of the gradient vector and limiting it within a set threshold to prevent gradient explosion.

[0019] Furthermore, gradient clipping is an optimization technique used in deep learning. Its main purpose is to prevent the gradient explosion problem. When the error backpropagates through multiple layers, the gradient may become very large, resulting in excessive weight updates, which in turn affects the stability and convergence of the training process. The specific method is: Assuming there is a model parameter vector θ and its corresponding gradient vector g, the basic operation of gradient clipping is as follows: Calculate the L2 norm of the gradient: ; if If the threshold clip_valu is exceeded, the gradient is scaled: clip_valu; in, is the clipped gradient.

[0020] Preferably, in step S6, the result evaluation is performed based on the coefficient of certainty in the national standard for hydrological information forecasting (GB / T 22482-2008); the degree of agreement between the runoff forecast process and the measured process is evaluated using the coefficient of certainty as an indicator to verify the rationality of the model; the calculation formula of the coefficient of certainty is as follows: ; in the formula is the coefficient of certainty, is the measured value; The mean of the measured values; n is the sequence length.

[0021] The beneficial effects of the present invention are as follows: By incorporating a set of remote sensing parameters for snow cover, the present invention makes up for the shortcomings of traditional hydrological models in factors affecting runoff in alpine regions, significantly improving the accuracy of snowmelt runoff forecasts. Combined with the EMD of runoff time series data, the model can better capture the local characteristics and trends of runoff changes, enhancing the accuracy and stability of the forecast. In addition, a BiLSTM model based on the attention mechanism was built to further enhance the ability to capture key information in time series data and improve the reliability of the forecast. The model was fine-tuned based on the characteristics and seasonal changes of alpine regions to ensure its applicability and robustness in complex environments. By combining EMD, BiLSTM and attention mechanisms, the model has achieved remarkable results in complex snowmelt runoff forecasting scenarios, solving the limitations of traditional models in dealing with multi-source factors, and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is the overall flow chart of the present invention; Figure 2 It is the overall technical flow chart of the present invention; Figure 3 This is a schematic diagram of the BiLSTM principle constructed by the present invention; Figure 4 This is the IMF component graph of the EMD at Dajin Station; Figure 5 This is a three-dimensional display of the EMD at Dajin Station; Figure 6 This is the final prediction result chart of Dajin Station. DETAILED DESCRIPTION

[0023] Example 1: like Figure 1 As shown in the figure, a runoff forecasting method based on the EMD-BiLSTM model with multiple snowmelt sources includes: S1, snowmelt factor extraction: Based on the digital elevation model (DEM) data of the hydrological station, the low-lying areas are eliminated through the depression filling algorithm, and the hydrological basin analysis method is used to determine the scope of the hydrological station's control basin; The long-term cloudless snow remote sensing dataset was used to extract the snowmelt factors of the basin, including snow cover area, snow depth, and total snow water equivalent within the basin. S2, Construction and preprocessing of multi-source factor dataset: Analyze multi-source data, collect and construct prediction data sets for key factors; Fill missing values, process outliers, and normalize data to ensure data consistency and integrity; S3, non-stationary treatment of runoff series: Aiming at the non-stationary nature of runoff time series, the empirical mode decomposition (EMD) method is used to decompose the runoff signal into multiple intrinsic mode components (IMFs). The instantaneous frequency and amplitude information are extracted through the Hilbert-Huang transform to generate a runoff change signal with strong interpretability. S4, construction of BiLSTM model: Based on the preprocessed multi-source snowmelt factor data and the EMD-decomposed runoff series data, a bidirectional long short-term memory network BiLSTM model was constructed. Add an attention mechanism layer to the model to dynamically assign information weights at different time steps, enhancing the model's ability to extract temporal information; S5, model training and optimization: During model training, Bayesian optimization techniques are used to optimize hyperparameters, combined with gradient clipping techniques to prevent gradient explosion and ensure the stability and robustness of model training. S6, Result Analysis and Application: The constructed model is used to predict runoff for the data in the test set, and the predicted results are evaluated against the measured results.

[0024] Preferably, in step S1, the depression filling algorithm is to correct the elevation of low-lying areas in the digital elevation model (DEM) data to ensure the connectivity of the water flow path in the basin.

[0025] Furthermore, in step S1, based on the DEM data of the observation hydrological station, combined with depression filling and hydrological basin analysis methods, the control basin range of the hydrological station is obtained; using the constructed long-term cloudless snow remote sensing dataset, the snow cover, snow depth, snow water equivalent and other parameters of the hydrological control basin of the prediction station are obtained; among them, the long-term cloudless snow remote sensing dataset is derived from the public remote sensing snow datasets such as MODIS, FY3 and ERA5 on the domestic and foreign networks.

[0026] Preferably, the prediction data set of key factors includes previous runoff, basin snow cover area, basin snow depth, total snow water equivalent, precipitation, sunshine duration, basin average temperature, maximum temperature and minimum temperature.

[0027] Preferably, in step S2, the normalization process of the multi-source factor dataset adopts the Z-score normalization method, and the normalization formula is: ; in, is the normalized data, ranging from 0 to 1. is the sequence data value, is the mean of the sequence data, is the standard deviation of the sequence data.

[0028] Preferably, in step S3, the empirical mode decomposition (EMD) process includes: S301, determining the extreme points of the runoff sequence and adding Gaussian white noise to eliminate signal discontinuities; S302, constructing upper and lower envelopes of the runoff sequence using a cubic spline interpolation method, and calculating a mean envelope; The calculation formula of the mean envelope is as follows: ; S303, subtract the mean envelope from the runoff series data to obtain the intrinsic mode function IMF: ; if If the standard deviation of the cycles in the sequence is less than the preset threshold, it is considered that the first intrinsic mode function has been separated Otherwise, continue iterating until the conditions are met; repeat the above process, adding white noise sequences with different normal distributions each time, until all IMF components that meet the conditions are obtained; S304: Perform Hilbert-Huang transform on the residual and each IMF component to extract their instantaneous frequency and amplitude information to reflect the time-frequency characteristics of the runoff series. The residual formula is as follows, and the calculation formulas for the Hilbert-Huang transform, instantaneous frequency, and amplitude information are as follows: ; ; ; ; in, is the residual, is the instantaneous frequency Hilbert transform of the analytical signal, represents the Hilbert transform, is the instantaneous amplitude, is the instantaneous frequency.

[0029] Preferably, in step S4, a bidirectional long short-term memory (BiLSTM) model is constructed based on the preprocessed dataset. The BiLSTM model can simultaneously utilize both forward and backward information in the input sequence, improving the accuracy of time series prediction. Furthermore, to enhance the model's ability to extract important information, an attention mechanism layer is introduced. This mechanism further improves the model's predictive capabilities by assigning weights to information at different time steps. The BiLSTM model includes two LSTM layers, one for processing the forward and one for processing the backward time series.

[0030] Furthermore, the BiLSTM model structure is as follows: Input layer: The input layer receives sequence data, usually represented as a matrix of shape (T, D), where T is the number of time steps and D is the feature dimension of each time step.

[0031] BiLSTM layer: 1. Forward LSTM processes the input sequence from forward to backward to generate the forward hidden state , the forward LSTM output is 2. Backward LSTM processes the input sequence from back to front to generate backward hidden state , the backward LSTM output is , merge the hidden states at each time step to generate a comprehensive hidden state ,in .

[0032] Preferably, the method of introducing the attention mechanism is: Calculate each time step t and calculate the attention weight , the attention weight is calculated by an attention score function , and then normalized by the softmax function: ; ; Among them, c is the context vector, which can be the hidden state of the previous time step or other relevant information; finally, the context vector c is passed to the output layer to generate the final output y; the output layer is a fully connected layer used for regression to obtain the predicted value.

[0033] Preferably, in step S5, the Bayesian optimization technique is used to determine the optimal hyperparameter combination of the BiLSTM model, including the learning rate, the number of LSTM units, and the number of training rounds, and the optimization goal is to minimize the validation set loss.

[0034] Furthermore, Bayesian optimization techniques were used during training to optimize hyperparameters and further enhance model performance. Bayesian optimization intelligently searches for optimal hyperparameter combinations, improving the model's generalization and performance. Furthermore, gradient clipping effectively prevents gradient explosion and ensures a stable training process. By combining Bayesian optimization with gradient clipping, the model demonstrates enhanced robustness and stability in complex data.

[0035] Furthermore, the Bayesian optimization algorithm is used to optimize the hyperparameters in the LSTM model to determine the optimal model structure and parameters. The main steps are as follows: a) Select a set of initial hyperparameters and evaluate their performance on the objective function using accuracy and loss values; b) Build a surrogate model and use Gaussian Process Regression (GPR) to fit known hyperparameter-performance pairs; c) Select the next evaluation point and use the acquisition function to select the next evaluation point. The acquisition function calculation formula is as follows: ; in, is the probability improvement function, that is, the acquisition function. For the best performance currently, A small positive number.

[0036] Preferably, in step S5, the gradient clipping is performed by calculating the L2 norm of the gradient vector and limiting it within a set threshold to prevent gradient explosion.

[0037] Furthermore, gradient clipping is an optimization technique used in deep learning. Its main purpose is to prevent the gradient explosion problem. When the error backpropagates through multiple layers, the gradient may become very large, resulting in excessive weight updates, which in turn affects the stability and convergence of the training process. The specific method is: Assuming there is a model parameter vector θ and its corresponding gradient vector g, the basic operation of gradient clipping is as follows: Calculate the L2 norm of the gradient: ; if If the threshold clip_valu is exceeded, the gradient is scaled: clip_valu; in, is the clipped gradient.

[0038] Preferably, in step S6, the result evaluation is performed based on the coefficient of certainty in the national standard for hydrological information forecasting (GB / T 22482-2008); the degree of agreement between the runoff forecast process and the measured process is evaluated using the coefficient of certainty as an indicator to verify the rationality of the model; the calculation formula of the coefficient of certainty is as follows: ; in the formula is the coefficient of certainty, is the measured value; The mean of the measured values; n is the sequence length.

[0039] Example 2: This embodiment provides a specific implementation process of a runoff forecasting method based on the EMD-BiLSTM model of multi-source snowmelt factors, and its flow chart is shown in Figure 2 , the steps are as follows: Step 1: Snowmelt factor extraction. Based on the DEM data of the observation hydrological station, a depression-filling algorithm is used to eliminate depressions in the hydrological control area to ensure the connectivity of the water flow path. Next, a hydrological basin analysis method is used to calculate and delineate the control basin of the hydrological station. This includes flow direction analysis and runoff accumulation calculations to determine the flow direction and catchment area of ​​each pixel, thereby obtaining an accurate basin boundary.

[0040] After determining the extent of the hydrologically controlled watershed, we used a long-term, cloud-free snow remote sensing dataset to extract snow cover parameters for that watershed. Specifically, we first extracted the distribution of snow cover within the watershed using remote sensing data, and then calculated key parameters such as average snow depth and snow water equivalent using empirical models or remote sensing inversion algorithms. These snow cover factors are core variables influencing snowmelt runoff and provide essential data support for subsequent model predictions.

[0041] Step 2: Construct and preprocess a multi-source factor dataset. First, collect data on key multi-source factors affecting runoff, including previous runoff, snow cover area in the watershed, snow depth, total snow water equivalent, precipitation, sunshine duration, and temperature (maximum, minimum, and average). Data from different sources are processed, missing values ​​are filled, and data is normalized to ensure consistency.

[0042] There are a large number of missing values ​​in the runoff data. This example uses the multiple interpolation method in SPSS. The specific steps are as follows: a) Interpolation: Use Markov Chain Monte Carlo simulation to generate m complete data sets, and interpolate the missing observation rows in the incomplete data sets. The filling value of each interpolation is drawn from a certain distribution; b) Analysis: Each complete dataset was analyzed separately; c) Merge: The results of each imputed dataset are merged to generate the final statistical inference, taking into account the uncertainty of data filling.

[0043] There are a large number of outliers (such as 999999 and NaN) in the weather station observation data. These outliers need to be treated as missing values ​​and corrected using linear interpolation. The processing steps are as follows: Delete records containing outliers; use multiple imputation or linear interpolation methods to fill in outliers to ensure data consistency.

[0044] Since meteorological and hydrological data usually show normal distribution characteristics, this embodiment uses the Z-score standardization method to process the prediction factors. Z-score is a commonly used data standardization technology; in a multi-indicator evaluation system, different indicators may have different effects in the analysis due to differences in their nature, order of magnitude and units. Directly using raw data may cause indicators with larger values ​​to have too high weights in the analysis, while the role of indicators with smaller values ​​is weakened; in order to ensure the scientificity and reliability of the analysis results, the data needs to be standardized. The standardization formula is as follows: ; Among them, is the normalized data, ranging from 0 to 1, is the sequence data value, is the sequence data mean, and is the sequence data standard deviation.

[0045] Step 3: Addressing the nonstationarity of the runoff series. To address the widespread nonstationarity in runoff time series, the EMD method is used to decompose the complex runoff signal into multiple intrinsic mode components (IMFs). First, all extreme points in the runoff series are identified, and Gaussian white noise is added to eliminate signal discontinuities. Next, the upper and lower envelopes are constructed using cubic spline interpolation, and the mean envelope is calculated: ; Subtract the mean envelope from the original runoff series data to obtain the first intrinsic mode function IMF1: ; If the periodic standard deviation of the sequence is less than the preset threshold, the iteration stops; otherwise, the decomposition continues until all IMF components are obtained. After that, the Hilbert-Huang transform is performed on the residual and IMF components to extract their instantaneous frequency and amplitude information and obtain the time-frequency characteristics. The residual formula is as follows: ; ; ; ; in is the residual, is the instantaneous frequency Hilbert transform of the analytical signal, represents the Hilbert transform, is the instantaneous amplitude, is the instantaneous frequency.

[0046] Step 4: Construction of BiLSTM model: like Figure 3 As shown in the figure, a bidirectional long short-term memory (BiLSTM) model was constructed using preprocessed snowmelt data, meteorological data, and empirical mode decomposition (EMD) runoff time series data as model inputs. An attention mechanism layer was introduced into the model design to dynamically assign weights to different time steps. Key parameters of the attention-based BiLSTM model, including the time step, number of hidden layers, learning rate, and training rounds, were initialized to ensure effective model training and optimization.

[0047] The Short-Term Memory (LSTM) model is a special type of Recurrent Neural Network (RNN) designed for time series forecasting. Unlike traditional RNNs, LSTMs incorporate memory cells, replacing standard hidden cells. Each memory cell consists of a state variable, a forget gate, an input gate, and an output gate. By controlling the flow of information and memory between these gates, LSTMs effectively mitigate vanishing and exploding gradients and can capture long-term dependencies and complex patterns in time series. LSTMs are particularly well-suited for processing complex hydrological time series, adapting to diverse data distributions and further enhancing the model's expressiveness through layer stacking.

[0048] Different from traditional RNN, LSTM adds a cell state (C t ) is a key variable used to store long-term memory information, and the cell state is passed through the forget gate (F t )、Input Gate(IN t ) and the output gate (O t ) to update and output. The forget gate controls the cell state information that needs to be removed from the previous moment (t−1) to the current moment (t); the input gate determines the new information that needs to be added; and the output gate is responsible for outputting the current moment's state information. The final cell state combines the current moment's input, information from each gate, the previous moment's hidden layer state, and the cell state information. A hydrological model suitable for runoff simulation in permafrost areas was constructed based on LSTM. The model structure is shown in the figure, where I, II, and III in the LSTM unit represent the forget gate, input gate, and output gate, respectively.

[0049] The LSTM model used in this embodiment includes an input layer, two LSTM layers, and an output layer. The hyperparameter settings of the model are as follows: Figure 4As shown; in order to reduce the risk of overfitting of the model, the specific calculation process of the model is as follows: Input layer: The input layer receives sequence data, usually represented as a matrix of shape (T, D), where T is the number of time steps and D is the feature dimension of each time step.

[0050] BiLSTM layer: 1. Forward LSTM processes the input sequence from forward to backward to generate the forward hidden state , the forward LSTM output is 2. Backward LSTM processes the input sequence from back to front to generate backward hidden state , the backward LSTM output is , merge the hidden states at each time step to generate a comprehensive hidden state ,in .

[0051] Add attention mechanism: 1 Calculate each time step , calculate the attention weight , the attention weight is usually calculated through an attention score function , and then normalized by the softmax function: ; ; Here, c is the context vector, which can be the hidden state of the previous time step or other relevant information. Finally, the context vector c is passed to the output layer to generate the final output y; the output layer is a fully connected layer used for regression to obtain the predicted value.

[0052] Step 5: Model optimization and training. To improve model performance, Bayesian optimization techniques were used during training to optimize the hyperparameters of the BiLSTM model. Bayesian optimization uses intelligent search strategies to optimize hyperparameter combinations and improve the model's generalization capabilities. The main steps are as follows: A) Initially select a set of hyperparameter combinations and evaluate their performance (accuracy, loss value); b) Using Gaussian process regression (GPR) to build a surrogate model and fit the known hyperparameter-performance pairs; c) Use the get function to select the next hyperparameter evaluation point. The calculation formula of the get function is as follows: ; in is the probability improvement function, that is, the acquisition function. For the best performance currently, A small positive number.

[0053] Step 6: Result Analysis and Application. By establishing a good prediction model and inputting the test data, runoff is predicted. Based on the coefficient of certainty in the National Standard for Hydrological Information Forecasting (GB / T 22482-2008), the degree of agreement between the runoff forecast process and the measured process can be used as an indicator to evaluate the prediction results and verify the rationality of the model.

[0054] The calculation formula of the coefficient of certainty is as follows: ; in the formula is the coefficient of certainty; is the measured value; The mean of the measured values; n is the sequence length.

[0055] Example 3: The EMD-BiLSTM model based on multi-source snowmelt factors in this invention is applied to the runoff forecast of the Dajin Station watershed. The calculation process can be divided into six parts. The specific steps of the implementation case are as follows: Step 1: Snowmelt Factor Extraction. Using DEM data from the Dajin Station watershed, a depression-filling algorithm was used to eliminate low-lying areas, ensuring water connectivity and defining the watershed boundaries. Using long-term cloudless snow cover remote sensing data, key snowmelt factors within the Dajin Station watershed were extracted: snow cover area, snow depth, and total snow water equivalent.

[0056] Step 2: Construction and preprocessing of the multi-source factor dataset. Data on previous runoff, snow cover area in the watershed, snow depth, total snow water equivalent, precipitation, sunshine duration, and temperature (maximum, minimum, and average) at Dajin Station were selected, and missing values ​​filled and normalized.

[0057] Step 3: Non-stationary processing of runoff series. The EDM method is used to decompose the non-stationary problem of the original runoff series data. The decomposition results are shown in Figure 4 , the characteristics of the IMF component clearly reflect the time-frequency characteristics of the runoff series.

[0058] Step 4: Build the BiLSTM model. Set key parameters, including time step, number of hidden layers, learning rate, training rounds, etc. An attention mechanism layer is added to enhance the model's ability to capture key time steps, and the Softmax function is used to calculate the attention weight.

[0059] Step 5: Model optimization and training. The model hyperparameter optimization process uses Bayesian optimization technology to optimize the following hyperparameters: Learning rate range: 0.0001, 0.010.0001, 0.010.0001, 0.01; Number of LSTM units: adjustable from 32 to 128; Optimization goal: minimize the validation set loss; During the training process, in order to prevent overfitting, gradient clipping technology is used, and the clipping threshold is set to: clip_value = 1.0 to ensure gradient stability and prevent gradient explosion.

[0060] Step 6: Result analysis and application. The predicted results were compared with the measured runoff data, and the coefficient of certainty in the national standard for hydrological information forecast (GB / T 22482-2008) was used as the evaluation index. In the Dajin Station watershed, the coefficient of certainty of the final prediction result is DC = *; the final prediction result is as follows Figure 5 and Figure 6 As shown in the figure, the predicted results are highly consistent with the actual measurement results, which verifies the rationality and application prospects of the model.

Claims

1. A runoff forecasting method based on the EMD-BiLSTM model of multi-source factors of snowmelt, characterized by: include: S1, snowmelt factor extraction: Based on the digital elevation model (DEM) data of the hydrological station, the low-lying areas are eliminated through the depression filling algorithm, and the hydrological basin analysis method is used to determine the scope of the hydrological station's control basin; The long-term cloudless snow remote sensing dataset was used to extract the snowmelt factors of the basin, including snow cover area, snow depth, and total snow water equivalent within the basin. S2, Construction and preprocessing of multi-source factor dataset: Analyze multi-source data, collect and construct prediction data sets for key factors; Fill missing values, handle outliers and normalize the data; S3, non-stationary treatment of runoff series: Aiming at the non-stationary nature of runoff time series, the empirical mode decomposition (EMD) method is used to decompose the runoff signal into multiple intrinsic mode components (IMFs). The instantaneous frequency and amplitude information are extracted through the Hilbert-Huang transform to generate a runoff change signal with strong interpretability. S4, construction of BiLSTM model: Based on the preprocessed multi-source snowmelt factor data and the EMD-decomposed runoff series data, a bidirectional long short-term memory network BiLSTM model was constructed. Add an attention mechanism layer to the model to dynamically assign information weights at different time steps, enhancing the model's ability to extract temporal information; S5, model training and optimization: During model training, Bayesian optimization techniques are used to optimize hyperparameters, combined with gradient clipping techniques to prevent gradient explosion and ensure the stability and robustness of model training. S6, Result Analysis and Application: The constructed model is used to predict runoff on the data in the test set, and the predicted results are evaluated against the measured results.

2. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S1, the depression filling algorithm corrects the elevation of low-lying areas in the digital elevation model (DEM) data to ensure the connectivity of water flow paths within the basin.

3. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: The prediction data set of key factors includes antecedent runoff, basin snow cover area, basin snow depth, total snow water equivalent, precipitation, sunshine duration, basin mean temperature, maximum temperature and minimum temperature.

4. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S2, the normalization process of the multi-source factor dataset adopts the Z-score normalization method, and the normalization formula is: ； in, is the normalized data, ranging from 0 to 1. is the sequence data value, is the mean of the sequence data, is the standard deviation of the sequence data.

5. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S3, the EMD process includes: S301, determining the extreme points of the runoff sequence and adding Gaussian white noise to eliminate signal discontinuities; S302, constructing upper and lower envelopes of the runoff sequence using a cubic spline interpolation method, and calculating a mean envelope; S303, subtracting the mean envelope from the runoff series data to obtain the intrinsic mode function IMF; S304: Perform Hilbert-Huang transform on the residual and each IMF component to extract their instantaneous frequency and amplitude information to reflect the time-frequency characteristics of the runoff series. The residual formula is as follows, and the calculation formulas for the Hilbert-Huang transform, instantaneous frequency, and amplitude information are as follows: ； ； ； ； in, is the residual, is the instantaneous frequency Hilbert transform of the analytical signal, represents the Hilbert transform, is the instantaneous amplitude, is the instantaneous frequency.

6. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S4, the BiLSTM model includes two LSTM layers, which are used to process forward and backward time series respectively.

7. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 6 is characterized in that: The method of introducing the attention mechanism is: Calculate each time step t and calculate the attention weight , the attention weight is calculated by an attention score function , and then normalized by the softmax function: ； ； Among them, c is the context vector, which can be the hidden state of the previous time step or other relevant information; finally, the context vector c is passed to the output layer to generate the final output y; the output layer is a fully connected layer used for regression to obtain the predicted value.

8. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S5, the Bayesian optimization technique is used to determine the optimal hyperparameter combination of the BiLSTM model, including the learning rate, the number of LSTM units, and the number of training rounds, with the optimization goal being to minimize the validation set loss.

9. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1 is characterized in that: In step S5, the gradient clipping is performed by calculating the L2 norm of the gradient vector and limiting it to a set threshold to prevent gradient explosion.

10. The runoff forecasting method based on the EMD-BiLSTM model of snowmelt multi-source factors according to claim 1, characterized in that: In step S6, the results are evaluated based on the deterministic coefficient of the national standard for hydrological information forecasting. The degree of agreement between the runoff forecast process and the measured process is evaluated using the deterministic coefficient as an indicator to verify the rationality of the model. The calculation formula for the deterministic coefficient is as follows: ； in the formula is the coefficient of certainty, is the measured value; The mean of the measured values; n is the sequence length.

Images