Runoff prediction method based on stl-svmd decomposition and adaptive integration

Abstract

The application discloses a runoff prediction method based on STL-SVMD decomposition and adaptive integration, and comprises the following steps: collecting runoff and multi-source meteorological data; constructing a trend item and a seasonal item responding to long-term climate driving by using a generalized STL method; constructing a multivariate SVMD target function containing a driving alignment item and a physical penalty item, and decomposing runoff residuals into mode subsequences with clear physical correlation on the premise of meeting causality constraints; performing closed-loop screening on the modes based on information retention rate and physical consistency double indexes; constructing a component-level integrated prediction model containing monotonicity, annual integral and extreme value sensitivity constraints; and finally, reconstructing runoff prediction results by using a multi-objective loss function and a time-varying adaptive Bayesian optimization strategy to dynamically adjust model parameters. The application improves the physical consistency and precision of runoff prediction.

Description

Technical Field

[0001] This application belongs to the field of hydrological prediction technology, and in particular to a runoff forecasting method based on STL-SVMD decomposition and adaptive integration. Background Technology

[0002] Runoff forecasting is fundamental to watershed water resources planning and management, flood control, drought relief, and the optimized scheduling of water conservancy projects, and is of great significance for ensuring water security. Especially under the dual influence of climate change and human activities, runoff processes in complex plateau watersheds such as the Yellow River source area exhibit significant non-stationarity and multi-scale variability. Traditional forecasting methods are insufficient to meet increasingly refined management needs. Therefore, developing high-precision runoff forecasting technologies with clear physical mechanisms and adaptability to non-stationary environments has become a research hotspot in the interdisciplinary field of hydrology and artificial intelligence.

[0003] Currently, runoff forecasting methods are mainly divided into two categories: process-driven physical-hydrological models and data-driven machine learning models. Physical models, such as the distributed physical-hydrological model SWAT, rely on complex parameter calibration, while machine learning models, such as support vector machines (SVM) and long short-term memory networks (LSTM), make predictions by mining statistical patterns from historical data. To address the non-stationarity of runoff sequences, current mainstream research tends to adopt a decomposition-ensemble framework. This involves using signal processing techniques such as seasonal trend decomposition (STL), empirical mode decomposition (EMD), or variational mode decomposition (VMD) to decompose the original runoff sequence into several relatively stationary subsequences, establish prediction models for each subsequence, and then superimpose and reconstruct the sequence.

[0004] However, existing decomposition-integration forecasting frameworks still face several deep-seated technical bottlenecks when dealing with complex hydrological processes. These bottlenecks mainly manifest as a disconnect between the signal decomposition process and hydrological physical mechanisms, as well as insufficient model adaptability to environmental changes. Specifically, existing signal decomposition algorithms are typically purely mathematically driven blind decompositions, lacking guidance from external meteorological driving factors and causal constraints. This results in decomposed modes that, while orthogonal in the frequency domain, lack physical interpretability and often fail to meet real-time forecasting requirements due to the use of future information. Simultaneously, downstream forecasting models often lack physical constraints, easily producing predictions that violate water balance or common-sense flow principles. Furthermore, static hyperparameter optimization mechanisms cannot adapt to the drift in hydrological response relationships caused by climate change, leading to a significant decline in predictive performance when facing extreme high-water periods or long-term environmental evolution. Summary of the Invention

[0005] The purpose of this invention is to provide a runoff forecasting method based on STL-SVMD decomposition and adaptive integration to solve the aforementioned problems in the existing technology.

[0006] According to one aspect of this application, a runoff forecasting method based on STL-SVMD decomposition and adaptive ensemble includes:

[0007] Collect runoff time series and multi-source meteorological data for the target watershed;

[0008] Using a generalized multiscale decomposition method that incorporates meteorological driving factors as constraints, multiscale decomposition is performed on runoff time series to separate the trend term series, seasonal term series, and residual term series in response to long-term climate change. The residual term series is further decomposed to extract residual term modal subsequences.

[0009] Candidate feature sets were constructed based on trend term sequences, seasonal term sequences, residual term modal subsequences, and multi-source meteorological data, and key input factor sets were obtained through screening.

[0010] The key input factor set is input into the preset prediction model. After parameter optimization, the predicted values ​​of the trend term, seasonal term, and residual term modal are obtained respectively.

[0011] The predicted values ​​of the trend term, seasonal term, and residual modal term are reconstructed to generate runoff prediction results for the target watershed.

[0012] Beneficial effects: Through the above technical solutions, this invention effectively solves the problems of unclear physical meaning of the decomposition process, poor real-time performance, and difficulty in adapting the model to non-stationary climate change in traditional methods, thereby improving the physical consistency and accuracy of runoff forecasting. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the overall process of the runoff forecasting method based on STL-SVMD decomposition and adaptive integration.

[0014] Figure 2 This is a flowchart illustrating the process of selecting the key input factor set.

[0015] Figure 3 This is a flowchart illustrating the process of predicting the trend term, seasonal term, and residual modal forecast values.

[0016] Figure 4 This is a schematic diagram of the iterative process for solving the multivariate decomposition objective function. Detailed Implementation

[0017] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0018] Example 1 details the complete implementation framework from multi-source data acquisition, data preprocessing, multi-scale decomposition, feature engineering, model prediction to result reconstruction, providing the basic architecture and data flow definition for subsequent more specific algorithm implementations and optimization strategies, such as... Figure 1 As shown.

[0019] Step 101: Collect runoff time series and multi-source meteorological data for the target watershed.

[0020] In this step, runoff time series refers to the time series formed by continuous observation of measured flow data at the outlet section of the target watershed or key hydrological control stations. Specifically, monthly or daily runoff data Q can be obtained from the historical records of hydrological stations. obs This sequence reflects the comprehensive results of runoff generation and confluence processes in the watershed. Multi-source meteorological data refers to observational records of multiple meteorological elements affecting the water cycle processes in the watershed. Specifically, it typically includes precipitation (P), average temperature (T), and other meteorological data. avg Highest temperature T max Minimum temperature T min Relative humidity (RH) and sunshine duration (S) sun And wind speed W speed The data can come from multiple meteorological stations distributed within the watershed, or from gridded reanalysis data released by the meteorological bureau. In practice, to ensure the spatial representativeness of the data, a weighted average of the data from multiple stations within the watershed is usually performed to obtain the average meteorological series for the watershed area. For example, for precipitation, the Thiessen polygon method or inverse distance weighting method can be used to convert point precipitation at each station into isal precipitation. To eliminate the impact of dimensional differences on subsequent model training, all collected time series are usually normalized, for example, by using the min-max normalization method to map the data to the interval between 0 and 1.

[0021] In some alternative implementations, in addition to the conventional meteorological elements mentioned above, interannual climate indices reflecting large-scale climate backgrounds can be collected, such as the El Niño Southern Oscillation Index (ENSO), the Arctic Oscillation Index (AO), or the Pacific Decadal Oscillation Index (PDO), as long-term climate drivers, which helps the model capture the changing patterns of runoff on interdecadal scales. Furthermore, remote sensing data, such as the Normalized Difference Vegetation Index (NDVI) or Snow Cover Rate (SCA), can be introduced as auxiliary input data reflecting changes in underlying surface conditions.

[0022] Step 102: Using a generalized multiscale decomposition method that incorporates meteorological driving factors as constraints, multiscale decomposition is performed on the runoff time series to separate the trend term series, seasonal term series, and residual term series in response to long-term climate change. The residual term series is further decomposed to extract residual term modal subsequences.

[0023] In this step, multi-scale decomposition aims to decouple the complex non-stationary runoff sequence into components with different physical meanings and frequency characteristics, thereby reducing the difficulty of prediction. First, the original runoff sequence Q is decomposed using the first-level decomposition. obs The runoff is decomposed into three main parts: the trend term sequence T, the seasonal term sequence S, and the residual term sequence R. The trend term sequence T characterizes the direction and rate of change in runoff over a long-term scale, typically related to changes in the underlying surface caused by climate change or human activities. The seasonal term sequence S characterizes the periodic fluctuations in runoff throughout the year with the changing seasons, mainly controlled by annual meteorological conditions. The residual term sequence R contains the remaining components after removing the trend and seasonal components, primarily reflecting short-term hydrological events with strong randomness, such as sudden torrential rains and floods or droughts. Mathematically, it usually satisfies an additive model, i.e., Q... obs =T+S+R. For the residual term sequence R, which has complex and unpredictable frequency components, a second-level fine decomposition is performed. Through a further decomposition algorithm, the residual term sequence R is decomposed into several residual term mode subsequences (IMFs) arranged from low to high frequency. k Each modal subsequence represents the component of the residual term in a predetermined frequency band. For example, a low-frequency mode may correspond to the slow fluctuations of groundwater baseflow, while a high-frequency mode may correspond to the rapid response of surface runoff.

[0024] In some optional implementations, the first-level decomposition can employ the Seasonal Trend Decomposition (STL) algorithm, utilizing Loess-weighted local regression to estimate trend and seasonal components, exhibiting strong robustness. The second-level decomposition can employ Variational Mode Decomposition (VMD) or its improved algorithm, Successive Variational Mode Decomposition (SVMD). Compared to Empirical Mode Decomposition (EMD), VMD has a more robust mathematical foundation and can effectively avoid mode aliasing. Furthermore, to adapt to the characteristics of different watersheds, the number of decomposition levels or modes can be fixed or adaptively determined based on sequence complexity.

[0025] Step 103: Construct a candidate feature set based on the trend term sequence, seasonal term sequence, residual term modal subsequence and multi-source meteorological data, and screen to obtain the key input factor set.

[0026] In this step, constructing a candidate feature set aims to uncover potential predictors that may influence changes in each runoff component. Specifically, for each component to be predicted, such as a trend term, seasonal term, or a residual mode, a set of related candidate features can be constructed. Candidate features include not only the component's own historical lag values, such as observations from the past 1 to 24 months, used to capture the autocorrelation of the time series, but also current and historical lag values ​​from multi-source meteorological data, used to capture the lagged effects of meteorological drivers on the runoff response. Furthermore, derived features, such as cumulative precipitation and moving average temperature, can be constructed to reflect the cumulative effects of hydrological processes. After constructing a high-dimensional candidate feature set containing a large number of features, feature filtering is required to remove redundant and irrelevant variables and prevent model overfitting. The filtering process aims to identify the features that contribute most to the prediction objective from the candidate feature set, forming a set of key input factors. This process can be based on statistical indicators, such as the Pearson correlation coefficient, or on machine learning algorithms, such as feature importance assessment using random forests.

[0027] In some alternative implementations, feature selection can be performed in two steps: first, linear correlation analysis is used to eliminate features with extremely weak correlation to the target variable; then, a nonlinear model is used to further select features that contribute to improving prediction accuracy. Furthermore, different sets of key input factors can be constructed for different prediction components. For example, when predicting trend terms, more emphasis may be placed on long-term meteorological accumulation; while when predicting high-frequency residual modes, more emphasis may be placed on current extreme precipitation events.

[0028] Step 104: Input the set of key input factors into the preset prediction model. After parameter optimization, predict the trend term, seasonal term, and residual term modal prediction values ​​respectively.

[0029] In this step, a divide-and-conquer strategy is employed, building separate prediction models for each of the decomposed components. Specifically, independent prediction sub-models are configured for the trend term sequence, the seasonal term sequence, and each residual term modal sub-sequence. The input to each sub-model is a selected set of key input factors specific to that component. The prediction model can be a single machine learning model, such as Support Vector Machine (SVM) or Random Forest (RF), or a deep learning model, such as Long Short-Term Memory (LSTM), or a combination thereof. To ensure optimal model performance, the hyperparameters need to be optimized. The optimization process typically involves dividing the dataset into training and validation sets using historical data. Hyperparameters, such as learning rate, number of iterations, and tree depth, are adjusted by training the model on the training set and evaluating the error on the validation set until the model performs optimally on the validation set. The parameter-optimized model is then used to predict future time steps, outputting the predicted trend term value T. pred Seasonal forecast value S pred And the predicted values ​​of each residual mode (IMF) k_pred .

[0030] In some alternative implementations, different types of model structures can be selected for components with different frequency characteristics. For example, for trend terms with gentle changes, a model with strong modeling ability for long-term dependencies can be used; for high-frequency residual terms with drastic fluctuations, a model sensitive to nonlinear abrupt changes can be used. Furthermore, parameter optimization can employ automated algorithms such as grid search, random search, or Bayesian optimization to improve optimization efficiency and accuracy.

[0031] Step 105: Reconstruct the predicted values ​​of the trend term, seasonal term, and residual modal term to generate the runoff prediction results for the target watershed.

[0032] In this step, reconstruction is the process of recombining the predictions of each component to restore the total runoff prediction. Based on the additive model assumption during decomposition, the reconstruction process typically involves algebraically adding the predictions of all components at their corresponding time points. The specific formula can be expressed as: Q pred =T pred +S pred +IMF k_pred ; where Q pred This step yields the final runoff prediction result. It ensures that the physical meaning of the component predictions can be regressed to the total runoff prediction, resulting in the future runoff values ​​for the target watershed.

[0033] In some alternative implementations, considering the differences in prediction accuracy among the components, the reconstruction process can also employ a weighted summation method, that is, giving greater weight to components with higher prediction accuracy and less weight to components with lower prediction accuracy, in order to further improve the robustness of total runoff prediction. Furthermore, after reconstruction, the prediction results can be post-processed, for example, by correcting negative values ​​to ensure that the predicted runoff values ​​conform to non-negative physical laws.

[0034] Example 2 describes in detail a runoff forecasting scheme based on STL-SVMD decomposition, energy frequency screening, dual feature selection, and component independent modeling. Taking the monthly runoff forecasting at the Tangnaihai hydrological station in the Yellow River source area as an example, the specific implementation details and parameter settings of the scheme are shown.

[0035] Step 201: The STL seasonal-trend decomposition algorithm is used to process the runoff time series, decoupling to obtain the trend term sequence representing long-term changes, the seasonal term sequence representing periodic fluctuations, and the residual term sequence containing random disturbances. The residual term sequence is input into the successive variational mode decomposition model, and the residual term mode subsequences with different center frequencies are extracted successively by constructing a variational constraint problem and solving it iteratively using the alternating direction multiplier method.

[0036] In this step, the preprocessed runoff time series Q from the Tangnaihai hydrological station is first analyzed. obs Perform STL decomposition. STL decomposition is based on the Loess locally weighted regression smoothing method, which contains an inner loop and an outer loop. The inner loop is mainly responsible for smoothing the trend term and calculating the seasonal term, while the outer loop is used to assign smaller weights to outliers to enhance robustness. In the specific implementation, the seasonal window parameter s is set. window For odd numbers of 13 or greater, the trend window parameter t window The value is set based on the sequence length and variation characteristics, for example, 21. After decomposition, it obtains the result satisfying Q. obs =T+S+R relationship has three components. For the high-frequency and non-stationary residual sequence R, successive variational mode decomposition (SVMD) is used for secondary processing. Unlike traditional variational mode decomposition (VMD), which requires a pre-set number of modes K, SVMD estimates one mode and its center frequency from the residual signal each time through successive extraction. Specifically, the algorithm constructs an augmented Lagrangian function containing a quadratic penalty term and Lagrange multipliers, transforming the constrained variational problem into an unconstrained problem, and then uses the alternating direction multiplier method (ADMM) for iterative solution. In the Lth iteration, the algorithm finds the Lth mode u. L and its center frequency ω L This makes the mode compact in the frequency domain and minimizes reconstruction error. This process continues until the energy of the residual signal is less than a preset threshold or the maximum number of modes is reached, resulting in the residual term mode subsequence IMF. k .

[0037] In some optional implementations, to improve the decomposition effect of SVMD, the residual terms R can be pre-whitened or denoised before decomposition. Furthermore, the penalty factor α of SVMD... * The convergence tolerance ε can be adjusted based on the signal-to-noise ratio of the residual term, for example, α. * The possible value is 2000, and ε can be 10. -7 .

[0038] Step 202: Calculate the modal energy density integral and Hilbert spectrum center frequency of each residual modal subsequence; remove subsequences whose modal energy density integral is lower than the preset energy threshold or whose Hilbert spectrum center frequency deviates from the preset physical frequency band, and retain the remaining high-energy effective modal subsequences to construct a candidate feature set.

[0039] In this step, to address the curse of dimensionality and noise introduction caused by excessive modalities generated during decomposition, a screening mechanism based on physical indices is introduced. First, the IMF (Integrated Motion Factor) of each modal subsequence is calculated. k Energy density integral E k Specifically, this can be achieved through the IMF. k The total energy, E, is obtained by integrating or summing the squares of the time series. k =∑(IMF k (t) 2 The Hilbert transform was used to apply the IMF to each modal subsequence. k The signal is processed to obtain its analytical signal, the instantaneous frequency is calculated, and its time average value is taken as the center frequency f of the Hilbert spectrum of this mode. center_k Furthermore, to quantify the information content of a mode, the energy entropy H can be calculated. en First, calculate the energy probability distribution p of the modes. j H is the proportion of the instantaneous energy at point j to the total energy, calculated using the Shannon entropy formula. en H en =-∑(p j ·ln(p j Set the preset energy threshold E. thresh For example, 2% of the total energy or determined based on the inflection point of the energy entropy curve; simultaneously, a reasonable physical frequency band range is set. If the energy density integral E of a certain mode... k Less than E thresh If the noise level is too high, or if it is essentially meaningless white noise, it will be discarded. Only valid modal subsequences (IMFs) that meet the specified conditions will be retained. selected .

[0040] In some alternative implementations, the screening strategy can be combined with correlation analysis, which calculates the correlation coefficient between the modes and the original runoff sequence in addition to energy and frequency, and eliminates spurious modes with correlations significantly lower than the threshold. In the example of Tangnaihai Station, after screening, the top few modes in the energy concentration, such as IMF2, IMF3, and IMF6, are usually retained as effective components.

[0041] Step 203: Filter to obtain the set of key input factors, such as... Figure 2 As shown, the process includes the following steps: calculating the Pearson correlation coefficient between each feature in the candidate feature set and the target runoff sequence, and removing features with absolute correlation coefficients lower than a preset correlation threshold; inputting the retained features into a random forest regression model for training, quantitatively evaluating the feature importance of each feature in the model's decision tree split, and extracting head features based on the importance ranking results; selecting features from the head features that simultaneously meet the Pearson correlation coefficient and feature importance requirements to form a key input factor set that balances linear correlation and nonlinear contribution.

[0042] In this step, a high-dimensional candidate feature set is constructed based on the selected effective modes, trends, and seasonal components, combined with multi-source meteorological data. A dual screening strategy is employed to refine the features.

[0043] The first round of screening is based on linear statistical indicators; specifically, it calculates x for each candidate feature. i The Pearson correlation coefficient r between the target variable y (i.e., the corresponding runoff component) and the target variable y (i.e., the runoff component) xy Set the correlation threshold ρ thresh For example, 0.4. If |r xy If | < 0.4, then the feature is considered to have too weak a linear correlation with the target variable and is therefore removed.

[0044] The second round of selection is based on the contribution of the nonlinear model. Specifically, the feature set selected in the first round is input into the random forest model for training. The random forest consists of multiple decision trees. During training, the reduction in impurity (e.g., the reduction in Gini index) brought by each feature in the node splits of all trees can be calculated, and this is used as the feature importance. i The quantitative metrics are used to sort features from highest to lowest importance and truncate the top features, such as retaining features with a cumulative importance of 90% or the top 50% of features.

[0045] The intersection of the two screening results, i.e. the features that simultaneously meet the requirements of linear correlation and nonlinear importance, is used to form the key input factor set.

[0046] In some alternative implementations, the correlation threshold and feature retention ratio can be dynamically adjusted based on the results of cross-validation. For example, for some highly nonlinear components, the threshold for the Pearson correlation coefficient can be appropriately lowered, relying more on the selection results of the random forest.

[0047] Step 204: The predicted values ​​for the trend term, seasonal term, and residual modal term are obtained respectively, such as... Figure 3 As shown, the process includes the following steps: For the trend term sequence, seasonal term sequence, and each effective modal subsequence, a prediction sub-model containing an XGBoost module and an LSTM module is independently constructed; the XGBoost module is used to mine the nonlinear regression features between meteorological factors and runoff components from the key input factor set, and the LSTM module is used to capture the long and short time series dependence features in the key input factor set; the prediction values ​​of the trend term, seasonal term, and residual modal prediction values ​​of the corresponding components are output by each prediction sub-model to achieve decoupled prediction of hydrophysical processes at different frequencies.

[0048] In this step, a hybrid architecture predictive sub-model is constructed to fully leverage the strengths of different algorithms. For each component to be predicted—the trend term, the seasonal term, or the effective residual mode—the predictive model includes an XGBoost module and an LSTM module.

[0049] The XGBoost module, based on the gradient boosting decision tree principle, excels at processing tabular data and can efficiently uncover complex nonlinear interactions and regression features between meteorological factors and runoff components. The LSTM module, a recurrent neural network, effectively addresses the vanishing gradient problem in long-sequence training by introducing gating mechanisms (input gate, forget gate, output gate), and is adept at capturing long-term and short-term temporal dependencies in time-series data.

[0050] In practical prediction, the outputs of the two modules can be weighted and fused, or a stacking ensemble strategy can be used, where the outputs of the two modules are used as input to the secondary meta-model to obtain the final predicted value of that component. These parallel or serial combinations achieve complementary extraction of the time-dimensional and regression-dimensional features of the runoff process.

[0051] In some alternative implementations, the weights of the two modules can be adjusted based on the characteristics of the components. For example, for seasonal terms with strong periodicity, the LSTM module may have a higher weight; for high-frequency residual terms that are greatly affected by abrupt changes in precipitation, the XGBoost module's ability to capture nonlinear responses may be more critical.

[0052] Step 205: Select the root mean square error or Nash-Sutcliffe efficiency coefficient as a single evaluation index, and construct an objective function to evaluate the generalization error of the prediction sub-model based on the pre-collected validation set data; define the hyperparameter search space of the XGBoost module and LSTM module in the prediction sub-model, and use the Bayesian optimization algorithm to perform probabilistic surrogate model construction and acquisition function calculation within the hyperparameter search space, and carry out iterative optimization guided by minimizing the objective function; obtain the optimal hyperparameter combination that makes the objective function reach the extreme value, and load the optimal hyperparameter combination into the corresponding prediction sub-model, and lock the model state used to perform the final runoff prediction.

[0053] In this step, to address the issues of low efficiency and difficulty in finding the global optimum through manual hyperparameter tuning, a Bayesian optimization algorithm is introduced to automatically optimize the hyperparameters of the XGBoost and LSTM modules. First, a single evaluation metric is determined, typically the root mean square error (RMSE) or the Nash efficiency coefficient (NSE). The objective function is either the RMSE value (to be minimized) or 1-NSE (to be minimized) on the validation set. The search space for hyperparameters is further defined. For XGBoost, optimization variables include the learning rate (e.g., 0.01 to 0.3), the maximum tree depth (e.g., 3 to 10), the subsample rate (e.g., 0.5 to 1.0), and the number of weak learners (e.g., 100 to 1000). For LSTM, optimization variables include the learning rate (e.g., 0.001 to 0.01), the number of hidden layer neurons (e.g., 32 to 256), the batch size (e.g., 16 to 128), and the dropout rate (e.g., 0.1 to 0.5). Bayesian optimization constructs a Gaussian process (GP) as a surrogate model for the objective function. It uses a sampling function (such as the desired improvement in efficiency) to balance exploration (sampling in undetermined regions) and exploitation (sampling in known optimal regions), intelligently selecting the next set of hyperparameters for evaluation. After several iterations, the algorithm converges and returns the hyperparameter combination that optimizes the objective function. The optimal parameters are then assigned to the model, completing training and solidifying the model state for subsequent testing and practical applications.

[0054] In some alternative implementations, to prevent overfitting, a regularization term can be added to the objective function, or the average error of cross-validation can be used as the optimization objective. Compared to the unoptimized model, the Bayesian-optimized combined model achieves a significant improvement in NSE during the validation period at Tangnaihai Station, thus improving prediction accuracy.

[0055] Example 3 describes in detail an improved multi-scale decomposition method, namely the multi-scale co-decomposition algorithm STL-SVMD driven and physically constrained. By introducing multi-source meteorological driving forces and physical constraints in the decomposition process, it solves the problems of unclear physical meaning, difficulty in adapting to climate change, and potential future data leakage in traditional pure signal decomposition methods. The mathematical principles and implementation process of generalized STL decomposition and multivariate SVMD decomposition are elaborated in detail.

[0056] Step 301: Calculate long-term cumulative or moving average meteorological elements based on multi-source meteorological data, and use the cumulative or moving average meteorological elements and time index together as independent variables to construct a trend term sequence that responds to long-term climate change; extract interannual climate indices reflecting large-scale climate background from multi-source meteorological data as modulation factors, decompose seasonal changes into a baseline seasonal model and a dynamic offset controlled by the modulation factor, and construct a seasonal term sequence that evolves with the climate background.

[0057] In this step, the traditional STL algorithm is generalized and improved to make it responsive to external drivers. For the trend term sequence T... * The construction no longer relies solely on time t as a single independent variable. First, it is based on the collected precipitation P and temperature T. air Calculate the cumulative value P over long-term timescales (e.g., 12 months, 24 months, or longer) using the data. acc or moving average T smooth The derived variables reflect the long-term aridity / wetness and energy background of the watershed. The time index t is compared with the long-term driving index X. agg Using all these as inputs, the trend term sequence T is obtained by fitting a multivariate locally weighted regression (such as a locally weighted regression surface, LoessSurface). * The trend term constructed in this way not only contains smoothing features over time, but also adaptively adjusts to long-term warming or aridification trends. For the seasonal term sequence S... * The construction of the model incorporates interannual climate indices (such as the ENSO index) as modulating factors M. clim The seasonal term is decomposed into a fixed baseline seasonal pattern S. base (Month-dependent only) and dynamic offset Δ modulated by climate factors S Specifically, Δ S It can be represented as M clim The function, for example, through a simple linear regression or regression tree model g(M) clim To estimate the final seasonal term sequence S. * =S base +Δ SThe above-mentioned treatment method makes seasonal fluctuations no longer a rigid pattern of repetition year after year, but rather allows for amplitude and phase shifts based on the climate background of the year (such as El Niño or La Niña years), which is more in line with actual hydrological patterns.

[0058] In some alternative implementations, human activity indicators (such as water consumption and land use change index) can be introduced as additional drivers of the trend term to reflect the long-term impact of human activities on the runoff baseline.

[0059] Step 302: Select key meteorological driving variables related to hydrological processes from multi-source meteorological data, and concatenate the residual term sequence with the key meteorological driving variables along the channel dimension to form a multivariate input vector; construct a multivariate decomposition objective function containing driving alignment terms based on the multivariate input vector, configure the driving alignment terms to maximize the correlation between the modal components generated by decomposition and the corresponding key meteorological driving variables, and guide the modal components to align with the physical driving process; decompose the multivariate input vector into multiple pairs of modal components by iteratively solving the multivariate decomposition objective function, and extract residual term modal subsequences from the channels of the corresponding residual term sequence.

[0060] In this step, to overcome the blindness of univariate decomposition, a multivariate SVMD algorithm is constructed. First, a multivariate input vector Z(t) is constructed. This includes the residual term sequence R. * In addition to Z(t), key meteorological driving variables, such as precipitation series P(t) and temperature series T(t), are spliced ​​along the channel dimension. At this point, Z(t) is a multidimensional vector sequence, with its first component being the runoff residual and subsequent components being meteorological drivers. A new decomposition objective function J is further constructed. Besides the bandwidth constraint term (minimizing the modal bandwidth) and reconstruction error term in traditional VMD, a driving alignment term J is explicitly added. drive This term is used to maximize the runoff mode u. k The correlation between J and the corresponding meteorological driving components. Specifically, J drive It can be defined as: J drive =-∑(w k ·Corr(u k ,P') 2 ); where Corr represents the correlation coefficient, w k These are the weighting coefficients, and P' is the precipitation driver after filtering within the corresponding frequency band. By minimizing this objective function, the decomposed modes u are forced to... k Not only is it narrow-band in the frequency domain, but it is also highly correlated with physical drivers such as precipitation in the time domain. This allows each mode to have a more explicit physical interpretation; for example, a certain mode may be identified as a fast response component driven by precipitation.

[0061] In some alternative implementations, the driving variables can be preprocessed, for example, by taking into account confluence lags and using lag-processed precipitation sequences in correlation calculations. Additionally, evaporation ET sequences can be introduced to extract negative response modes associated with the evaporation process.

[0062] Step 303, the multivariate decomposition objective function further includes a physical penalty term, which is configured to impose a numerical penalty on the decomposition results that violate hydrophysical laws. The hydrophysical laws include at least: non-negativity constraints on reconstructed runoff flow, long-term water balance constraints after superposition of decomposition components, or monotonicity constraints on the changing trends between the trend term sequence and key meteorological driving variables.

[0063] In this step, to ensure the physical rationality of the decomposition results, a physical constraint penalty term J is added to the objective function J. phys This term consists of a weighted average of multiple sub-penalty terms. First, there's the non-negativity constraint on flow: since actual river flow cannot be negative, a non-negativity penalty term P is defined. neg In each iteration, the reconstructed total runoff Q is calculated. recon (Trend + Seasonal + Residual Modes) For all time points less than 0, calculate the sum of squares of their amplitudes as a penalty value. Further, there is a long-term water balance constraint: over longer timescales (e.g., several years), the total outflow runoff of the basin should not exceed the limit of total precipitation minus total evaporation. Define the water balance penalty term P. bal The deviation between cumulative runoff and cumulative effective precipitation is calculated, and a penalty is imposed if the deviation exceeds a reasonable range. Finally, there is a trend monotonicity constraint: for the trend term sequence T... * Its direction of change should maintain a certain physical consistency with long-term driving variables (such as cumulative precipitation). For example, if cumulative precipitation increases significantly, the trend term should not show an unexplained sharp decline. Define the monotonicity penalty term P. mono This imposes a penalty on trend changes that violate monotonic relationships. By introducing the above physical constraints, the algorithm automatically avoids mathematically feasible but physically absurd solutions when searching the modal solution space, thus improving the reliability of the decomposition results.

[0064] In some alternative implementations, the physical penalty term can be constructed using a linear rectified function (ReLU), such as P. neg =∑(ReLU(-Q recon ) 2 The weighting coefficient γ of the physical penalty term is used to control the strength of the physical constraints, and usually needs to be adjusted according to the characteristics of the watershed.

[0065] Step 304, iteratively solve the multivariate decomposition objective function, such as... Figure 4As shown, the process includes the following steps: applying causal constraints during the solution process, using a one-sided filter kernel or causal transformation operator with one-sided support characteristics to process the multivariate input vector; ensuring that when calculating the modal component at any current time, only the multivariate input vector data of the current time and historical time are used, and the participation of data from future time times is masked, so as to generate a residual modal subsequence that meets the requirements of real-time forecasting.

[0066] In this step, to make the algorithm applicable to real-time forecasting scenarios, the problem of future data leakage in traditional VMD or Hilbert transform must be eliminated. Traditional methods typically use bilateral filters or Fourier transforms based on the entire sequence, causing the decomposition result at time t to be affected by data after time t. In this embodiment, strict causality constraints are imposed on the SVMD algorithm. Specifically, when constructing the frequency or time domain operators of the variational problem, a one-sided filter kernel h(τ) is forced to be used. This filter kernel satisfies the condition: when τ < 0, h(τ) ≡ 0. When performing convolution operations, the output at the current time depends only on the past and present inputs. Alternatively, a causal Hilbert transform (such as a transform based on the Kramers-Kronig relationship) can be used to compute the analytic signal. In iteratively solving for the mode u k When (t), the update restriction rule only applies to data points with a time index less than or equal to t. Using this method, the generated residual modal subsequence is based on historical information. When new data arrives, there is no need to recalculate the entire historical sequence; only incremental updates are required, thus meeting the requirements of online real-time forecasting.

[0067] In some alternative implementations, a sliding window strategy can be used to achieve near real-time causal decomposition, that is, at each time t, a window of data of length L from the past is extracted for local decomposition, and only the value at the very end of the window is taken as the modal estimate at the current time.

[0068] Example 4 details how to introduce quantitative information retention rate and physical consistency indicators in the multi-scale decomposition process, establish a closed-loop feedback mechanism from screening to decomposition parameter adjustment, solve the problems of open-loop operation of traditional decomposition algorithms, difficulty in detecting loss of effective information or violation of physical logic, and ensure that the decomposition results meet the standards at both the mathematical statistics and hydrophysical levels.

[0069] Step 401: The selected and retained residual modal subsequences are superimposed to generate a reconstructed residual sequence. The comprehensive information retention rate index is calculated based on the reconstructed residual sequence and the residual sequence before decomposition. The comprehensive information retention rate index comprehensively represents the energy retention rate, the degree of temporal correlation and the spectral similarity.

[0070] In this step, the reconstruction operation is first performed, which converts the filtered effective modal subsequences (IMFs) into IMFs. selectedBy accumulating point by point, the reconstructed residual sequence R is obtained. recon To quantitatively assess whether the screening process has lost too much useful information, a composite evaluation index is defined, namely the information retention rate comprehensive index I. info . Specifically, I info It consists of three weighted components: energy retention rate (ER), time-domain correlation coefficient (ρ). time And spectral similarity PSI. Energy retention rate ER is calculated as the reconstructed residual sequence R. recon Energy and residual sequence R * The ratio of their energies, that is:

[0071] ER=∑(R recon 2 ) / ∑(R * 2 Time-domain correlation coefficient ρ time Calculate the reconstructed residual sequence R recon With the residual sequence R * The Pearson correlation coefficient between the waveforms is used to measure their consistency. Spectral similarity (PSI) is determined by comparing the power spectral density (S) of the two waveforms. recon (f) and S * (f) Obtain, for example, the overlap area or correlation of two spectra within the effective frequency band. The final comprehensive information retention rate index I. info =w E *ER+w ρ *ρ time +w P *PSI, where w E w ρ w P These are the weighting coefficients, and their sum is 1. For example, they can be 0.4, 0.3, and 0.3 respectively.

[0072] In some alternative implementations, the comprehensive index of information retention rate I can be used. info It is set as a piecewise function, and when the energy retention rate (ER) is extremely low, it is directly judged as unqualified. Furthermore, in order to capture local information loss, a sliding window can be used to calculate the comprehensive index I of local information retention rate. info The minimum value is taken as the final evaluation result to ensure that the information retention meets the requirements in all time periods.

[0073] Step 402: Compare the comprehensive information retention rate index with the preset retention rate threshold. If the preset retention rate threshold is not met, dynamically adjust the energy threshold or physical frequency band range of modal screening, and re-execute the screening and reconstruction steps until the comprehensive information retention rate index meets the requirements.

[0074] In this step, an inner feedback loop is established for the modal screening process. A preset retention rate threshold T is set.info For example, 0.90. When the calculated I info Less than T info If this occurs, it indicates that the current screening strategy is too aggressive, filtering out too many valid signals. At this point, the algorithm automatically triggers an adjustment mechanism: on the one hand, it can lower the preset energy threshold E. thresh For example, lowering the total energy percentage threshold from 2% to 1% allows more low- and medium-energy modes to enter the effective set; on the other hand, it can broaden the physical frequency band range to accept modes at the edge of the band. After adjusting the parameters, the mode selection and reconstruction calculations are re-executed, and the comprehensive information retention rate index I is re-evaluated. info This continues until the target is met. This ensures that the signals entering subsequent modeling are statistically complete.

[0075] In some alternative implementations, the adjustment strategy can employ a step-by-step search, adjusting by a fixed step size each time until the condition is met. Alternatively, a binary search method can be used to quickly locate the comprehensive information retention rate index I. info The required critical threshold parameters are determined to improve computational efficiency.

[0076] Step 403: Construct a physical consistency index to quantitatively assess the degree of conformity of the reconstructed residual series with respect to hydrological physical laws after superimposing the trend term series and the seasonal term series.

[0077] In this step, in addition to retaining information at the statistical level, a physical-level rationality assessment is also introduced, namely, constructing a physical consistency index I. phys This indicator is based on the reconstructed total runoff Q. recon Calculate, where Q recon =T * +S * +R recon I phys Primarily composed of nonnegative scores S neg Water balance score S bal And monotonicity score S mono Composition. Specifically, the non-negative score S neg Used to penalize negative flow, the calculation formula can be S. neg =exp(-α neg ·∑(max(0,-Q recon )) 2 ), where α neg The adjustment coefficient is 1 when there are no negative values. Water balance score S bal The score is calculated based on the ratio deviation between long-term cumulative runoff and cumulative effective precipitation (precipitation minus evaporation); the smaller the deviation, the higher the score. Monotonicity score S mono Evaluation of trend term sequence T *The score is determined by whether the sign of the correlation with long-term driving variables (such as cumulative precipitation) matches physical expectations; if the sign is opposite, the score drops sharply. The final physical consistency index is I. phys This is the weighted sum of the three scores mentioned above.

[0078] In some alternative implementations, the physical consistency index may also include an extreme value rationality component, i.e., checking the reconstructed total runoff Q. recon To check whether the high water level during the rainy season exceeds a certain multiple of historical observation records, such as 1.5 times, to prevent physically impossible abnormally high values ​​from occurring.

[0079] Step 404: Perform joint verification of information and physical dual indicators. When the comprehensive index of information retention rate or the physical consistency index is lower than the preset standard, trigger the decomposition-side closed-loop optimization mechanism. Jointly backtrack and adjust the penalty term weight, modality number setting or initial value of center frequency in the multivariate decomposition objective function, and use the adjusted parameters to re-perform multivariate decomposition on the residual term sequence until the generated decomposition result simultaneously meets the information retention requirement and the physical consistency requirement.

[0080] In this step, a top-level decomposition-side closed-loop optimization mechanism was constructed. The system simultaneously monitors the comprehensive information retention rate index I. info Physical consistency index I phys Two indicators. If I info Meets the standard but I phys If the value is too low, it indicates that while the decomposition retains information, the physical meaning of the generated modes is chaotic. For example, the trend term may absorb too many high-frequency fluctuations, leading to a loss of monotonicity. In this case, the algorithm backtracks to the SVMD decomposition stage, adjusting the weight γ of the physical penalty term in the objective function to increase its value and strengthen physical constraints; or adjusting the setting of the number of modes K to forcibly merge certain modes with aliased frequencies. If the comprehensive information retention index I... info If the value is too low, it may be necessary to increase the number of modes K or decrease the bandwidth penalty parameter to release more signal details. After adjusting the parameters, restart the entire multivariate decomposition, filtering, and reconstruction process. The dual closed-loop mechanism ensures that the final output modes are both signal-fidelity-preserving and physically interpretable, solving the problem of garbage entering and exiting the signal.

[0081] In some alternative implementations, the parameter tuning for closed-loop optimization can employ a reinforcement learning strategy, integrating the information retention rate index I. info Physical consistency index I phys As a reward function, the agent learns the optimal parameter adjustment direction through trial and error, thereby achieving automated decomposition parameter tuning.

[0082] Example 5 elaborates on the physical constraint prediction model architecture designed for different frequency components, demonstrating how hydrophysical laws can be embedded into the structural design of a machine learning model to form a gray-box prediction system.

[0083] Step 501: Construct a long-term trend branch model and apply monotonicity constraints or water balance constraints to the key input factors that represent the characteristics of long-term driving elements, so as to ensure that the response direction between the predicted value of the trend term and the long-term driving elements conforms to the hydrophysical laws.

[0084] In this step, for the trend term sequence T obtained from the decomposition... * A long-term trend branching model was established. The main inputs to this model include long-term cumulative precipitation P. acc Accumulated temperature T acc And slowly changing driving factors such as land use change index. To prevent the model from learning spurious correlations, such as increased precipitation leading to a decreased runoff trend, monotonicity constraints are imposed on the model structure. Specifically, if a gradient boosting tree model, such as Extreme Gradient Boosting (XGBoost) or LightGBM, is used, monotonicity constraint parameters for features can be set. For example, setting the constraint value corresponding to the cumulative precipitation feature to +1 forces the branch direction of that feature to be positively correlated with the target value during splitting. If a neural network model is used, non-negativity constraints can be introduced into the network weights, or a penalty term for the derivative sign can be added to the loss function. Furthermore, water balance constraints can be introduced, restricting the output trend value from deviating from the theoretical baseline of precipitation minus evaporation over a long period. Through these constraints, the predicted value T output by the trend branching model... pred It always maintains physical rationality in terms of direction and magnitude.

[0085] In some alternative implementations, the trend branch model can employ a combination of linear regression and nonlinear residual correction. The linear component captures the basic water balance relationship, while the nonlinear component handles fine-tuning. The output amplitude of the nonlinear component is strictly limited to ensure the dominance of physical laws.

[0086] Step 502: Construct a seasonal branch model and impose an annual integral constraint on the seasonal forecast values ​​to limit the accumulation and approach to zero of the seasonal forecast values ​​in a single hydrological year, so as to ensure that the seasonal forecasts do not change the long-term water balance.

[0087] In this step, for the seasonal sequence S *A seasonal branching model is established. The model's inputs include the month index, interannual climate indices (such as ENSO), and the current weather conditions. Since the physical definition of the seasonal term is a periodic fluctuation within the year, its main function is the redistribution of water volume within the year, and it should not result in an increase or decrease in net water volume. Therefore, a zero constraint on the intra-annual integral is introduced into the model. In practical implementation, a network layer can be designed, whose output is the seasonal forecast sequence S for the next 12 months. pred_1 To S pred_12 Add a normalization operation to the network output: calculate the mean S of the 12 values. mean Then output S final_j =S pred_j -S mean S pred_j S is the predicted value of the seasonal term for the j-th month. final_j This represents the final seasonal term after zero-mean normalization. After processing, S... final The intra-annual sum of the sequence is always equal to zero. Alternatively, a seasonal penalty term L can be added to the loss function. season =(∑(S pred )) 2 This forces the model to learn a zero-mean fluctuation pattern. It ensures that the seasonal branch only reflects the alternation of wet and dry seasons and does not interfere with the trend branch's prediction of long-term water volume.

[0088] In some alternative implementations, the seasonal branching model can be constructed based on Fourier series, where the predictor variables are not direct runoff values, but rather dynamically changing Fourier coefficients (amplitude and phase). Since the integrals of sine and cosine waves over their periods are zero, this parameterization method naturally satisfies integral constraints.

[0089] Step 503: Construct a high-frequency residual branch model and introduce an attention mechanism to assign high weights to the time periods corresponding to high-energy modes or extreme precipitation events in the key input factors, forming an extreme value sensitive branch that can capture sudden events during the wet season or dry season.

[0090] In this step, the effective modal subsequences (IMFs) after filtering are analyzed. selectedA high-frequency residual branch model is established to capture flood processes caused by sudden heavy rainfall or snowmelt. Traditional LSTM models tend to fit the mean state and are prone to peak smoothing during the high-water season. To address this issue, an attention mechanism is introduced into the model. Specifically, a context vector is constructed, which dynamically calculates attention weights based on the numerical values ​​of input features (such as daily precipitation intensity and previous soil moisture). When extreme precipitation or a surge in previous modal energy is detected in the input, the attention network assigns extremely high weights to that moment, allowing the model to focus on these abrupt signals rather than being overwhelmed by a large amount of data from the normal water season. Furthermore, the output layer of this branch model can be designed with an asymmetric structure, allowing positive (flood) signals to generate large activation values ​​while appropriately suppressing negative fluctuations to simulate the rapid rise and slow fall characteristics of hydrological processes.

[0091] In some alternative implementations, the high-frequency residual branch can employ a quantile regression network, which not only outputs deterministic predicted values ​​but also predictive intervals corresponding to different quantiles (such as the 90th quantile), directly providing the probability range of the occurrence of the high-water season and offering richer information for flood control decisions.

[0092] Example 6 details how to improve the training loss function and optimization strategy to enable the prediction model to simultaneously consider overall accuracy and extreme events, and to have the ability to adapt to environmental changes.

[0093] Step 601: The parameter optimization process of the prediction model is guided by a multi-objective loss function, which is a weighted synthesis of the overall prediction error term, the extreme value prediction error term, and the physical constraint penalty term.

[0094] The overall prediction error term calculated based on the full sample data is used to evaluate the model's fitting accuracy under normal hydrological conditions; the extreme value prediction error term calculated based on the selected high-water or low-water samples is used to enhance the model's ability to capture extreme hydrological events; and the physical constraint penalty term calculated based on the degree of non-negativity violation of the prediction results, the degree of long-term water balance deviation, or the degree of monotonicity violation is used to constrain the physical rationality of the model's output.

[0095] In this step, a composite multi-objective loss function L is defined. multi This is used to guide the joint training or parameter optimization of the aforementioned component-level models. multi It consists of three parts: L multi =α g *L global +β e *L extreme +γ p *L phys Among them, the first item L globalThese are conventional errors calculated based on the entire sample, such as mean squared error (MSE) or 1-NSE, used to ensure the model's ability to fit the overall hydrological process. The second term, L... extreme This is an error term designed for extreme value samples. Specifically, firstly, based on historical flow rankings, the thresholds for high-water seasons (e.g., the top 10% quantile) and low-water seasons (e.g., the bottom 10% quantile) are determined, and the extreme value sample set is then selected. L extreme Only the model's weighted error on these samples is calculated; for example, samples from the high-water season are given a weight five times that of ordinary samples. The third term, L... phys This is a physical constraint penalty term, which calculates the degree to which the predicted result violates physical laws, including the sum of squares of negative flows and long-term water imbalance. During training, the algorithm no longer simply pursues the minimum RMSE, but rather seeks a Pareto optimal solution among these three objectives.

[0096] In some alternative implementations, the weighting coefficient α g β e γ p It can be dynamically changed. For example, during flood season forecasting, the weight β of the extreme value term can be temporarily increased. e This forces the model to prioritize prediction accuracy during the wet season; in the dry season or water resource allocation scenarios, the weight γ of the physical penalty term is increased. p This ensures the accuracy of water balance calculations.

[0097] Step 602: During the operation of the prediction model, the generated historical prediction residual sequence is rolled over and the statistical characteristics of the residual sequence within the sliding time window are calculated. When the statistical characteristics change exceeds the preset significance level, it is determined that the hydrological pattern of the watershed has drifted, and local parameter optimization based on recent data and multi-objective loss function is triggered.

[0098] In this step, to address the concept drift problem caused by climate change and underlying surface alterations, a time-varying adaptive optimization mechanism BO is designed. * During operation, the system maintains a sliding window of length W (e.g., 36 months) and monitors the predicted residual sequence e within this window in real time. t Perform statistical tests periodically (e.g., every 6 months), such as the Kolmogorov-Smirnov test (KS test), to compare the current window residual distribution with the historical baseline distribution to see if there are significant differences. Simultaneously, monitor the multi-objective loss L. multi The recent mean. If the p-value of the KS test is less than 0.05 (significant difference) or L multi If the hyperparameters increase significantly, it indicates that the model is no longer suitable for the current hydrological environment. At this point, a local parameter optimization process is triggered: a local search space is defined centered on the current hyperparameters, and data from the most recent N years (e.g., 5 years) is used as a new training set, with L...multi Given the objective function, initiate the Bayesian optimization algorithm to find a new optimal combination of hyperparameters.

[0099] Step 603: Dynamically adjust the prediction model using the updated hyperparameters generated by local parameter optimization to adapt to the evolution of hydrological response relationships caused by climate change or human activities.

[0100] In this step, once a new hyperparameter combination with better performance is found, the system will perform an online update. This process involves more than just simple parameter replacement; it also includes a smooth transition of model weights. To avoid oscillations in prediction results caused by abrupt parameter changes, the model parameters can be updated using an exponential moving average (EMA), i.e., the new model parameters θ... new =(1-λ)·θ old +λ·θ opt ;θ old For the old model parameters, θ opt To optimize the new model parameters, λ is the update weight of EMA. This mechanism allows the model to smoothly adjust its structural parameters as watershed characteristics slowly evolve, much like a biological adaptation to its environment, maintaining high predictive accuracy throughout its entire lifecycle.

[0101] In some alternative implementations, if a sudden change is detected, such as the construction and impoundment of a large reservoir, the adaptive mechanism can trigger deeper updates, including re-performing feature selection and model retraining, rather than simply adjusting hyperparameters.

[0102] Example 7: Through the actual application case of the Tangnaihai Hydrological Station in the Yellow River source area, this paper verifies in detail the actual effect of the proposed multi-scale co-decomposition and adaptive integrated prediction method of driving and physical constraints from multiple dimensions such as quantitative index comparison, physical consistency verification, and the effectiveness of adaptive mechanism. It provides specific experimental settings, comparative model configurations, and detailed analysis results, and solves the technical problems of insufficient non-stationary signal analysis, lack of physical mechanism, and poor environmental adaptability in the runoff forecasting of complex plateau watersheds in existing technologies.

[0103] Step 701: Set experimental data and evaluation indicators. In this embodiment, the Tangnaihai Hydrological Station, a key control station in the upper reaches of the Yellow River, is selected as the research object. Monthly runoff observation data from January 1960 to December 2020 (61 years) and meteorological driving data within the control basin, including precipitation, temperature, evaporation, and relative humidity, are collected from this station. To verify the generalization ability of the model, the dataset is divided into two parts: the period from 1960 to 2000 is the calibration period, where the data is used for model parameter calibration and Bayesian optimization; the period from 2001 to 2020 is the model validation period. The evaluation indicators selected are the Nash efficiency coefficient (NSE), root mean square error (RMSE), and mean absolute error (MAE), which are commonly used in hydrology. Among them, NSE measures the degree of agreement between the predicted process line and the observed process line; the closer it is to 1, the better the effect. RMSE and MAE measure the magnitude of the deviation of the predicted values; the smaller the value, the higher the accuracy. In addition, the relative error R during the high-water season is introduced. peak The average deviation R during the dry season low These are used to evaluate the model's ability to capture extreme flood and drought events, respectively.

[0104] Step 702 involves configuring comparative models and analyzing quantitative results. Three representative comparative models are constructed: the first group is the distributed physical hydrological model SWAT, representing the benchmark of traditional mechanistic models; the second group is a single XGBoost model, representing a pure data-driven benchmark without multi-scale decomposition; and the third group is the standard STL-SVMD-LSTM combined model, representing existing technologies that incorporate conventional decomposition but lack physical constraints and driver alignment. The proposed multi-scale co-decomposition and adaptive ensemble prediction method model for driver and physical constraints is denoted as STL-SVMD-ML. During the validation period, the SWAT model achieved an NSE of 0.82 and an RMSE of 331.50m. 3 The NSE of the single XGBoost model is only 0.78, and the RMSE is 360.20m, which is limited by the uncertainty of parameter calibration. 3 The standard STL-SVMD-LSTM model, with a flow rate of [value] / s, struggles to simultaneously fit both the long-term trend and high-frequency fluctuations of runoff. The standard STL-SVMD-LSTM model improves performance through decomposition techniques, achieving an NSE of 0.86 and reducing the RMSE to 196.11m. 3 / s, but its performance still fluctuates in years of abrupt climate change. The STL-SVMD-ML model proposed in this invention achieves better prediction accuracy, with an NSE of 0.92 and an RMSE further reduced to 155.20m during the validation period. 3 / s, MAE is 107.46m 3 / s. Compared with the standard decomposition model, this invention shows a significant reduction in RMSE, demonstrating the effectiveness of introducing a driving alignment term in the decomposition stage and a physical constraint loss function in the prediction stage.

[0105] Step 703 verifies physical consistency and dry season stability. For dry season prediction, a key concern in water resource management, the model's performance during the winter receding water phase is analyzed. Purely data-driven models (such as XGBoost) tend to output negative values ​​or violently oscillating values ​​that violate physical laws during the dry season because the model learns noise rather than physical laws in the low-value range. The model of this invention introduces a physical penalty term in the decomposition stage and explicitly adds non-negativity constraints and water balance constraints to the loss function in the prediction stage. During the dry season from 2000 to 2002, the predicted output of the model of this invention was consistently greater than zero, and the receding water curve conformed to the physical characteristics of exponentially decaying baseflow, without exhibiting non-physical random oscillations. Physical consistency index I phys The evaluation results show that the physical violation rate of the model of the present invention is close to zero. That is to say, the physical constraint mechanism proposed in this invention restricts the degree of freedom of the data-driven model to the range allowed by hydrological physical laws, thereby improving the credibility of the prediction results.

[0106] Step 704: Verify the adaptability of the adaptive optimization mechanism to non-stationary environments. Since the beginning of the 21st century, the runoff characteristics in the Yellow River source area have changed significantly due to climate warming and humidification. To verify the effectiveness of the time-varying adaptive optimization mechanism, an ablation experiment was conducted. The static optimization model (where hyperparameters are determined once during training and not adjusted again) and the adaptive optimization model (using the sliding window detection and local update mechanism of this invention) were compared. The static model's predictive performance showed a declining trend year by year after 2010, with the average NSE dropping to 0.85 from 2014 to 2020. The adaptive model of this invention, by monitoring changes in residual statistical characteristics, triggered local parameter updates twice in 2009 and 2013, automatically adjusting the tree depth of XGBoost and the forget gate bias of LSTM, adapting to the new rainfall-runoff response relationship. The adaptive model maintained an average NSE of 0.91 from 2014 to 2020, indicating that the mechanism can effectively address the non-stationarity of hydrological sequences caused by climate change and extend the lifespan of the predictive model.

[0107] Example 8 details how to incorporate multi-source meteorological driving elements as constraints into the generalized multi-scale decomposition framework during the driving constraint decomposition process, thus solving the problem that traditional decomposition methods have unclear physical meanings and cannot respond to external driving changes.

[0108] Step 801: Calculate long-term cumulative or moving average meteorological elements based on multi-source meteorological data, and use the cumulative or moving average meteorological elements and time index together as independent variables to construct a trend term sequence in response to long-term climate change.

[0109] In this step, the driving constraint decomposition first requires constructing a trend term sequence T that can reflect long-term climate driving factors. * Unlike traditional STL decomposition, which only uses time t as an independent variable for local regression, this embodiment introduces a long-term driving index X. agg As an additional constraint dimension, specifically, based on the collected precipitation sequence P(t), an exponentially decaying weighted cumulative precipitation index is constructed, calculated using the following formula: X agg (t) = ∑P(t-i')*exp(-i' / τ); where τ is the decay time constant, usually taken as 12 months, characterizing the length of the memory effect of precipitation on runoff, and the variable i' is the summation index, representing the number of time steps (lag order) for backtracking. The design of the exponential decay weight takes into account the hydrophysical law that the contribution of recent precipitation to the current runoff trend is greater than that of long-term precipitation. At the same time, the 12-month moving average T of the temperature series T(t) is calculated. smooth (t) is used to characterize the long-term thermal background of the watershed. The time index t and the cumulative precipitation index X are used... agg (t) and average temperature T smooth (t) is used as a three-dimensional input variable, and a trend term, namely T, is constructed using a multivariate locally weighted regression algorithm. * (t)=LOESS 3D (t,X agg (t),T smooth (t)). The bandwidth parameter for local regression is set to 24 months in the time dimension and 20% of the data range in the driving dimension to balance smoothing effects with local adaptability. The trend component constructed in this way not only retains the smoothing characteristics of the time dimension but also dynamically responds to changes in cumulative precipitation and temperature trends, giving the trend component a clear climate-driven physical meaning.

[0110] In some alternative implementations, the construction of long-term driving indicators can also consider the previous cumulative state of soil moisture. For example, the 90-day moving average of the soil moisture index (SMI) can be introduced as a fourth input dimension to reflect the moderating effect of watershed water storage on baseflow. Furthermore, for high-altitude watersheds significantly affected by snowmelt runoff, the cumulative value of snow depth (SD) can be incorporated into the driving indicator system. In addition, the form of the cumulative function can be adjusted according to watershed characteristics, for example, using power function decay instead of exponential decay, or using a piecewise linear weighting function to assign different cumulative weights at different time scales.

[0111] Step 802: Extract interannual climate indices reflecting large-scale climate background from multi-source meteorological data as modulation factors, and decompose seasonal variation into a baseline seasonal model and a dynamic offset controlled by the modulation factors.

[0112] In this step, the seasonal term sequence S* The construction of the model employs a strategy combining a baseline model and dynamic modulation. First, the baseline seasonal model S is obtained by calculating the multi-year average of historical data. base (m), where m represents the month, ranging from 1 to 12. In the specific calculation, for each month m, S... base (m) can be calculated as the median of runoff anomalies for all historical years in that month to enhance robustness to outliers. Interannual climate indices are introduced as modulating factors, including the El Niño Southern Oscillation Index (ENSO). index (t), Arctic Oscillation Index (AO) index (t) and the Pacific Decadal Oscillation Index (PDO) index (t). Dynamic offset Δ S (t) is calculated using a modulation function g, which can be a linear regression model or a regression tree model. In the case of a linear model,

[0113] Δ S (t)=β1·ENSO index (t)·S base (m(t))+β2·AO index (t)·sin((2π·m(t)) / 12), where β1 and β2 are regression coefficients obtained through training with historical data. The final dynamic seasonal term sequence S * (t)=S base (m(t))+Δ S (t). This design allows seasonal fluctuations to be adaptively adjusted in amplitude and phase according to the large-scale climate background of the year. For example, in a strong El Niño year, the peak summer runoff may occur earlier and be larger.

[0114] In some alternative implementations, the modulation function g can employ a nonlinear model, such as a neural network or support vector regression, to capture the complex nonlinear relationship between climate indices and seasonal variations. Furthermore, regional climate indicators, such as the East Asian Monsoon Index (EAMI) or the Tibetan Plateau Snow Index, can be introduced to better reflect the climate-driven characteristics of a specific watershed. For watersheds with less pronounced seasonality, the baseline seasonal model can be set as a zero vector, relying on climate indices to drive the construction of seasonal components.

[0115] Step 803: Construct a multivariate decomposition objective function that includes driving alignment terms, so that the residual modes generated by the decomposition maintain a high correlation with the corresponding meteorological driving variables.

[0116] In this step, we define the sequence R applicable to the residual terms. * The objective function for multivariate SVMD decomposition is given by: J = α·J band +β·J drive+γ·J phys J band J is a frequency band sparsity constraint term. drive To drive the alignment item, J phys This is a physical constraint penalty term.

[0117] Band sparsity constraint term J band =∑ k ∑ ω |(ω-ω k ) 2 ·U k (ω)| 2 ;

[0118] Where, ω k U is the center frequency of the k-th mode. k (ω) represents the frequency domain of this mode. This ensures that each mode has a compact bandwidth in the frequency domain.

[0119] Drive Alignment Item J drive =-∑ k ∑ m w k_m ·(Corr(u k ,X m )) 2 ;

[0120] Among them, u k X is a time-domain modal sequence. m w is the m-th driving variable, i.e., the m-th meteorological driving sequence after being filtered through the corresponding frequency band. k_m Corr(u) represents the weighting coefficient. k ,X m The ) represents the correlation coefficient between the two, usually referring to the Pearson correlation coefficient. Minimizing the square of the negative correlation coefficient effectively maximizes the correlation between the mode and the drive. Physical constraint penalty term J phys It contains three sub-items:

[0121] Non-negative punishment P neg =∑ t (max(0,-Q recon (t))) 2 ; where Q recon This represents the total runoff after reconstruction.

[0122] Water balance penalty P bal =|∑ t R(t)-∑ t P eff (t)| / ∑ t P eff (t); where R(t) is the runoff at time t, P eff(t) represents the effective precipitation at time t.

[0123] Monotonicity penalty P mono The value is 1 when the correlation coefficient between the trend term and the cumulative driver is negative, and 0 otherwise.

[0124] The parameters α, β, and γ control the relative importance of the three items, and can be selected as α=1.0, β=0.5, and γ=0.3.

[0125] In some alternative implementations, the driving alignment term can take into account time delay effects, i.e., calculate u k (t) and X m The correlation of (t-τ), where τ is the lag time, can be determined through cross-correlation analysis to identify the optimal lag. Furthermore, the weighting coefficient w... k_m The weighting can be adaptively set based on the frequency characteristics of the modes. For example, high-frequency modes have a larger weighting for precipitation, while low-frequency modes have a larger weighting for temperature. The specific form of the physical constraint terms can also be adjusted according to the characteristics of the watershed. For example, evaporation constraints can be added in arid areas, while temperature threshold constraints can be added in glacial watersheds.

[0126] Example 9 details how to screen effective mode sets through triple constraints of energy density, physical band mapping, and driving correlation, thereby solving the problems of dimensionality curse and noise introduction caused by redundant components generated by traditional mode decomposition.

[0127] Step 901: Calculate the modal energy density integral and energy entropy of each residual modal subsequence, and remove low-energy modes whose energy density integral is lower than a preset energy threshold.

[0128] In this step, the IMFs of all modal subsequences obtained from SVMD decomposition are processed. k Perform energy assessment. Modal energy density E k The calculation formula is E k =∑ t (IMF k (t)) 2 This reflects the total signal energy carried by that mode. To determine an effective energy threshold, the total energy E of all modes is first calculated. total =∑ k E k ; Calculate the energy ratio for each mode k =E k / E total .

[0129] Meanwhile, in order to evaluate the information complexity of the modes, the energy entropy H is calculated. en The specific calculation process is as follows: First, obtain the instantaneous amplitude 'a' of the mode. k (t)=sqrt((IMF k(t)) 2 +(H(IMF k (t))) 2 H represents the Hilbert transform operator.

[0130] Calculate the energy probability distribution p k =(a k 2 ) / (∑a k 2 ).

[0131] Calculate the energy entropy H using the Shannon entropy formula. en =-∑ k (p k ·ln(p k )).

[0132] Set energy threshold E thresh =0.02·E total This means that each effective mode must account for at least 2% of the total energy. When the E of a certain mode... k <E thresh And H en <H thresh When the preset entropy threshold (usually 2.0) is reached, the mode is considered to mainly contain noise components and is therefore removed. This dual-screening mechanism considers both energy magnitude and information complexity, avoiding the monotonic noise components that might be retained by simple energy screening.

[0133] In some alternative implementations, the energy threshold can be determined adaptively. For example, a cumulative curve of modal energy percentage can be plotted, the inflection point of the curve's slope can be identified, and the energy value corresponding to the inflection point can be used as the threshold. Furthermore, the signal-to-noise ratio (SNR) can be introduced as a supplementary screening metric, calculating the statistical difference between each mode and Gaussian white noise, and eliminating modes with an SNR below 3 dB. For specific applications, such as flood forecasting, the energy threshold for high-frequency modes can be appropriately reduced to retain more rapid response information.

[0134] Step 902: Calculate the instantaneous frequency of each mode using Hilbert transform, perform frequency band mapping based on typical time scales of hydrophysical processes, and retain modes that fall within physically meaningful frequency bands.

[0135] In this step, the frequency characteristics of each mode are determined and physically interpreted using Hilbert spectral analysis.

[0136] IMF for residual modal subsequences k The analytic signal z is obtained by performing the Hilbert transform. k (t):

[0137] z k(t)=IMF k (t)+i·H(IMF k (t)); where H represents the Hilbert transform operator and i is the imaginary unit.

[0138] Instantaneous phase φ k (t)=arctan(Im(z k (t)) / Re(z k (t)));Im(...) and Re(...) represent the imaginary and real parts of a complex number, respectively.

[0139] instantaneous frequency f k (t)=(1 / 2π)·dφ k (t) / dt.

[0140] Center frequency f center_k The time average of the instantaneous frequency is taken.

[0141] Step 903: Calculate the cross-correlation coefficients between each mode and key meteorological driving variables, and screen out physical driving modes that are significantly correlated with the driving factors.

[0142] In this step, the correlation between each mode and the external drive is evaluated. For the modes selected through energy and frequency screening, the IMF is... k Calculate the cross-correlation functions between the cross-correlation function and key driving variables such as precipitation series P(t), temperature series T(t), and evaporation ET(t). Cross-correlation function R km The formula for calculating (τ) is:

[0143] ;

[0144] Where τ is the lag time, X m Let σ be the m-th driving variable. IMF_k , σ X_m Let R be the corresponding standard deviation. Find the standard deviation that makes R0... km (τ) The optimal lag with the largest absolute value. opt Record the corresponding maximum correlation coefficient ρ max_km Set the correlation threshold ρ thresh =0.3, and a significance test is performed to calculate the p-value. When |ρ max_km |>ρ thresh Furthermore, when p < 0.05, the mode is considered to have a clear physical relationship with the driving variable. A mode is included in the valid mode set as long as it satisfies the above conditions with any driving variable. This ensures that the retained modes are not artificial products of mathematical decomposition, but rather physical response components with a clear driving source.

[0145] In some alternative implementations, nonlinear correlation measures, such as mutual information (MI) or maximum information coefficient (MIC), can be used to capture the nonlinear dependence between modes and drivers. Furthermore, multivariate correlation indices can be constructed to consider the joint correlation between modes and multiple driving variables. For driving relationships with threshold effects, such as runoff response only occurring when precipitation intensity exceeds a certain threshold, conditional correlation analysis can be used to calculate correlations only during periods when predetermined conditions are met.

[0146] Example 10 details the specific calculation process of the multivariate decomposition objective function, which includes frequency band constraints, drive alignment, and physical penalties, and its iterative solution algorithm based on the alternating direction multiplier method.

[0147] Step 1001: Construct a complete multivariate SVMD objective function and clarify the mathematical expressions and physical meanings of each term.

[0148] In this step, the complete objective function J of multivariate SVMD is defined as:

[0149] J=αJ band +βJ drive +γJ phys +J recon The parameters α, β, and γ control the relative importance of the three terms, respectively.

[0150] Among them, J recon To reconstruct the error term, This ensures the completeness of the decomposition.

[0151] Band sparsity constraint term J band The expanded form can be:

[0152] ;

[0153] Among them U k (ω) represents the mode u k The Fourier transform of (t). This term guarantees the narrowband characteristics of each mode by penalizing the spread of the spectrum.

[0154] Drive Alignment Item J drive The expanded form can be:

[0155] ; ;

[0156] Among them, Corr(u k ,X m ) represents the correlation coefficient between the two, σ(u) k ), σ(X m ) represents their respective standard deviations.

[0157] Physical constraint penalty term J phys The expanded form can be: J phys =λ1·P neg +λ2·P bal +λ3·P mono ;

[0158] The ReLU function is defined as max(0,x).

[0159] .

[0160] The parameters can be set as follows: α=1.0, which controls the sparsity intensity of the frequency band; β=0.5, which controls the driving alignment intensity; γ=0.3, which controls the physical constraint intensity; and λ1=10, λ2=5, and λ3=2, which control the weights of the three physical sub-items respectively.

[0161] In some alternative implementations, a regularization term can be added to the objective function to prevent overfitting, such as an L2 regularization term ∑. k ||u k ||2 2 Or L1 regularization term ∑ k ||u k 1. Weight w in the driving alignment term k_m It can be adaptively set according to modal frequency and drive type, for example, w k_m =exp(-|f k -f m_typical |); where f m_typical The typical response frequency driving m is given. Furthermore, time-varying weights can be introduced, increasing the driving alignment weight of high-frequency modes during flood season and increasing the constraint weight of low-frequency modes during dry season.

[0162] Step 1002: The constrained optimization problem is transformed into an iterative solution of a series of subproblems using the Alternating Direction Multiplier Method (ADMM).

[0163] In this step, the ADMM algorithm framework is used to solve the complex constrained optimization problem described above. First, an auxiliary variable ν is introduced. k And dual variable (Lagrange multiplier) λ k This transforms the original problem into an augmented Lagrangian form. The augmented Lagrangian function L is defined. aug :

[0164] ;

[0165] Where 〈·,·〉 denote the inner product, and ρ is the penalty parameter. ADMM decomposes the global optimization into three alternating sub-steps: the u update step, the ω update step, and the dual variable λ. k Update steps.

[0166] In the u update step, the center frequency ω is fixed. k and dual variable λ k For each mode u k Solution:

[0167] u k_new =arg min u (J recon +β·J drive +A(u k ;ν k ,λ k )); A(u k ;ν k ,λ k (This is to augment the Lagrange term.)

[0168] This subproblem can be solved in the frequency domain using a Wiener filter in a closed-form solution:

[0169] ;

[0170] Among them, U input (ω) is the input spectrum derived from the reconstruction error and other coupling terms; U k_new (ω) is u k_new Fourier transform.

[0171] In the ω update step, the mode u is fixed. k Update center frequency:

[0172]

[0173] Dual variable λ k The update uses standard gradient ascent: λ new =λ old +ρ·(u k_new -ν k_new The iterative process continues until the convergence criterion is met: the u values ​​of two consecutive iterations are equal. k The change norm is less than ε1 and the constraint violation is less than ε2.

[0174] ;

[0175] Where ε1 and ε2 are each set to 10. -6 and 10 -7 .

[0176] In some alternative implementations, an adaptive penalty parameter ρ can be used, dynamically adjusted based on the relative magnitudes of the original and dual residuals, to accelerate convergence. Furthermore, a momentum term or Nesterov acceleration can be introduced to improve the algorithm's convergence speed. For large-scale problems, a stochastic or mini-batch version of ADMM can be used, updating only a subset of modes in each iteration to reduce computational complexity.

[0177] Step 1003: Apply causal constraints to ensure that the decomposition process uses only historical data.

[0178] In this step, to meet the causality requirements of real-time forecasting, strict time constraints are imposed on the Alternating Direction Multiplier Method (ADMM) solution process. In the frequency domain filtering operation, the standard bilateral filter is replaced with a single-sided filter. Specifically, a causal filtering kernel h is defined. causal (τ), satisfying:

[0179]

[0180] In calculating mode u k When (t), the convolution operation is modified as follows:

[0181] It is guaranteed that the output at time t depends only on the input at time t and before.

[0182] For the Hilbert transform, a causal Hilbert transform is adopted, based on the Kramers-Kronig relation (KK relation for short):

[0183] ; where PV represents the Cauchy principal value integral, used to deal with the singularity of the integrand at s=t.

[0184] In the ADMM iterative update, for the modal value u at time t k The calculation of (t) involves only historical data points with indices less than or equal to t. While this causal constraint slightly reduces the mathematical optimality of the decomposition, it ensures that the algorithm can be applied to online real-time scenarios, requiring only incremental updates when new data arrives, without needing to recalculate the entire historical sequence.

[0185] In some alternative implementations, a sliding window strategy can be used to achieve quasi-causal decomposition. A buffer window of length L' (e.g., 120 time steps) is set. At each time t, local SVMD decomposition is performed using data in the interval [t-L'+1,t]. The decomposition result at the end of the window is taken as the output at time t. Furthermore, a forgetting factor can be introduced to apply exponentially decaying weights to historical data, making recent data have a greater impact on the current decomposition and improving adaptability to non-stationary changes.

[0186] Example 11: A detailed explanation of the decomposition-side closed-loop feedback mechanism based on the dual indicators of information retention rate and physical consistency, including the quantitative calculation method of the indicators and the adaptive parameter adjustment strategy.

[0187] Step 1101: Construct a comprehensive index for information retention rate to quantitatively evaluate the information integrity after modality screening.

[0188] In this step, we define the comprehensive index I of information retention rate. info Used to evaluate the degree to which the effective mode set retains information from the original residual signal. This index is composed of three weighted components: I info =w E ·ER+w ρ ·ρ time +w P • PSI; weights set to w E =0.4、w ρ =0.35、w P =0.25.

[0189] The energy retention rate ER is calculated as the reconstructed residual sequence R. recon Total energy and original residual (residual term sequence) R * The ratio of total energy: ER = (∑ t (R recon (t)) 2 ) / (∑ t (R * (t)) 2 ).

[0190] Time-domain correlation coefficient ρ time ρ is calculated using Pearson correlation. time =Cov(R recon ,R * ) / (σ(R recon )·σ(R * )), where Cov represents covariance and σ represents standard deviation.

[0191] Spectral similarity index (PSI) is calculated by comparing power spectral densities. First, a fast Fourier transform is performed on the two sequences to obtain the spectra F. recon (f) and F * (f) Calculate the power spectrum S recon (f)=F recon (f)·conj(F recon (f)), S * (f) Similar. Spectral similarity PSI = 1 - D spec D spec Let be the normalized spectral distance, where:

[0192] Dspec =(∫|sqrt(S recon (f))-sqrt(S * (f))|df) / (∫sqrt(S * (f))df). To capture local information loss, the local information retention rate I of the sliding window is calculated. info_local (t), with a window length set to 60 time steps, the final comprehensive information retention rate index I info Take the weighted average of the global value and the local minimum:

[0193] I info_final =0.7·I info_global +0.3·min(I info_local (t)).

[0194] In some alternative implementations, higher-order statistics can be introduced to assess the degree of information preservation, such as the preservation of skewness and kurtosis. Specifically, statistical similarity S is defined. stat :S stat =(S skew +S skurt ) / 2;

[0195] ; ;

[0196] Among them, S skew For skewness similarity, S skurt For kurtosis similarity, skew(.) represents the skewness function and skurt(.) represents the kurtosis function.

[0197] Furthermore, for residual signals with obvious periodicity, the retention rate of the main periodic components can be calculated separately to ensure that key periodic information is not lost. The information retention rate can also be calculated by frequency band, assigning different weights to different frequency bands. For example, higher weights can be given to the low-to-medium frequency bands that contain the main hydrological information.

[0198] Step 1102: Construct physical consistency indicators to assess the hydrophysical rationality of reconstructing runoff.

[0199] In this step, the physical consistency index I is defined. phys Used for quantification and reconstruction of total runoff Q recon The physical plausibility of Q, among which recon =T * +S * +R recon This metric consists of three physical constraint scores:

[0200] I phys =0.4·S neg +0.35·Sbal +0.25·S mono .

[0201] Non-negative score S neg Assess the severity of negative flow:

[0202] S neg =exp(-10·r neg ), where r neg For negative proportions, r neg =count(Q recon <0) / length(Q recon ), length(Q recon S is the signal length. When there are no negative values, S... neg =1, the more negative the value, the lower the score index.

[0203] Water balance score S bal Assessing long-term water conservation:

[0204] S bal =exp(-5·e bal ), where e bal The relative balance error is calculated as: e bal =|∑(Q recon )-∑(Q obs )| / ∑(Q obs ).

[0205] Monotonicity score S mono Assess the consistency between trends and drivers:

[0206] Calculate the trend term sequence T * With long-term driving indicator X agg correlation coefficient r trend When r trend When S > 0, mono =1 indicates that it conforms to the physical expectation that the trend will increase with the increase of cumulative precipitation; when r trend When <0, S mono =0.5-0.5·|r trend The more negative the correlation, the lower the score.

[0207] In addition, an extreme value rationality check is introduced to calculate the reconstructed total runoff Q. recon The maximum value and the historical maximum value Q max_obs The ratio r peak r peak =max(Q recon ) / Q max_obs When r peak Additional penalty is applied when the value is greater than 2, I phys '=I phys·exp(-(r peak -2)), I phys 'This is a metric for physical consistency after applying additional penalties.'

[0208] In some alternative implementations, the physical consistency assessment can incorporate seasonality testing. The monthly average pattern of reconstructed runoff is calculated. recon Compared with historical monthly average patterns obs The similarity is assessed using cosine similarity or dynamic time warped distance (DTW). Furthermore, it can be used to examine whether the autocorrelation function (ACF) of the reconstructed sequence maintains reasonable decay characteristics, preventing non-physical long-range correlations. For specific watersheds, physical constraints can be added, such as temperature threshold constraints for glacial watersheds and baseflow stability constraints for karst watersheds.

[0209] Step 1103: Implement closed-loop parameter adjustment driven by dual indicators until the decomposition results simultaneously meet information and physical requirements.

[0210] In this step, a complete closed-loop optimization process is established. A preset information retention rate threshold T is set. info =0.9 and physical consistency threshold T phys =0.85. The pseudocode for the closed-loop adjustment algorithm is as follows:

[0211] Initialize the number of iterations iter=0, and the maximum number of iterations max iter =10.

[0212] When iter < max iter And (I info <T info or I phys <T phys The loop executes when I... Within the loop, SVMD decomposition is first performed to obtain the mode set, mode filtering is performed to obtain the effective set, the reconstructed signal is calculated, and two metrics are evaluated. If I... info <T info This indicates excessive information loss, necessitating adjustments to prioritize information retention.

[0213] Update the energy screening threshold, E thresh_new =E thresh • 0.8, lower the screening criteria; update the upper limit of modality count, K max_new =K max +1, allowing more modes to be retained; update the frequency band constraint weight α. new =α·0.9, relax the requirement for sparse frequency bands.

[0214] If I phys <T physThis indicates insufficient physical rationality, necessitating a physical constraint-oriented adjustment: updating the physical penalty weight γ. new =γ·1.3, strengthen physical constraints; update driving alignment weight β new =β·1.2, enhancing the correlation with the driver; adjusting the modal merging threshold to merge modes with similar frequencies to reduce spurious components. After each parameter adjustment, re-execute the decomposition and evaluation, iter new =iter+1. The loop exits when both metrics meet the target or the maximum number of iterations is reached. This dual-metric closed-loop mechanism ensures the quality of the decomposition results in both mathematical and physical dimensions.

[0215] In some alternative implementations, the parameter tuning strategy can be optimized using reinforcement learning. The information retention rate is integrated into the overall index I. info Physical consistency index I phys The combination forms a reward function R'=I info +I phys -P ip , where P ip The penalty term includes the number of iterations and the magnitude of parameter change. A policy network is trained, taking the current metric value and parameter state as input, and outputting parameter adjustment actions. Furthermore, a historical record of parameter adjustments can be maintained, and a Gaussian process model can be constructed based on historical data using Bayesian optimization methods to predict metric values ​​under different parameter combinations, guiding the parameter search direction.

[0216] This invention achieves deep integration of signal processing and physical mechanisms by constructing a multi-scale co-decomposition framework based on driving forces and physical constraints, thus solving the problem of unclear physical meaning in traditional decomposition methods. By designing a component-level physical constraint integrated prediction model and a multi-objective loss function, it overcomes the defect of pure data-driven models that are prone to producing non-physical prediction results. By introducing a time-varying adaptive Bayesian optimization mechanism, it solves the problem of model adaptability in non-stationary climatic environments.

[0217] Specifically, the present invention employs a multi-scale co-decomposition method involving both driving forces and physical constraints. By introducing long-term climate factors into the STL decomposition and explicitly adding driving alignment terms and physical penalty terms to the objective function of the SVMD decomposition, the decomposed modes are forced to maintain a high correlation with meteorological driving forces such as precipitation and temperature and conform to hydrological laws, thus solving the problem of the signal decomposition process being disconnected from the physical mechanism.

[0218] This invention addresses the real-time application challenges caused by non-causal decomposition by imposing a one-sided filter kernel constraint on the decomposition algorithm, thus shielding it from future data involvement. To address the lack of physical constraints and insufficient adaptability of the model, this invention constructs a component-level physically constrained prediction model and introduces a time-varying adaptive Bayesian optimization mechanism. By incorporating water balance, non-negativity, and extremum weighting terms into the loss function, and dynamically updating hyperparameters using a sliding window detection mechanism, it effectively solves the problems of prediction results violating physical realities and concept drift caused by the model's inability to adapt to climate change. While maintaining overall accuracy, it improves the ability to capture extreme events and the system's robustness.

[0219] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for runoff prediction based on STL-SVMD decomposition and adaptive integration, characterized in that, The method comprises the following steps: collecting runoff time series and multi-source meteorological data of a target basin; performing multi-scale decomposition on the runoff time series by using a generalized multi-scale decomposition method with meteorological driving elements introduced as constraint conditions, separating a trend item sequence, a seasonal item sequence and a residual item sequence responding to long-term climate change, and further decomposing the residual item sequence to extract residual item modal sub-sequences; constructing a candidate feature set based on the trend item sequence, the seasonal item sequence, the residual item modal sub-sequences and the multi-source meteorological data, and screening to obtain a key input factor set; inputting the key input factor set into a preset prediction model, and after parameter optimization, respectively predicting a trend item prediction value, a seasonal item prediction value and a residual item modal prediction value; reconstructing the trend item prediction value, the seasonal item prediction value and the residual item modal prediction value to generate a runoff prediction result of the target basin; wherein, performing multi-scale decomposition on the runoff time series further comprises: calculating long-term scale cumulative or moving average meteorological elements from the multi-source meteorological data, and using the cumulative or moving average meteorological elements and time index as independent variables to construct the trend item sequence responding to long-term climate change; extracting an interannual climate index reflecting large-scale climate background from the multi-source meteorological data as a modulation factor, decomposing seasonal changes into a baseline seasonal pattern and a dynamic offset controlled by the modulation factor, and constructing a seasonal item sequence evolving with the climate background; further decomposing the residual item sequence, comprising: selecting key meteorological driving variables related to hydrological processes from the multi-source meteorological data, concatenating the residual item sequence and the key meteorological driving variables along the channel dimension to form a multivariate input vector; constructing a multivariate decomposition objective function containing a driving alignment item based on the multivariate input vector, configuring the driving alignment item to maximize the correlation between the modal components generated by the decomposition and the corresponding key meteorological driving variables, and guiding the modal components to align with the physical driving process; by iteratively solving the multivariate decomposition objective function, the multivariate input vector is decomposed into multiple pairs of modal components, and the residual item modal sub-sequences are extracted from the corresponding residual item sequences; the multivariate decomposition objective function further contains a physical penalty term, which is configured to impose a numerical penalty on the decomposition results that violate hydrological physical laws, and the hydrological physical laws at least include: the non-negativity constraint of the flow of the reconstructed runoff, the long-term water balance constraint after the decomposition components are superimposed, or the monotonicity constraint of the change trend between the trend item sequence and the key meteorological driving variables.

2. The method of claim 1, wherein, performing multi-scale decomposition on the runoff time series, comprising: processing the runoff time series by using the STL seasonal-trend decomposition algorithm to decouple the trend item sequence representing long-term changes, the seasonal item sequence representing periodic fluctuations, and the residual item sequence containing random disturbances; inputting the residual item sequence into a successive variational modal decomposition model, constructing a variational constraint problem and iteratively solving it by using the alternating direction multiplier method to successively extract residual item modal sub-sequences with different center frequencies.

3. The method of claim 2, wherein, Before constructing the candidate feature set, further comprising: calculating the modal energy density integral and the Hilbert spectrum center frequency of each residual item modal sub-sequence; Subsequences with modal energy density integral below a preset energy threshold or Hilbert spectrum center frequency deviating from a preset physical frequency band are removed, and the remaining high-energy effective modal subsequences are retained for constructing a candidate feature set.

4. The method of claim 3, wherein, The key input factor set is screened, including: Pearson correlation coefficients of each feature in the candidate feature set and the target runoff sequence are calculated, and the correlation features with absolute correlation coefficients below a preset correlation threshold are removed; The retained features are input into a random forest regression model for training, the feature importance of each feature in the model decision tree splitting is quantitatively evaluated, and the top features are selected according to the importance ranking results; The features that meet both the Pearson correlation coefficient requirement and the feature importance requirement are selected from the top features to form a key input factor set that takes into account linear correlation and nonlinear contribution.

5. The method of claim 1, wherein, The trend item prediction value, the seasonal item prediction value and the residual item modal prediction value are predicted respectively, including: For the trend item sequence, the seasonal item sequence and each effective modal subsequence, a prediction sub-model containing an XGBoost module and an LSTM module is independently constructed; The XGBoost module is used to mine the nonlinear regression features between meteorological factors and runoff components from the key input factor set, and the LSTM module is used to capture the long-short time sequence dependence features in the key input factor set; The trend item prediction value, the seasonal item prediction value and the residual item modal prediction value of the corresponding component are output respectively by each prediction sub-model.

6. The method of claim 5, wherein, Parameter optimization, including: The root mean square error or Nash-Sutcliffe efficiency coefficient is selected as a single evaluation index, and a target function for evaluating the generalization error of the prediction sub-model is constructed based on the pre-collected validation set data; The search space of the hyperparameters of the XGBoost module and the LSTM module in the prediction sub-model is defined, and the Bayesian optimization algorithm is used to perform probability agent model construction and function calculation in the hyperparameter search space to minimize the target function for iterative optimization; The optimal hyperparameter combination that minimizes the target function is obtained, and it is loaded into the corresponding prediction sub-model to lock the model state for performing the final runoff prediction.

7. The method of claim 1, wherein, Iterative solution of the multivariate decomposition target function, including: During the solution process, causal constraints are applied, and a one-sided support kernel or causal transformation operator with one-sided support characteristics is used to process the multivariate input vector; When calculating the modal component at any current time, only the multivariate input vector data at the current time and historical times are used, and the participation of future time data is shielded, to generate a residual item modal subsequence that meets the real-time prediction requirements.

8. The method of claim 7, wherein, After extracting the residual item modal subsequence, it further includes: The filtered residual item modal subsequences are superimposed to generate a reconstructed residual sequence, and an information retention rate comprehensive index is calculated based on the reconstructed residual sequence and the residual item sequence before decomposition, which comprehensively represents the energy retention rate, time domain correlation degree and frequency spectrum similarity; The index is compared with a preset retention rate threshold, if the preset retention rate threshold is not met, the energy threshold or the physical frequency band range of modal screening is dynamically adjusted, and the screening and reconstruction steps are re-executed until the information retention rate comprehensive index meets the requirements.

Images