Runoff prediction methods, devices, and media that integrate physical priors and machine learning

Abstract

This invention discloses a method, apparatus, and medium for runoff prediction that integrates physical priors and machine learning. The runoff prediction method uses the basic watershed unit as the smallest spatial unit, first running multiple physical / conceptual hydrological models in parallel, extracting physical prior factors from the output of each model to construct a physical prior factor library; systematically feature-engineering the original time-series meteorological, remote sensing, topographic, and prior factors; training an integrated interpretable machine learning model using joint features, and calculating PBIAS reflecting water balance after training to verify physical consistency; calculating the comprehensive importance score of factors using an integrated interpretation framework and extracting response thresholds; subsequently, partitioning the basic watershed unit based on contribution vector clustering to achieve adaptive partitioning and runoff prediction based on dominant factors. Implementation results show that this method can improve the accuracy of daily runoff simulation and enhance the physical interpretability and engineering operability of the model under various climatic and topographic conditions.

Description

Technical Field

[0001] This invention belongs to the field of hydrology and water resources engineering, and involves watershed runoff prediction, interpretable machine learning and regional modeling technology. In particular, it relates to a method and system that uses the output of multiple physical / conceptual hydrological models as physical priors and meteorological and underlying surface factors as inputs into a machine learning model, and quantitatively extracts factor contributions and response thresholds based on the joint analysis of SHAP, PDP, ICE and PFI (HyPS-IF). Subsequently, it partitions the watershed basic units (HRUs) based on contribution vector clustering (K-means) to achieve adaptive selection and dynamic simulation according to the dominant factors. Background Technology

[0002] Traditional process-oriented hydrological models (such as SWAT and TOPMODEL) possess good physical interpretability, but suffer from problems such as complex parameter calibration, high data requirements, and difficulties in cross-domain transfer. Pure data-driven methods (including deep learning and traditional machine learning) can achieve high fitting accuracy in many cases, but often neglect physical priors, lack interpretability, and are prone to physically unreasonable predictions under small sample or extrapolation conditions. Existing research attempts to integrate physical models and data-driven models, but generally faces challenges such as: how to systematically extract and screen "verifiable" physical prior factors, how to introduce physical consistency checks while ensuring training efficiency, and how to transform interpretable results into engineeringable partitioning and model adaptation rules.

[0003] Meanwhile, existing runoff prediction schemes mostly rely on a single physical model or a single data-driven model, which makes it difficult to adaptively match the differences in runoff generation and confluence mechanisms under different climate zones, topography and underlying surface conditions. In complex watersheds or areas with large climate gradient changes, problems such as decreased simulation accuracy, insufficient physical consistency and limited regional applicability are likely to occur, making it difficult to meet the actual needs of engineering-based, long-sequence, multi-regional collaborative runoff prediction. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a method, device and storage medium for runoff prediction that integrates physical priors and machine learning to improve prediction accuracy.

[0005] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0006] This invention first provides a runoff prediction method that integrates physical priors and machine learning, comprising the following steps:

[0007] S1. Divide the target watershed into basic watershed units and collect multi-source data for each basic watershed unit;

[0008] S2. Based on the collected multi-source data, run various types of hydrological models for each basic watershed unit to extract candidate physical prior factors; among which, the different types of hydrological models include physical process hydrological models, semi-physical and semi-empirical hydrological models, and conceptual hydrological models;

[0009] S3. Perform single-factor perturbation sensitivity analysis on the extracted candidate physical prior factors, and select the core physical prior factors from the candidate physical prior factors to construct a joint training feature set.

[0010] S4. Construct a Stacking ensemble model; using the joint training feature set as input, train the Stacking ensemble model for daily runoff prediction, and calculate the percentage deviation index PBIAS after training to verify the physical consistency of the water balance in the prediction results.

[0011] S5. The trained Stacking ensemble model is subjected to interpretability analysis using an ensemble interpretation framework. The comprehensive importance score of each feature factor in the joint training feature set is calculated. The response threshold of the feature factors in the joint training feature set whose comprehensive importance score is higher than a preset threshold is extracted. Based on the interpretability analysis results of the comprehensive importance score, response threshold and single factor contribution intensity, an HRU-level contribution vector is constructed for each watershed basic unit.

[0012] S6. Based on the HRU-level contribution vectors of each watershed basic unit, the K-means method is used to perform clustering partitioning of the target watershed in units of watershed basic units;

[0013] S7. For each partition after clustering, evaluate within the candidate model pool consisting of the various hydrological models of different categories run in step S2 and the Stacking ensemble model trained in step S4, and determine the preferred model set for that partition; output the daily runoff simulation results based on the preferred model set determined for each partition.

[0014] The runoff prediction method of this invention, at the HRU granularity, firstly divides the target watershed into basic watershed units and collects multi-source data; then, it runs candidate physical / conceptual hydrological models for each basic watershed unit to extract candidate physical prior factors; next, it performs single-factor perturbation sensitivity analysis on the candidate physical prior factors, screens core physical prior factors, and constructs a joint training feature set; it uses the joint training feature set as input to construct and train a Stacking ensemble model, and verifies the physical consistency of water balance using PBIAS after training; subsequently, it performs interpretability analysis on the trained Stacking ensemble model, calculates the comprehensive importance score of each feature factor in the joint training feature set, and extracts the response threshold; further, it uses the K-means method to partition the runoff response based on the HRU-level contribution vector of each basic watershed unit; finally, for each partition, it determines the preferred model set from the candidate model pool and outputs the daily-scale runoff simulation results; in some embodiments, it can also implement re-screening, retraining, and iterative optimization of partitioning / model adaptation according to preset performance trigger conditions during the online operation phase.

[0015] The present invention also provides a flow prediction device that integrates physical priors and machine learning, including a processor and a memory; the memory stores a program or instructions, which are loaded and executed by the processor to implement the steps of the runoff prediction method provided above.

[0016] The present invention also provides a computer-readable storage medium storing a program or instructions that, when executed by a processor, implement the steps of the runoff prediction method described above.

[0017] Compared with the prior art, the beneficial technical effects of the present invention using the above technical solution are as follows:

[0018] The runoff prediction method integrating physical priors and machine learning provided by this invention offers the following advantages over existing technologies: This invention uses the watershed basic unit (HRU) as the smallest spatial granularity, jointly running multiple physical / conceptual hydrological models (such as GXAJ, SAC-SMA, VIC, SWAT, and TOPMODEL) in parallel. It systematically extracts and filters core physical prior factors from the outputs of each model through single-factor perturbation sensitivity analysis, thus introducing verifiable physical information as machine learning features into the modeling process. Based on systematic temporal feature engineering and Stacking-type interpretable ensemble learning, and by introducing water balance-based penalty terms or multi-task learning during training or post-processing to ensure physical consistency, it significantly improves the accuracy and physical reliability of daily-scale runoff simulation. Simultaneously, this invention calculates the comprehensive importance of factors and extracts response thresholds using multi-dimensional interpretation tools (SHAP, PFI, PDPs / ICE), and uses the K-means method to partition features based on the feature importance scores of different HRUs. This method has significant advantages in improving the accuracy of daily runoff simulation, enhancing model interpretability and cross-basin mobility, and can provide more reliable scientific support for water resource management, flood early warning and engineering decision-making. Attached Figure Description

[0019] Figure 1 shows the overall flowchart of the method of the present invention.

[0020] Figure 2 Map showing the HRU delineation results for the Shuimingya watershed.

[0021] Figure 3 Sensitivity ranking diagram of factors in the Shuimingya watershed.

[0022] Figure 4 Importance ranking chart of SHAP in the Shuimingya watershed.

[0023] Figure 5 PDP response curves for the six core factors.

[0024] Figure 6 Clustering results of the Shuimingya watershed. Detailed Implementation

[0025] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0026] In this invention, various aspects of the invention are described with reference to the accompanying drawings, in which numerous illustrative embodiments are shown. Embodiments of the invention are not limited to those depicted in the drawings. It should be understood that the invention is implemented through any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0027] Example 1

[0028] refer to Figure 1 This embodiment provides a runoff prediction method based on the fusion of physical priors and machine learning, including the following steps:

[0029] S1. Division of Basic Watershed Units and Data Preparation

[0030] Division principle: The target watershed is divided into several HRUs based on topography, river network, and management units, with each HRU area controlled between 10 and 50 km². 2 To ensure both hydrological uniformity and facilitate parallel computing, the optimal area for each HRU is 1–20 km². 2 For each HRU, a spatial index and attribute table (including HRU_ID, area, average slope, dominant soil type, main vegetation type, and adjacency relationship) are established. The granularity of HRU division can be adjusted according to engineering requirements. Smaller HRUs are beneficial for localization accuracy but increase computational load. This is accomplished using ArcMap software.

[0031] To verify the applicability and transferability of the method of this invention in different climate zones, this study covers four climate zones in the Chinese Climate Zoning Standard (GB / T 20107–2006): arid, semi-arid, semi-humid, and humid. Within each climate zone, representative small watersheds with topography, area, and land use were selected as validation samples (see Table 1). The selected samples form gradients in terms of multi-year average precipitation, slope, and intensity of human activities, which is beneficial for evaluating the ability of the method of this invention to characterize different runoff generation and confluence mechanisms and the cross-regional parameter / model transfer performance. Table 1 lists the basic attributes of the nine experimental watersheds (annual average precipitation, watershed area, average slope, and average elevation). Figure 2 The HRU delineation results for the Shuimingya watershed are shown as an example.

[0032] Table 1 Basic attributes of typical small watersheds in the study area

[0033]

[0034] Multi-source data collection: The following data are collected for each HRU: (1) Daily rainfall P, runoff Q, and evaporation data PET from hydrological stations, from which the previous day's rainfall P is derived. t-1 (2) Remote sensing meteorological data: daily maximum temperature T max Daily minimum temperature T min Daily average temperature T avg Dew point temperature (Td), wind speed (Wind), solar radiation (Solar), direct radiation (Thermal) (derived from ERA5-Land), relative humidity (RH) = Td / T avg Net radiation Rn = Solar-Thermal; Drought index AI is calculated by the Penman formula: AI = Potential evapotranspiration PET / Concurrent precipitation P; (3) Normalized Difference Vegetation Index (NDVI) (NOAA CDR); (4) DEM (resolution 30 m), and then the average slope Slope is calculated using ArcMap software. avg (5) Soil type and land use data.

[0035] The time-varying meteorological and hydrological continuous variables in the multi-source data are uniformly processed into daily-scale time series, and missing value imputation and outlier detection are performed; the static topographic continuous features and categorical features in the multi-source data are subjected to spatial consistency checks, missing item completion and coding preprocessing.

[0036] S2. Run the multiphysics model and extract candidate physics prior factors.

[0037] S21. Model Configuration and Parameterization:

[0038] Multiple hydrological models were configured in the Shuimingya watershed of the study area, including the physical process hydrological model SWAT, the semi-physical and semi-empirical hydrological model TOPMODEL, and the conceptual hydrological models GXAJ, SAC-SMA, and VIC. For each HRU, based on topography, soil, land use, meteorological conditions, and existing calibration experience, parameter values ​​were assigned, regionalized migrations or parameter calibrations were performed on each model, and the configuration of model input, parameter sets and operation schemes was completed.

[0039] S22. Model Running and Candidate Factor Extraction:

[0040] For each HRU, the selected hydrological model is run within a time frame (e.g., the Shuimingya watershed, 2003-2022). During the model run, hydrological process variables are extracted from the output of each model as candidate physical prior factors. For example, in the GXAJ model, upper soil moisture content WU, lower soil moisture content WL, and simulated runoff Q are extracted. XAJThe upper tension water (UZTWC), lower tension water (LZTWC), and basement flow (BASE) were extracted from the SAC-SMA model; the soil moisture (SM1, SM2, and SM3) of each layer were extracted from the VIC model; the soil water content (SW) was extracted from the SWAT model; and the soil saturation (S(t)) was extracted from the TOPMODEL model.

[0041] The extracted hydrological process variables are all output on a daily scale and used to construct a physical prior factor library.

[0042] S23. Characterization of Time-Varying Factors: Constructing Time-Window Statistical Features: Multi-scale sliding window statistical features are constructed for continuous input variables that have been uniformly processed into daily time series. Continuous input variables include time-varying meteorological and hydrological continuous variables from multi-source data, as well as candidate physical prior factors extracted in S22. Continuous input variables include precipitation. Potential evaporation Daily maximum temperature Daily minimum temperature Daily average temperature Dew point temperature Wind speed Solar radiation thermal radiation relative humidity Net radiation , and the physical prior factors extracted from S22; window length The preferred dates are 7 or 30 days.

[0043] For any continuous input variable Calculate the past separately The cumulative value, mean, standard deviation, maximum value, minimum value, and linear slope within a daily window are expressed as follows:

[0044]

[0045] in, For the past The cumulative value of the variable is continuously entered within the daily window; for The mean of variables is continuously entered within the daily window; for The standard deviation of the variable is continuously entered within the daily window; for The maximum value of continuously entered variables within the daily window; for The minimum value of a variable continuously entered within the daily window; for The linear slope of continuously input variables within the day window; , This represents the mean value within the window.

[0046] Static terrain features such as DEM, slope, and TWI are not included in the sliding window statistics, but are directly input into the Stacking ensemble model in S4 as static continuous features in the joint training feature set.

[0047] Constructing lagged cumulative precipitation characteristics: To characterize the lagged impact of previous precipitation on the current runoff response, lagged cumulative characteristics are constructed for precipitation variables:

[0048]

[0049] in, Indicates the first Daily rainfall, Indicates as of the date Days passed Daily cumulative rainfall.

[0050] S3. Single-factor perturbation sensitivity analysis

[0051] S31. Construct a set of candidate driving factors:

[0052] To identify the core physical prior factors that significantly contribute to runoff prediction, candidate physical prior factors extracted from the outputs of various hydrological models in S22 were used as the objects of single-factor perturbation sensitivity analysis. These candidate physical prior factors include, but are not limited to: upper soil moisture content WU, lower soil moisture content WL, and simulated runoff Q from the GXAJ model. XAJ ; Upper tension water UZTWC, lower tension water LZTWC and basement flow BASE in the SAC-SMA model; Soil moisture SM1, SM2 and SM3 in each layer in the VIC model; Soil water content SW in the SWAT model; Soil saturation S(t) in the TOPMODEL model.

[0053] S32. Perturbation Application and Sensitivity Calculation: For each candidate driver in the candidate driver set... A single-factor perturbation of ±δ% is applied, where δ is 5%–20%, preferably 10%; while keeping the remaining inputs, parameters, and operating conditions of the corresponding hydrological model from which the candidate physical prior factor is extracted unchanged; the baseline NSE values ​​under the baseline operating conditions are then obtained. Candidate physical prior factors NSE value after applying a positive perturbation and the NSE value after applying a negative perturbation. And calculate the relative change using the following formula. :

[0054]

[0055] In the formula, Candidate driving factors Caused by disturbance NSE The relative rate of change of the (Nash-Sutcliffe efficiency coefficient) To correspond to the NSE value of the hydrological model under baseline operating conditions, Candidate driving factors The model after applying a positive perturbation NSE , Candidate driving factors The model after applying a negative perturbation NSE ;

[0056] S33. Judgment Criterion: If Then the candidate physical prior factors The core physical prior factors were identified as the core physical prior factors of the HRU and included in its joint training feature set. These core physical prior factors were used in subsequent S34 to construct statistical and ratio-based derived features, and in the training sample construction of S4, they were used together with meteorological time-series continuous variables, static topographic continuous features, and category coding features as inputs to the Stacking ensemble model. The core physical prior factors of the Shuimingya watershed are as follows: Figure 3 As shown.

[0057] S34. Constructing Physical Prior Derived Features: To enhance the temporal representation ability and physical correlation expression ability of candidate physical prior factors, further derived features are constructed for the candidate physical prior factors extracted in S22. The derived features include statistical derived features and ratio derived features.

[0058] Constructing physical prior derived features: To enhance the ability to represent time series and express physical correlations, statistical and ratio-based derived features are constructed only for core physical prior factors determined by single-factor perturbation sensitivity analysis. Candidate physical prior factors that are not selected as core factors do not participate in the construction of derived features.

[0059] Let the physical prior factors be... Constructing statistical derived features: past Daily moving average Coefficient of variation and quantile difference .

[0060] past Daily moving average Coefficient of variation and quantile difference They are defined as follows:

[0061]

[0062] in, For the past Standard deviation within the daily window, and They are the 75th percentile and the 25th percentile, respectively. To prevent extremely small positive numbers with a denominator of zero, it is preferable to use 10. -6 Throughout the text, when dealing with denominator stabilization, the notation is consistently used. It represents a very small positive number.

[0063] For two physical prior factors that have a physical relationship and Further construct ratio-based derived features:

[0064]

[0065] in, Physical prior factors and Ratio-type derived features;

[0066] Preferably, the ratio-based derived features include at least the following factor pairs with clear physical relationships: the ratio of upper soil moisture content WU to lower soil moisture content WL in the GXAJ model (WU / WL); the ratio of upper tension water UZTWC to lower tension water LZTWC in the SAC-SMA model (UZTWC / LZTWC); and the soil moisture ratios of adjacent soil layers SM1 / SM2 and SM2 / SM3 in the VIC model.

[0067] Optionally, other ratio-based derived features with clear physical correlations can be constructed based on the availability and physical significance of candidate physical prior factors, including but not limited to SM1 / SM3, subsurface baseflow BASE, and simulated runoff Q. XAJ The ratio of the shallow / deep water storage distribution relationship, as well as other ratios that characterize the shallow / deep water storage distribution relationship or the fast / slow flow distribution relationship.

[0068] Feature encoding and normalization are performed: one-hot encoding is applied to categorical features, and Z-score standardization is applied to continuous features.

[0069] Among them, the categorical features include discrete, non-ordered variables such as soil type, land use / cover type, and dominant vegetation type; if the first Each category feature contains If there are several categories, then convert them to... indivual Indicator variable:

[0070] in, For the first The category features in the th _____ One-hot encoding results for each category For this sample in the 1st The values ​​of each category feature This represents the total number of categories for this feature type. Continuous features include continuous meteorological time-series variables, static topographic continuous features, candidate physical prior factors, and statistical and ratio-based derived features of candidate physical prior factors. Continuous features are standardized as follows:

[0071]

[0072] in, These are the original eigenvalues. and These are the original eigenvalues. The mean and standard deviation on the training set, These are the standardized feature values.

[0073] S4 builds and trains the Stacking ensemble model:

[0074] The Stacking ensemble model comprises base learners and meta-learners. Base learners include two or more of XGBoost, LightGBM, RandomForest, and CatBoost, while meta-learners are linear regression or Lasso. The Stacking ensemble model is trained using a joint training feature set as input and diurnal runoff on the target date as output. After training, the percentage bias index (PBIAS) is calculated to verify the physical consistency of the predicted water balance. This includes:

[0075] S41. Training Sample Construction: A sliding window approach is used to construct the sample set. The input is the joint training feature vector within the past window, and the output is the daily runoff on the target date. The joint training feature vector includes: time-varying meteorological and hydrological continuous variables from multi-source data and their sliding window statistical features, precipitation lag cumulative features, static topographic continuous features, category coding features, and core physical prior factors retained after screening in S3, along with their statistical and ratio-based derived features.

[0076] S42. Base Learner Training: Train XGBoost, CatBoost, LightGBM and Random Forest models respectively. Use 5-fold cross-validation to obtain the optimal hyperparameters for each base learner (XGBoost, CatBoost, LightGBM and Random Forest models). Use early stopping to prevent overfitting.

[0077] The optimal parameters for each model are as follows:

[0078] XGBoost:

[0079] Preferred parameters include:

[0080] The maximum tree depth (max_depth) is 6-10, the learning rate (learning_rate) is 0.01-0.1, and the number of base learners (n_estimators) is 200-1000. In this embodiment, the maximum tree depth is 6, the learning rate is 0.03, and the number of base learners is 400. Furthermore, the following settings are used: subsample ratio (subsample) is 0.9, feature sampling ratio (colsample_bytree) is 0.8, L2 regularization coefficient (reg_lambda) is 1.0, L1 regularization coefficient (reg_alpha) is 0.0, the objective function is squared error regression (reg:squarederror), the random seed (random_state) is RANDOM_STATE, and the number of parallel threads (n_jobs) is -1.

[0081] LightGBM:

[0082] Preferred parameters include:

[0083] The number of leaf nodes (num_leaves) is 31-127, the learning rate (learning_rate) is 0.01-0.1, and the number of base learners (n_estimators) is 200-1000.

[0084] In this embodiment, the number of leaf nodes is 31, the learning rate is 0.03, and the number of base learners is 500. Furthermore, the following settings are made: the maximum tree depth (max_depth) is 6, the minimum number of leaf node samples (min_child_samples) is 30, the sample sampling ratio (subsample) is 0.8, the feature sampling ratio (colsample_bytree) is 0.8, the random seed (random_state) is RANDOM_STATE, and the log output level (verbosity) is -1.

[0085] RandomForest:

[0086] Preferred parameters include:

[0087] The number of trees (n_estimators) is 200-1000, and the maximum tree depth (max_depth) can be empty. In this embodiment, the number of trees is 300, and the maximum tree depth is 12; further settings include: minimum number of samples per leaf node (min_samples_leaf) is 1, random seed (random_state) is RANDOM_STATE, and the number of parallel threads (n_jobs) is -1.

[0088] CatBoost:

[0089] Preferred parameters include:

[0090] The number of iterations is 500-2000, and the learning rate is 0.01-0.1.

[0091] In this embodiment, the number of iterations is set to 500, the learning rate is 0.05, and the following settings are further defined: tree depth is 6, L2 leaf node regularization coefficient (l2_leaf_reg) is 3, sample sampling ratio (subsample) is 0.8, output log switch (verbose) is False, and random seed (random_seed) is RANDOM_STATE.

[0092] S43. Stacking Training (Stacked Generalization Training): XGBoost, LightGBM, RandomForest, and CatBoost are selected as base learners, and linear regression is used as the meta-learner. The training set is divided into 5 mutually exclusive subsets. For any k-th fold, each base learner is trained using the samples from the other 4 folds, and predictions are made for the k-th fold sample to obtain the corresponding out-of-fold prediction value. After traversing all 5 folds, the out-of-fold prediction results of each training sample under each base learner are obtained, and concatenated according to the learner dimension to form a meta-feature matrix.

[0093] Using the meta-feature matrix as input and the corresponding observed runoff as output, the meta-learner is fitted and trained to learn the fusion weight relationship of each base learner. The meta-learner preferably uses ordinary least squares linear regression, which by default includes an intercept term.

[0094] Subsequently, each base learner was retrained using all the training samples to obtain a set of base learners for deployment prediction.

[0095] For any sample to be predicted, first input the retrained base learners to obtain the predicted values ​​of each base learner, then concatenate the predicted values ​​output by each base learner into a meta-feature vector and input it into the meta-learner to output the final runoff prediction result.

[0096] The simulation results of each model in the Shuimingya watershed are shown in Table 2. The Stacking ensemble model has the highest accuracy and is effective.

[0097] Table 2. Accuracy Indicators of Machine Learning and Stacking Ensemble Models in the Shuimingya Watershed

[0098]

[0099] S5. HyPS-IF: Interpretability analysis, factor composite importance, and response threshold extraction

[0100] S51. Interpretation Toolkit: Calculates SHAP (sample-level SHAP value), PFI (permutation importance), PDP (partial dependency), and ICE (individual condition expectation) for the Stacking model. The SHAP results for the Shuimingya watershed are shown in the figure below. Figure 4 As shown.

[0101] S52. HRU Level Contribution Vector and Feature Factor Comprehensive Importance Score: For the first... Construct HRU-level contribution vectors for each HRU. Its elements are the representative contributions of each feature factor in the joint training feature set to the HRU sample set.

[0102] Preferably, the mean of the absolute values ​​of the SHAP contribution values ​​corresponding to all samples within the HRU is used to construct the SHAP contribution value. in, denoted as the mean absolute value of the SHAP contribution of the m-th feature factor included in the joint training feature set and used as input to the Stacking ensemble model on the h-th HRU, used to measure the strength of the contribution of the factor to the model prediction result; m is the number of feature factors.

[0103] Based on this, SHAP, PFI, PDP, and ICE are normalized and a multidimensional explanatory index vector is constructed. The overall importance score of HRU is then calculated through weighted aggregation.

[0104]

[0105] in, For the joint training feature set, the first The overall importance score of each feature factor. , For the SHAP and PFI values ​​normalized to [0,1] respectively, For the joint training feature set, the first The goodness of fit of the PDP for each feature factor. Weights are a measure of ICE curve consistency. Based on the validation set, it is automatically determined and The preferred weight value is (Can be dynamically adjusted according to verification accuracy) The sum is 1, and SHAP is the main method.

[0106] S53. Response Threshold Extraction: For feature factors whose overall importance score in the joint training feature set is higher than a preset threshold, the response threshold is extracted from their PDP curve and ICE curve. Specifically, the PDP curve of the feature factor is first fitted with a smooth spline, and the absolute value of the second derivative of the fitted PDP curve at each feature value is calculated. The 95th percentile of the absolute value sequence of the second derivative is used as the curvature screening threshold, and the feature value positions corresponding to the absolute value of the second derivative being greater than or equal to the curvature screening threshold are determined as the first candidate response threshold. Further, the 25th percentile and 75th percentile of the feature factor sample value distribution are determined as the second candidate response threshold. The first and second candidate response thresholds are merged to form a candidate response threshold set. The 95th percentile is a statistic based on the absolute value sequence of the second derivative of the PDP curve, used to screen for locations of curvature abrupt changes; the 25th and 75th percentiles are statistics based on the distribution of sample values ​​for this feature factor, used to supplement the representative thresholds reflecting the data distribution. These two are statistics from different objects. The PDP / ICE curves of the high comprehensive importance feature factor in the Shuimingya watershed are shown below. Figure 5 As shown. Figure 5 In the PDP curve, the overall trend of the influence of the corresponding characteristic factor value changes on the predicted runoff response is reflected, while the ICE curve reflects the consistency and heterogeneity of individual responses of different samples. When the PDP curve shows a significant turning point in a certain value range, and the corresponding ICE curve shows high consistency in that range, the values ​​near that range can be used as important candidates for the response threshold.

[0107] S6. K-means clustering partitioning based on contribution vectors

[0108] In this invention, "partitioning" refers to the process of merging multiple HRUs with similar runoff response mechanisms into clustered regions based on the similarity of their HRU-level contribution vectors. HRUs are the basic units of a watershed, and partitioning is the result of clustering one or more HRUs; each HRU belongs to only one partition.

[0109] S61. Using HRUs as clustering samples, use the K-means method to calculate the HRU-level contribution vector for each HRU. Clustering is performed to group HRUs with similar contribution patterns into the same partition, thereby obtaining data-driven and interpretable runoff response partitions.

[0110] S62. K value selection: The Elbow method is used to determine the candidate optimal number of clusters based on the inflection point of the intra-cluster sum of squares (SSE) as a function of the cluster number. Based on the candidate optimal cluster number, upper and lower limits are set according to the engineering application requirements, and the stability and interpretability of the partitioning results are manually verified.

[0111] S63. Cluster Stability Test: Perform bootstrapping sampling on the clustering input and repeatedly execute K-means clustering, calculating the ARI (Aggregate Membership Integrity Index) among different repeated results; when the ARI is lower than a preset threshold, adjust the number of clusters. Alternatively, adjust the composition of the clustering input vector.

[0112] S64. Partition Output: Output the partition number, partition centroid, HRU list within the partition, and partition statistical characteristics for each partition. The partition centroid is the mean vector of the contribution vectors of each HRU level within the partition, used to characterize the representative runoff response mechanism of the partition; the partition statistical characteristics include sample size, average runoff coefficient, coefficient of variation, and importance ranking of feature factors.

[0113] S7. Unified Model Selection and Parameterization After Partitioning

[0114] S71. Candidate Model Pool: Retains multiple candidate models (GXAJ, XGBoost, CatBoost, LightGBM, RandomForest, etc.).

[0115] S72. Unified evaluation within the region: For each region, summarize historical samples within the region and perform unified cross-validation evaluation on candidate models (test the NSE, KGE, RMSE and PBIAS of each candidate model respectively), and determine the preferred model set within the region accordingly.

[0116] S73. Example of unified model strategy: If a candidate model has a significantly higher average NSE and PBIAS within the historical window within a partition, it is selected as the preferred model; if multiple candidate models have similar performance, they are weighted and stacked according to their historical performance to output the final prediction.

[0117] S74. Parameter migration principle: Use the parameter weighted average of the partition centroid or the K nearest neighbor HRU in the cluster as the initial guess of the unobserved HRU, and then use the lightweight optimization SCE-UA at the partition level for fine-tuning.

[0118] S75. Output: For each region, output daily runoff forecasts, fusion / preferred model identifiers, model parameter sets, forecast confidence intervals, and HyPS-IF interpretation reports.

[0119] Table 3 presents the cross-validation results of candidate models in different zones of the Shuimingya watershed.

[0120] Table 3. Selection and Performance Comparison of Zoning Models for the Shuimingya Watershed

[0121]

[0122] S8. Online monitoring, triggering rules, and iterative optimization

[0123] Monitoring metrics: Calculate NSE, RMSE, KGE, PBIAS and confidence interval coverage for each partition on a monthly or preset periodic basis; record historical trends and save versions.

[0124] Example of triggering conditions: Automatic iteration is triggered when NSE < γ for M consecutive months (preferably γ = 0.7, M = 2) or RMSE increases by > p% (preferably p = 20%) compared to the baseline in any month. This can be configured according to production requirements.

[0125] Experimental results demonstrate that this invention can optimize the combined use of hydrological models by deeply exploring the response relationships between environmental factors, physical priors, and runoff. Through comprehensive analysis of multiple typical watersheds, the method of this invention not only improves the flexibility and adaptability of the models but also ensures high-accuracy hydrological simulation results even under conditions of limited data and variable environmental conditions. In the future, this method is expected to play an important role in flood forecasting, watershed management, and ecological protection, providing scientific basis and technical support for the sustainable use of water resources.

[0126] Example 2

[0127] This embodiment provides a target watershed runoff prediction device, including a processor and a memory; the memory stores a program or instructions, which are loaded and executed by the processor to implement the steps of the target watershed runoff prediction method provided above.

[0128] Example 3

[0129] This embodiment provides a computer-readable storage medium storing a program or instructions that, when executed by a processor, implement the steps of the target watershed runoff prediction method described above.

[0130] While the present invention has been described above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A runoff prediction method integrating physical priors and machine learning, characterized in that, Includes the following steps: S1. Divide the target watershed into basic watershed units and collect multi-source data for each basic watershed unit; S2. Based on the collected multi-source data, run various types of hydrological models for each basic watershed unit to extract candidate physical prior factors; among which, the different types of hydrological models include physical process hydrological models, semi-physical and semi-empirical hydrological models, and conceptual hydrological models; S3. Perform single-factor perturbation sensitivity analysis on the extracted candidate physical prior factors. Select core physical prior factors from the candidate physical prior factors to construct a joint training feature set. Apply a single-factor perturbation of ±δ% to each candidate physical prior factor and calculate the relative change ΔNSE between the baseline NSE of the corresponding hydrological model under the baseline operating conditions and the NSE after applying the perturbation. When ΔNSE ≥ 5%, the candidate physical prior factor is determined to be a core physical prior factor, where δ is 5% to 20%. Based on the determined core physical prior factors, statistical derived features and ratio derived features are constructed. One-hot coding is performed on the categorical features to obtain categorical coding features. Z-score standardization is performed on the meteorological time series continuous variables, static topographic continuous features, candidate physical prior factors, and the statistical and ratio-type derived features of the candidate physical prior factors. The above-mentioned features after standardization are used to construct a joint training feature vector, and the joint training feature vector constitutes a joint training feature set. The joint training feature vector includes time-varying meteorological and hydrological continuous variables and their sliding window statistical features, precipitation lag accumulation features, static topographic continuous features, category coding features, as well as the core physical prior factors for judgment and the statistical and ratio-based derived features constructed from the core physical prior factors. S4. Construct a Stacking ensemble model; using the joint training feature set as input, train the Stacking ensemble model for daily runoff prediction, and calculate the percentage bias index (PBIAS) after training to verify the physical consistency of the water balance in the prediction results; the calculation method for the percentage bias index (PBIAS) is as follows: in, The first Daily runoff forecast and observed values, where n is the total number of daily-scale samples involved in the PBIAS calculation; S5. The trained Stacking ensemble model is subjected to interpretability analysis using an ensemble interpretation framework. The comprehensive importance score of each feature factor in the joint training feature set is calculated. The response threshold of the feature factors in the joint training feature set whose comprehensive importance score is higher than a preset threshold is extracted. Based on the interpretability analysis results of the comprehensive importance score, response threshold and single factor contribution intensity, an HRU-level contribution vector is constructed for each watershed basic unit. S6. Based on the HRU-level contribution vectors of each watershed basic unit, the K-means method is used to perform clustering partitioning of the target watershed in units of watershed basic units; S7. For each partition after clustering, evaluate within the candidate model pool consisting of the various hydrological models of different categories run in step S2 and the Stacking ensemble model trained in step S4, and determine the preferred model set for that partition; output the daily runoff simulation results based on the preferred model set determined for each partition.

2. The runoff prediction method according to claim 1, characterized in that, In step S2, based on the collected multi-source data, multiple hydrological models of different categories are run for each basic watershed unit to extract candidate physical prior factors. This includes: uniformly processing the time-varying meteorological and hydrological continuous variables in the multi-source data into daily-scale time series, and preprocessing the static topographic continuous features and category features in the multi-source data; running the selected hydrological model for each basic watershed unit within a preset time range, and extracting hydrological process variables from the output of each hydrological model as candidate physical prior factors during the model's operation; and constructing multi-scale sliding window statistical features using the time-varying meteorological and hydrological continuous variables and the candidate physical prior factors as continuous input variables; wherein the continuous input variables include precipitation. Potential evaporation Daily maximum temperature Daily minimum temperature Daily average temperature Dew point temperature Wind speed Solar radiation thermal radiation relative humidity Net radiation , , and the extracted candidate physical prior factors; For any continuous input variable, calculate the past... The sliding window statistical characteristics within the daily window include cumulative value, mean, standard deviation, maximum value, minimum value, and linear slope, with W taking values ​​of 5, 15, and 30. Construct precipitation lag cumulative features for precipitation variables in continuous input variables.

3. The runoff prediction method according to claim 1, characterized in that, In step S4, the Stacking ensemble model includes a base learner and a meta learner; the base learner includes two or more of XGBoost, CatBoost, LightGBM and Random Forest, and the meta learner is linear regression or Lasso; K-fold cross-validation is used to optimize the parameters of the base learner and the meta learner, and an early stopping strategy is used for learners that support iterative training to prevent overfitting, where K is 3 to 10.

4. The runoff prediction method according to claim 1, characterized in that, In step S5, the method for extracting the response threshold of feature factors whose comprehensive importance score in the joint training feature set is higher than a preset threshold is as follows: The PDP curve of the feature factor is fitted with a smooth spline. The absolute value of the second derivative of the fitted PDP curve at each feature value is calculated. The 95th percentile of the absolute value sequence of the second derivative is used as the curvature screening threshold. The feature value position corresponding to the absolute value of the second derivative being greater than or equal to the curvature screening threshold is determined as the first candidate response threshold. The 25th and 75th percentile values ​​of the feature factor sample value distribution are then determined as the second candidate response thresholds; the first and second candidate response thresholds are merged to form a candidate response threshold set.

5. The runoff prediction method according to claim 4, characterized in that, In step S5, the comprehensive importance score of each feature factor in the joint training feature set is: in, For the joint training feature set, the first The overall importance score of each feature factor. , For the SHAP and PFI values ​​normalized to [0,1] respectively, As a factor The goodness of fit of PDP, Weights are a measure of ICE curve consistency. Based on the validation set, it is automatically determined and .

6. The runoff prediction method according to any one of claims 1-5, characterized in that, The process after step S7 also includes: S8. Provide online accuracy feedback based on preset performance trigger conditions and implement iterative optimization of re-screening, retraining, and partitioning and model adaptation when triggered; In step S8, the online accuracy triggering condition is: when NSE < γ for M consecutive months or the RMSE of any month increases by more than p% relative to the baseline, automatic iterative optimization is triggered; after optimization, it needs to be verified by at least 30 days of actual test data and the NSE during the verification period must be ≥0.75 before the optimization result can be solidified; p=20%, γ = 0.7, M = 2.

7. A runoff prediction device, characterized in that, It includes a processor and a memory; the memory stores a program or instructions which are loaded and executed by the processor to implement the steps of the runoff prediction method as described in any one of claims 1 to 6.

8. A computer-readable storage medium, characterized in that, The readable storage medium stores a program or instructions that, when executed by a processor, implement the steps of the runoff prediction method as described in any one of claims 1 to 6.

Images