River ice thickness interpretable prediction and early warning method based on multi-scale hydrothermal characteristics

Abstract

This invention relates to the field of cold-region river ice monitoring, prediction, and early warning technology, and discloses an interpretable prediction and early warning method for river ice thickness based on multi-scale hydrothermal features. The method includes: collecting hourly air temperature, ice surface temperature, water level elevation, and ice thickness data; constructing a multi-scale hydrothermal feature matrix based on air temperature, ice surface temperature, and water level, including real-time observation features, rolling statistical features of air temperature, rolling mean features of ice surface temperature, cumulative freezing time features across multiple time scales, and temperature change rate features; training a random forest full-feature model using this feature matrix as input; calculating the mean absolute value of each feature's SHAP value and sorting them from largest to smallest, thus classifying the features into factors, auxiliary factors, and negligible factors. The optimal feature subset is determined through recursive feature elimination cross-validation and forward feature selection, and a simplified forecast model is obtained through retraining. ALE curves are generated for the factors, and physical control parameters such as the freezing energy critical threshold, saturation threshold, water level suppression inflection point, and optimal ice-inducing temperature are extracted by identifying points where the first derivative, second derivative, and peak points of the curves are zero. Based on these physical control parameters, hierarchical early warning rules are established for cold accumulation state classification, rapid thickening co-existence, high water level instability and thermal-hydrological offsetting, and multi-scale cold accumulation synergistic enhancement. Residual analysis, autocorrelation testing, interval error statistics, and model validation consistency verification are performed on the simplified forecast model. This invention achieves interpretability and automatic extraction of physical parameters in river ice thickness prediction, improving the targeting and transferability of early warnings.

Description

Technical Field

[0001] This invention relates to the field of river ice monitoring, prediction and early warning technology in cold regions, specifically to a method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics. Background Technology

[0002] The evolution of ice conditions in cold-region rivers directly affects the safe operation of water conservancy projects, ice disaster prevention, and water resource allocation during the ice season. Ice thickness, as a core characteristic parameter of ice condition evolution, is the foundation for accurate prediction of ice forecasts and disaster warnings. Currently, ice thickness prediction mainly relies on two types of methods: empirical formulas and physical process models. Empirical formulas establish statistical relationships between ice thickness and a few environmental factors such as temperature and number of freezing days through regression analysis. They are simple in form, but these relationships are usually based on observational data of specific river sections. When moving to river sections with different climatic conditions or river morphologies, the formula parameters need to be recalibrated, resulting in poor adaptability. Physical process models, based on thermodynamic and hydraulic principles, describe the physical processes of ice growth and melting. Theoretically, they have high universality, but the models involve a large number of parameters. Some parameters, such as the ice surface heat exchange coefficient and water heat flux, are difficult to obtain directly through conventional observations and need to be determined by empirical values ​​or inversion methods. The model calibration process is complex, limiting its application in practical operations.

[0003] In recent years, machine learning methods have been gradually applied in river ice thickness prediction. Data-driven models such as random forests and support vector machines can automatically learn the nonlinear mapping relationship between environmental factors and ice thickness, and their prediction accuracy is usually better than that of traditional empirical formulas. However, these methods have obvious limitations. The internal mechanisms of the models are not transparent, making it difficult to explain the independent contribution of each environmental factor to the evolution of ice thickness. The model prediction results cannot be correlated with physical processes such as ice surface heat exchange and cold accumulation. The models themselves do not produce any physically meaningful parameters, which makes it difficult for ice condition forecasters to understand and trust the model prediction results, and also fails to provide a physical basis for decision-making for ice jam disaster early warning. In addition, existing early warning methods are mostly based on a single threshold or simple empirical rules, and lack an effective joint judgment mechanism when the influence of multiple environmental factors on ice thickness is inconsistent. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides an interpretable prediction and early warning method for river ice thickness based on multi-scale hydrothermal characteristics. This method solves the problems of existing river ice thickness prediction methods, such as difficulty in automatically extracting physically meaningful control parameters from models, reliance on human experience in obtaining physical parameters, poor parameter mobility between different river sections, and lack of a multi-factor joint early warning mechanism.

[0005] To achieve the above objectives, the present invention provides an interpretable prediction and early warning method for river ice thickness based on multi-scale hydrothermal characteristics, comprising the following steps: S1, Multi-parameter Continuous Monitoring Data Acquisition: A fixed-point multi-parameter integrated monitoring system is deployed in the target cold-region river section to collect hourly air temperature, hourly ice surface temperature, hourly water level elevation, and hourly ice thickness data. The sampling frequency is once per hour, with continuous operation for at least three complete ice-covering seasons to obtain continuous observation data covering the entire cycle of ice condition evolution. The raw data is cleaned, including outlier removal based on interquartile range and missing value filling using linear interpolation, and the cleaned data is then time-series aligned.

[0006] S2, a multi-scale hydrothermal feature system, is constructed based on hourly air temperature, hourly ice surface temperature, and hourly water level elevation data collected in S1, forming a 26-dimensional multi-scale hydrothermal feature matrix. This matrix includes five types of features: real-time observation features (real-time air temperature, real-time ice surface temperature, and real-time water level elevation); rolling statistical features of air temperature (rolling mean and rolling minimum values ​​for window lengths of 6, 12, 24, 72, and 168 hours); rolling mean features of ice surface temperature (rolling mean values ​​for window lengths of 6, 12, 24, and 72 hours); cumulative freezing time features across multiple time scales (cumulative freezing time for 1, 3, 7, 14, and 30 days); and temperature change rate features (1-hour, 24-hour, 1-hour, and 24-hour rates of change for air temperature and ice surface temperature). This feature matrix only includes external environmental variables and excludes lagged ice thickness variables to avoid interference from autoregressive information in the physical analysis.

[0007] S3, the full-feature random forest model training, uses the 26-dimensional feature matrix constructed in S2 as input and hourly ice thickness as the target variable to train a random forest regression model, obtaining a full-feature ice thickness prediction model. The dataset is divided into training and test sets in an 8:2 ratio. 500 decision trees are generated by sampling using Bootstrap, and the model output is the average of the prediction results of all decision trees.

[0008] S4, based on SHAP-based global importance ranking and factor decoupling, uses the TreeSHAP algorithm to calculate the SHAP value for each feature of each sample in the full-feature random forest model trained in S3. The SHAP value decomposes the model's prediction for each sample into the sum of the baseline prediction and the contribution of each feature; a positive value indicates that the feature makes the prediction higher than the baseline, and a negative value indicates that the prediction is lower than the baseline. For each feature, the mean of its absolute SHAP values ​​across all samples is calculated as the global importance index for that feature. All 26 features are ranked from largest to smallest global importance, and based on the ranking, features are divided into three levels: Factors are features ranked in the top 5 and whose global importance is more than twice that of the 6th ranked feature; Auxiliary factors are features ranked from 6th to 15th and whose global importance is greater than that of the last ranked feature; Negligible factors are features ranked from 16th to 26th and whose global importance is less than 10% of that of the 1st ranked feature.

[0009] S5, Feature Selection and Simplified Prediction Model Construction: Based on the SHAP global importance ranking results from S4, recursive feature elimination and cross-validation are performed. Starting from the full feature set, each iteration removes the feature with the lowest SHAP global importance. The remaining features are used to retrain the random forest model, and the model performance is evaluated using 5-fold cross-validation. The average coefficient of determination and average root mean square error are recorded after each iteration. The optimal number of features is selected when the average coefficient of determination reaches its maximum value and the average root mean square error reaches its minimum value. Simultaneously, forward feature selection is performed, based on the SHAP global importance ranking, starting with the most important feature, adding one feature at a time. The model is retrained and evaluated on the test set. Adding features stops when the increase in the coefficient of determination on the test set after adding a new feature is less than 0.005. This subset of features is called the simplified feature subset. The intersection of the feature subset corresponding to the optimal number of features and the simplified feature subset is taken as the final feature subset. The random forest regression model is retrained using this final feature subset to obtain the simplified prediction model.

[0010] S6. Physical control parameter extraction based on the geometric features of the ALE curve. Factors are selected as the analysis objects for the ALE curve, including 30-day cumulative freezing time, 14-day cumulative freezing time, water level elevation, 168-hour rolling minimum air temperature, and 72-hour rolling average ice surface temperature. A one-dimensional cumulative local effect curve is calculated, which describes the average local influence of a single feature on the model's predicted value, eliminating bias caused by the correlation between features. The ALE curve is then smoothed using a sliding window with a width of 5 intervals. The first and second derivatives of the smoothed ALE curve are calculated. The first derivative reflects the changing trend of the ALE curve, and the second derivative reflects the inflection point of the changing trend. Based on the calculation results of the first and second derivatives, three types of geometric feature points are identified: points where the first derivative is zero, corresponding to local extreme points of the ALE curve; points where the second derivative is zero, corresponding to inflection points of the changing trend of the ALE curve; and curve peak points, i.e., local maximum points on the ALE curve. The precise location of the geometric feature points in the feature value space is obtained through linear interpolation.

[0011] A mapping rule is established between physical control parameters and the geometric characteristics of the ALE curve. For the 30-day cumulative freezing time characteristic, the point where the first derivative of the ALE curve is zero and the ALE value changes from negative to positive is identified, and the characteristic value corresponding to this point is denoted as the freezing energy critical threshold. The point where the second derivative of the ALE curve changes from rising to falling and is zero is identified, and the characteristic value corresponding to this point is denoted as the freezing energy saturation threshold. Based on the freezing energy critical threshold, 50% of this value is calculated as the slow growth initiation threshold. For the water level elevation characteristic, the point where the first derivative of the ALE curve is zero and the ALE value reaches its maximum value is identified, and the characteristic value corresponding to this point is denoted as the water level inhibition inflection point. The water level interval where the ALE value is greater than 90% of the ALE value corresponding to the water level inhibition inflection point is identified and denoted as the optimal hydrological window, and its lower limit and upper limit are denoted as the lower limit and upper limit of the optimal hydrological window, respectively. For the 168-hour rolling minimum temperature feature, the point where the first derivative of the ALE curve is zero and the ALE value reaches its maximum value is identified; the feature value corresponding to this point is denoted as the optimal ice-inducing temperature. Two feature values ​​are identified when the ALE value decreases from its maximum value to 50% of its maximum value; the smaller value is denoted as the temperature suppression threshold, and the larger value as the temperature upper limit threshold. For the 14-day and 30-day cumulative freezing time features, the integral values ​​of their ALE curves within their respective feature value ranges are calculated, and the ratio of the two is used as a discriminant variable for the co-enhancing model. The ratio of the larger to the smaller value is used as a discriminant variable for the single-scale dominant model. The maximum absolute value of the 30-day cumulative freezing time feature ALE curve and the maximum absolute value of the water level elevation feature ALE curve are calculated, and the thermo-hydrological offsetting intensity coefficient is calculated based on these two values.

[0012] S7, Generation of Graded Early Warning Rules Based on Physical Parameters: Based on the set of physical control parameters extracted in S6, the following early warning rules are established: Graded Monitoring and Early Warning Rules for Cold Accumulation Status: Based on the comparison between the current 30-day cumulative freezing time and the slow growth initiation threshold, the freezing energy critical threshold, and the freezing energy saturation threshold, the cold accumulation status is divided into four states: severely insufficient cold energy, slow growth, rapid thickening, and thickening saturation, each corresponding to different monitoring frequencies and early warning colors. Joint Early Warning Rules for Rapid Thickening: When the current 30-day cumulative freezing time is greater than or equal to the freezing energy critical threshold, and the change in the 30-day cumulative freezing time over the past 24 hours is greater than zero, and the current 168-hour rolling minimum temperature is between the temperature suppression threshold and the optimal ice-promoting temperature, a joint early warning for rapid thickening is triggered. High water level instability warning and thermal-hydrological offset warning rules: The high water level instability warning is triggered when the current water level elevation is greater than or equal to the water level suppression inflection point; the thermal-hydrological offset warning is triggered when the current 30-day cumulative freezing time is greater than or equal to the freezing energy critical threshold, the current water level elevation is greater than or equal to the upper limit of the optimal hydrological window, and the thermal-hydrological offset intensity coefficient is greater than 0.3. Multi-scale cold accumulation synergistic enhancement warning rules: When the discriminant variable of the synergistic enhancement mode is greater than 1 and the discriminant variable of the single-scale dominant mode is between 0.5 and 2, it is determined that there is a multi-scale cold accumulation synergistic enhancement effect, triggering a red warning, increasing the monitoring frequency to once / 3 hours, and extending the warning lead time by 12 hours. The localization and migration mechanism of all warning thresholds is as follows: When this method is migrated to a new river section, steps S1 to S6 are repeated, the model is retrained using the observation data of the new river section, and the physical control parameters are extracted to obtain the localized warning thresholds applicable to the new river section.

[0013] S8, Model Validation Consistency Check: The simplified forecast model constructed in S5 is checked as follows: Calculate the test set residuals, plot the residual histogram and normal probability plot, calculate the residual mean and standard deviation, and check whether the residuals follow a normal distribution; use the Durbin-Watson statistic to test the first-order autocorrelation of the residual sequence, and determine that there is no significant autocorrelation when the statistic is between 1.5 and 2.5; divide the ice thickness into intervals with a step size of 0.1 meters according to the measured value, and calculate the root mean square error and mean absolute error in each interval; perform model validation consistency check, verifying whether the inequality constraints such as the slow growth initiation threshold being less than the freezing energy critical threshold, the freezing energy critical threshold being less than the freezing energy saturation threshold, the lower limit of the optimal hydrological window being less than the upper limit of the optimal hydrological window, the temperature suppression threshold being less than the optimal ice-promoting temperature, and the optimal ice-promoting temperature being less than the upper limit of the temperature threshold hold, and at the same time verifying whether the correspondence between the position of the geometric feature points of the ALE curve and the physical prediction is reasonable.

[0014] This invention provides an interpretable prediction and early warning method for river ice thickness based on multi-scale hydrothermal characteristics. It has the following beneficial effects: 1. This invention constructs a multi-scale hydrothermal feature system that includes real-time observation features, rolling temperature statistics, rolling mean ice surface temperature, cumulative freezing time features across multiple time scales, and temperature change rate features. This system integrates environmental information from different time scales into the model input, enabling the model to simultaneously perceive the impact of short-term temperature fluctuations and long-term cold accumulation on ice thickness evolution. Based on this, the importance of features is globally ranked according to SHAP values, and physical control parameters such as the critical threshold for freezing energy, the saturation threshold for freezing energy, the water level inhibition inflection point, and the optimal ice-inducing temperature are extracted by combining the geometric features of the ALE curve. This makes the model's prediction results interpretable and overcomes the shortcomings of traditional machine learning models, such as opaque prediction mechanisms and difficulty in extracting physical parameters from within the model.

[0015] 2. This invention determines the optimal feature subset by combining recursive feature elimination cross-validation with forward feature selection. This reduces feature redundancy and model complexity while ensuring prediction accuracy, and improves the model's computational efficiency and generalization ability. Based on the extracted physical control parameters, it establishes hierarchical monitoring and early warning rules for cold accumulation, joint early warning rules for rapid thickening, early warning rules for high water level instability and thermo-hydrological counterbalancing, and multi-scale cold accumulation collaborative enhancement early warning rules. This enables multi-level differentiated early warning for different ice condition development stages, improving the pertinence and accuracy of the early warning.

[0016] 3. This invention verifies the consistency between the simplified forecast model's prediction results and physical laws through model validation consistency check, ensuring that the model output conforms to the basic constraints of thermodynamics and hydraulics, thereby improving the model's credibility and reliability. All warning thresholds are automatically extracted from the observation data by the algorithm without relying on manual presets. When migrating to a new river section, repeating steps S1 to S6 can obtain localized warning thresholds, enabling this method to adapt to the differences in hydrothermal conditions in different river sections and exhibiting good portability. Attached Figure Description

[0017] Figure 1 This is the overall flowchart of the present invention; Figure 2 This is a flowchart of the monitoring data acquisition process of the present invention; Figure 3 This is a flowchart illustrating the construction process of the multi-scale hydrothermal feature system of the present invention. Figure 4 This is a flowchart of the full-feature model training process of the present invention; Figure 5 This is a flowchart of the global importance ranking and factor decoupling process of the present invention; Figure 6 This is a flowchart illustrating the construction process of the simplified forecast model of the present invention; Figure 7This is a flowchart of the physical control parameter extraction process of the present invention; Figure 8 A flowchart for generating the hierarchical early warning rules of this invention; Figure 9 This is a flowchart of the model verification consistency check for the present invention. Detailed Implementation

[0018] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. 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 are within the scope of protection of the present invention.

[0019] See attached document Figure 1 , Figure 1 This is a flowchart illustrating the overall process of a river ice thickness prediction and early warning method based on multi-scale hydrothermal feature parameters and interpretable machine learning, according to an embodiment of the present invention. The present invention provides a river ice thickness prediction and early warning method based on multi-scale hydrothermal feature parameters and interpretable machine learning, comprising the following steps: S1, Multi-parameter continuous monitoring data acquisition: Deploy a fixed-point multi-parameter integrated monitoring system in the target cold-region river section to collect hourly air temperature, hourly ice surface temperature, hourly water level elevation and hourly ice thickness data. The sampling frequency is 1 time / hour and the continuous operation time is no less than 3 complete ice-covered seasons. S2, Construction of Multi-Scale Hydrothermal Feature System: Based on the hourly air temperature, hourly ice surface temperature and hourly water level elevation data collected in S1, a 26-dimensional multi-scale hydrothermal feature matrix is ​​constructed. This feature matrix includes five types of features: real-time observation features, rolling statistical features of air temperature, rolling mean features of ice surface temperature, cumulative freezing degree hours features across multiple time scales, and temperature change rate features. S3, Random Forest Full Feature Model Training: Using the 26-dimensional feature matrix constructed in S2 as input and hourly ice thickness as the target variable, a random forest regression model is trained to obtain a full feature ice thickness prediction model. The dataset is divided into training and test sets in an 8:2 ratio. S4, Global Importance Ranking and Factor Decoupling Based on SHAP: For the full-feature random forest model trained in S3, calculate the absolute value of the SHAP mean of all 26 features, sort them from largest to smallest, quantify the independent contribution of each environmental factor to ice thickness evolution, and divide the features into three levels: factors, auxiliary factors and negligible factors based on the ranking results. S5, Feature Selection and Construction of Simplified Prediction Model: Based on the SHAP importance ranking results of S4, the recursive feature elimination cross-validation and forward feature selection methods are used to determine the optimal feature subset that maintains prediction accuracy with minimal feature redundancy. The random forest regression model is then retrained with this optimal feature subset to obtain the simplified prediction model. S6, Physical control parameter extraction based on the geometric features of the ALE curve: The cumulative local effect method is used to generate a one-dimensional ALE curve for the factors in the full feature model of S3. By identifying the points where the first derivative of the ALE curve is zero, the points where the second derivative is zero, and the peak points of the curve, physical control parameters such as the critical threshold of freezing energy, the saturation threshold of freezing energy, the water level inhibition inflection point, and the optimal ice-promoting temperature are extracted. S7, Generation of hierarchical early warning rules based on physical parameters: Based on the set of physical control parameters extracted in S6, establish hierarchical monitoring and early warning rules for cold accumulation status, joint early warning rules for rapid thickening, early warning rules for high water level instability and thermo-hydrological offsetting, and multi-scale cold accumulation synergistic enhancement early warning rules. S8, Model Validation Consistency Check: Perform residual distribution analysis, residual autocorrelation test, ice thickness interval error statistics and model validation consistency check on the simplified forecast model constructed in S5 to verify the model's prediction accuracy.

[0020] See attached document Figure 2 , Figure 2 This is a flowchart illustrating the overall process of a river ice thickness prediction and early warning method based on multi-scale hydrothermal characteristic parameters and interpretable machine learning, according to an embodiment of the present invention. The following provides a detailed description of the S1 multi-parameter continuous monitoring data acquisition step.

[0021] S101, a fixed-point multi-parameter integrated monitoring system is deployed in the target cold-region river section. This system includes four types of sensors: an air temperature sensor for hourly air temperature, an ice surface temperature sensor for hourly ice surface temperature, a water level sensor for hourly water level elevation, and an ice thickness measurement device for hourly ice thickness. The air temperature sensor is a platinum resistance temperature sensor installed in a Stevenson screen, with a measurement accuracy of ±0.1℃. The ice surface temperature sensor is a non-contact infrared temperature sensor, installed on a fixed support on the riverbank, vertically aligned with the center of the ice surface, with a measurement accuracy of ±0.3℃. The water level sensor is a pressure-type water level gauge, installed at a fixed position on the bottom of the riverbed, with a measurement accuracy of ±0.5cm. The ice thickness measurement device uses air-coupled ground-penetrating radar, installed on a fixed support on the riverbank, with the transmitting and receiving antennas arranged parallel to the ice surface, a detection depth range of 0~3m, and a measurement accuracy of ±1cm. For the specific working principle and data processing methods of the air-coupled ground-penetrating radar in the ice thickness measurement device, those skilled in the art can refer to relevant technical standards; these are well-known technologies in the field and will not be elaborated further here.

[0022] S102 is configured with a sampling frequency of once per hour, with all sensors collecting data synchronously to form an hourly observation sequence. It operates continuously for at least three complete ice-covering seasons to obtain continuous observation data covering the entire cycle of ice condition evolution. The collected data is uploaded to the data center for storage in real time via a wireless transmission module.

[0023] S103, clean the raw collected data. Outlier removal uses an interquartile range (IQR)-based method: for each variable, calculate its first quartile Q1 and third quartile Q3, defining the IQR as IQR = Q3 - Q1. Data points less than Q1 - 1.5 × IQR or greater than Q3 + 1.5 × IQR are marked as outliers and removed. Missing values ​​are handled using linear interpolation: for periods with no more than 6 consecutive missing values, linear interpolation is used to fill the missing value using two valid observations before and after it; for periods with more than 6 consecutive missing values, they are marked as invalid and removed from the dataset. After data cleaning, a valid sample set is formed, with each sample containing four synchronous observations: air temperature T_a (°C), ice surface temperature T_s (°C), water level elevation H (m), and ice thickness h (m). T_a, T_s, and H are input variables, and h is the target variable.

[0024] S104. Perform a time-series alignment check on the cleaned data. Since different sensors may have timestamp discrepancies, a sliding window matching method is used. Using the temperature sensor's timestamp as a reference, observations from other sensors whose deviations from the reference timestamp are within ±30 minutes are aligned to the reference timestamp. Data points outside this deviation range are marked as invalid. The aligned data sequence constitutes the input data for subsequent steps.

[0025] See attached document Figure 3 , Figure 3 This is a flowchart illustrating the construction of a multi-scale hydrothermal feature system according to an embodiment of the present invention. The following provides a detailed description of the steps involved in constructing the S2 multi-scale hydrothermal feature system.

[0026] S201, based on the hourly air temperature Ta, hourly ice surface temperature Ts, and hourly water level elevation H data collected and cleaned in step S1, constructs a 26-dimensional multi-scale hydrothermal feature matrix X. This feature matrix only includes external environmental variables and excludes the lagged variable of ice thickness h to avoid interference from autoregressive information in the physical analysis. The feature matrix contains 26 features in five categories, each corresponding to an independent physical quantity or statistical quantity.

[0027] S202, Construct real-time observation features, a total of 3. The real-time observation features are the original observation values ​​at the current time, including: real-time air temperature x1, which is the air temperature Ta at the current time; real-time ice surface temperature x2, which is the ice surface temperature Ts at the current time; and real-time water level elevation x3, which is the water level elevation H at the current time.

[0028] S203, constructing rolling statistical features of temperature, a total of 10. For the hourly temperature series Ta, each feature is calculated using a window length... Calculate the rolling mean and rolling minimum for the range {6, 12, 24, 72, 168} hours. The formula for calculating the rolling mean is: ; In the formula: For the current moment forward to the next moment Hourly temperature observations Window length (hours) For the temperature at the window The rolling average over an hour. The rolling minimum is denoted as min. t∈[1,w] This indicates the temperature at the window. The minimum value within an hour. The above features correspond to x4 to x13, specifically: x4 is the 6-hour rolling average of temperature, x5 is the 6-hour rolling minimum of temperature, x6 is the 12-hour rolling average of temperature, x7 is the 12-hour rolling minimum of temperature, x8 is the 24-hour rolling average of temperature, x9 is the 24-hour rolling minimum of temperature, x10 is the 72-hour rolling average of temperature, x11 is the 72-hour rolling minimum of temperature, x12 is the 168-hour rolling average of temperature, and x13 is the 168-hour rolling minimum of temperature.

[0029] S204, constructing four features for the rolling mean of ice surface temperature. For the hourly ice surface temperature sequence Ts, each feature is calculated using a window length... The rolling average for the intervals ∈{6, 12, 24, 72} is calculated using the following formula: ; In the formula: For the current moment forward to the next moment Hourly ice surface temperature observation values, For the ice surface temperature at the window The rolling average over the hour. The above features correspond to x14 to x17, specifically: x14 is the 6-hour rolling average of ice surface temperature, x15 is the 12-hour rolling average of ice surface temperature, x16 is the 24-hour rolling average of ice surface temperature, and x17 is the 72-hour rolling average of ice surface temperature.

[0030] S205 constructs a multi-timescale cumulative freezing hours (CFDH) feature, consisting of five features. The cumulative freezing hours are defined as the sum of the absolute values ​​of hourly negative air temperatures below 0°C within a specified timescale, characterizing the degree of cold accumulation. These features are defined by time scale. The cumulative freezing time is calculated for days ∈{1,3,7,14,30} using the following formula: ; In the formula: The time scale is (day). Indexed by hour. For the first Hourly temperature observation (°C) This means that a positive value is taken only when the temperature is below 0℃, and a value of 0 is taken when the temperature is above 0℃. Time scale The cumulative freezing time within a day, expressed in °C·h. The above features correspond to x18 to x22, specifically: x18 is the cumulative freezing time for 1 day, x19 is the cumulative freezing time for 3 days, x20 is the cumulative freezing time for 7 days, x21 is the cumulative freezing time for 14 days, and x22 is the cumulative freezing time for 30 days.

[0031] S206, construct four temperature change rate features. The temperature change rate reflects the short-term trend of air temperature and ice surface temperature, calculated using first-order difference. The 1-hour air temperature change rate x23 = Ta(t) - Ta(t-1), the 24-hour air temperature change rate x24 = Ta(t) - Ta(t-24), the 1-hour ice surface temperature change rate x25 = Ts(t) - Ts(t-1), and the 24-hour ice surface temperature change rate x26 = Ts(t) - Ts(t-24). Where Ta(t), Ta(t-1), and Ta(t-24) are the observed air temperatures at the current time, the previous hour, and the previous 24 hours, respectively; and Ts(t), Ts(t-1), and Ts(t-24) are the observed ice surface temperatures at the current time, the previous hour, and the previous 24 hours, respectively.

[0032] S207, Align the above five categories of 26 features by time to form a feature matrix X. X has a dimension of n×26, where n is the number of valid samples, and each row corresponds to a 26-dimensional feature vector at a given time. The feature matrix X and the corresponding target variable h (time-wise ice thickness) together constitute the input data for the subsequent step S3.

[0033] See attached document Figure 4 , Figure 4 This is a flowchart of the training process for a random forest full-feature model according to an embodiment of the present invention. The training steps of the S3 random forest full-feature model are described in detail below.

[0034] S301, Dataset Splitting. The feature matrix X constructed in step S2 and the corresponding target variable h (hourly ice thickness) are arranged in chronological order and split into training and test sets in an 8:2 ratio. The training set contains samples from the first 80% of the time series, and the test set contains samples from the last 20% of the time series. A fixed random seed is used to ensure the splitting results are repeatable. The training set is denoted as... The test set is denoted as ,in and The number of training samples and the number of test samples are respectively, satisfying the following conditions: : =8:2.

[0035] S302, Set the hyperparameters of the Random Forest model. Random Forest is an ensemble learning algorithm based on decision trees. It obtains the final prediction value by constructing multiple decision trees and averaging their predictions. The hyperparameters are set as follows: number of trees. =500, meaning 500 decision trees are built; max_depth=None, meaning there is no limit to the maximum depth of each decision tree, allowing it to grow until all leaf nodes are pure; min_samples_split=5, meaning the minimum number of samples required for each internal node to continue splitting is 5; min_samples_leaf=2, meaning the minimum number of samples contained in each leaf node is 2; max_features=26≈5, meaning 5 features are randomly selected from all 26 features for evaluation in each split; enable out-of-bag evaluation oob_score=True, meaning the model performance is estimated unbiasedly using out-of-bag samples (training samples not used by a decision tree).

[0036] S303, Training a Full-Feature Random Forest Regression Model. Using the training set feature matrix Xtrain as input and the training set target variable htrain as output, the random forest regression model is trained using the hyperparameters set in S302 to obtain the full-feature model f(X). During model training, 500 subsets are generated from the training set through Bootstrap sampling. A decision tree is trained on each subset, and the average of the predictions from all decision trees is used as the model output.

[0037] S304. Evaluate model performance using the test set. Input the test set feature matrix Xtest into the trained full-feature model f(X) to obtain the test set predictions. =f(Xtest). The following three metrics are used to evaluate model performance: Coefficient of Determination, ... The calculation formula is: ; In the formula: For the first Measured ice thickness values ​​for each test sample. For the first Predicted ice thickness for each test sample. = This is the average of the measured ice thickness values ​​for the test collection. The value of is in the range of (-∞, 1], and the closer it is to 1, the better the model fit. The formula for calculating the root mean square error (RMSE) is: ; The unit of RMSE is the same as that of ice thickness, in meters (m). A smaller RMSE value indicates a smaller prediction error. The formula for calculating the Mean Absolute Error (MAE) is: ; The unit of MAE is the same as that of ice thickness, which is meters (m). The smaller the value, the higher the prediction accuracy.

[0038] S305, Record the out-of-bag evaluation results. During model training, approximately 63.2% of the training samples for each decision tree are used to construct the tree, and the remaining approximately 36.8% of the samples are used as out-of-bag samples. After model training is complete, the out-of-bag evaluation results are calculated using all out-of-bag samples. Value, denoted as oob The calculation method is the same as in S304. The calculation method is the same, the difference being that out-of-bag samples are used instead of test set samples. oob It can be used as an unbiased estimate of the model's generalization performance and as a preliminary indicator of whether the model is overfitting.

[0039] See attached document Figure 5 , Figure 5 This is a flowchart of global importance ranking and factor decoupling based on SHAP according to an embodiment of the present invention. The following provides a detailed description of step S4, global importance ranking and factor decoupling based on SHAP.

[0040] S401, calculate the SHAP value for each feature of each sample in the full-feature model. SHAP (Shapley Additive Explanations) is an interpretability analysis method based on cooperative game theory. Its core idea is to decompose the model's prediction for each sample into the baseline prediction plus the sum of the contributions of each feature. For the random forest model, the TreeSHAP algorithm is used for efficient computation. This algorithm is specifically designed for tree models and can quickly and accurately calculate the SHAP value of each feature for each sample. For the i-th sample, the model prediction value... The relationship between the SHAP values ​​of each feature is expressed as follows: ; In the formula, For the full feature model, the first Predicted ice thickness for each sample This is the baseline predicted value, which is the mean of the predicted values ​​for all training samples. For the first The feature in the first SHAP values ​​on each sample. The value reflects the first The feature in the first The direction and magnitude of the contribution relative to the baseline predicted value in each sample. A value greater than 0 indicates that the feature causes the predicted value to be higher than the baseline. <0 indicates that the feature causes the predicted value to be lower than the baseline. =0 indicates that the feature does not contribute to the predicted value. For the specific implementation method of the TreeSHAP algorithm, those skilled in the art can refer to relevant literature, as it is a well-known technique in the field and will not be elaborated upon here.

[0041] S402, Calculate the global importance ranking index. For each feature, calculate the mean of its absolute values ​​of the SHAP values ​​across all samples, which serves as the global importance index for that feature. The calculation formula is as follows: ; In the formula, For the total number of samples (including the training and test sets, or using only the training set). For the first The feature in the first SHAP values ​​on each sample This indicates taking the absolute value. For the first The global importance value of each feature. A higher value indicates a more significant impact of that feature on the model output. (By...) The 26 features are sorted from largest to smallest to obtain a feature importance ranking table.

[0042] S403, based on the ranking results, performs factor hierarchical partitioning, dividing features into three levels: factors, auxiliary factors, and negligible factors. The partitioning criterion is: the first 5 features in the ranking table and... Features whose values ​​are more than twice that of the 6th ranked feature in the ranking list are defined as factors, and these factors play a decisive role in ice thickness prediction; features located between the 6th and 15th ranked features are considered as factors. Features whose values ​​are greater than the last value in the ranking list are defined as auxiliary factors. These auxiliary factors play a supplementary and moderating role in ice thickness prediction. Features from the 16th to the 26th position in the ranking list and... Features whose values ​​are less than 10% of the first value in the sorting table are defined as negligible factors. Negligible factors have a weak impact on ice thickness prediction and can be removed in subsequent steps.

[0043] S404 records the physical meaning of the factors. Based on the division results of S403, it typically includes 30-day cumulative freezing time x22, 14-day cumulative freezing time x21, water level elevation x3, 168-hour rolling minimum air temperature x13, and 72-hour rolling average ice surface temperature x17. These features correspond to cold accumulation parameters, hydrological condition parameters, and ice surface heat exchange parameters in ice thermodynamics theory, respectively. Their physical meanings are clear, and the subsequent step S6 will extract physical control parameters based on these factors. See attached document Figure 6 , Figure 6 This is a flowchart of feature selection and simplified forecast model construction according to an embodiment of the present invention. The following provides a detailed description of step S5, feature selection and simplified forecast model construction.

[0044] S501, based on the SHAP global importance ranking results obtained in step S4, performs Recursive Feature Elimination with Cross-Validation (RFECV) on all 26 features. RFECV is a feature selection method whose core process is as follows: starting from the full feature set, in each iteration, the feature with the lowest SHAP global importance in the current feature set is removed. The remaining features are used to retrain the random forest model, and the model performance is evaluated through 5-fold cross-validation. The model's coefficient of determination on the test set is recorded. Calculate the root mean square error (RMSE). Repeat the above iterative process until the feature set is empty. The calculation formula is: ; In the formula: This represents the number of validation samples for the current fold. For the first Measured ice thickness values ​​for a verification sample For the first Predicted ice thickness for one validation sample. This represents the mean of the measured ice thickness values ​​for the current verification sample. The formula for calculating RMSE is: ; Record the average of 5-fold cross-validation after each iteration. Value and average RMSE value. Select average. The optimal number of features is the number of features that corresponds to the point where the value reaches its maximum and the average RMSE value reaches its minimum. .Will The corresponding feature subset is denoted as the optimal feature subset Xopt.

[0045] S502, perform forward feature selection. Based on the SHAP global importance ranking from step S4, starting with the most important features, add one feature at a time. Retrain the random forest model using the current feature subset, evaluate the model performance on the test set, and record the test set results. Value. When a new feature is added, the test set... When the increment is less than the threshold ε=0.005, feature addition stops, and the feature subset at this point is denoted as the simplified feature subset Xsel. Test set Increment The calculation formula is: ; In the formula: For inclusion A subset of features of each feature For inclusion A subset of features of each feature To use feature subsets On the test set after training the model value, For the newly added After the test set has features The increment. When When the value is less than 0.005, the feature is considered to have limited contribution to improving the model's prediction accuracy and is therefore removed, while the previous feature is retained. Each feature is used as a subset of simplified features. .

[0046] S503, compare the optimal feature subset obtained in S501. The simplified feature subset obtained from S502 The intersection of the two is denoted as the final feature subset. .like and If the number of features is different, then choose the test set with fewer features. Not less than test set A subset of 0.95 times as .

[0047] S504, with the final feature subset Using the input as input and hourly ice thickness h as the target variable, the random forest regression model is retrained. The hyperparameter settings are consistent with S302, i.e., the number of trees n. estimators=500, maximum depth max_depth=None, minimum number of samples per internal node min_samples_split=5, minimum number of samples per leaf node min_samples_leaf=2, number of features considered in each split , among which | | represents the number of features in the final feature subset. After training, a simplified prediction model f is obtained. sel ( To evaluate the performance of the simplified forecast model using a test set, the test set is required. The RMSE and MAE values ​​must be ≥0.95 and not higher than 1.05 times the corresponding indices of the full-feature model. If the above conditions are met, the simplified forecast model will be used as the input model for subsequent steps S6 and S7; otherwise, return to S501, adjust the RFECV parameters or the forward feature selection threshold, and re-execute the feature selection process.

[0048] See attached document Figure 7 , Figure 7 This is a flowchart illustrating the extraction of physical control parameters based on the geometric features of an ALE curve according to an embodiment of the present invention. The following provides a detailed description of step S6, the extraction of physical control parameters based on the geometric features of an ALE curve.

[0049] S601, determine the factors for which ALE curves need to be generated. Based on the factor hierarchy division results in step S4, select factors as the ALE curve analysis objects. Factors include 30-day cumulative freezing time x22, 14-day cumulative freezing time x21, water level elevation x3, 168-hour rolling minimum air temperature x13, and 72-hour rolling average ice surface temperature x17. For each factor, perform operations S602 to S604 independently.

[0050] S602, calculate the one-dimensional Accumulated Local Effects (ALE) curve. The ALE curve describes the average local influence of a single feature on the model's predicted value, and can eliminate bias caused by correlation between features. For feature xj, its value range is divided into N equal intervals, where N ranges from 50 to 100. The boundary of the k-th interval is denoted as... -1 and Where k = 1, 2, ..., N, Features The minimum value, The maximum value of feature xj. ALE value. The calculation formula is: ; in, This is a range index, with values ​​from 1 to N; For the summation index, values ​​range from 1 to... ; For the first The number of samples within each interval; The sample index represents all features. Falling in each interval Samples within; For the first Features in each sample The possible values ​​of ; For the first In each sample, excluding features The value vectors of all other features; To feature Fixed as the upper bound of the interval Then, use the full feature model. For the Predicted values ​​for each sample; To feature Fixed as the lower bound of the interval Then, use the full feature model. For the Predicted values ​​for each sample; Features In taking values The ALE value at the location represents the feature. From the minimum value to The cumulative average impact of time on the model's predicted values.

[0051] S603, smoothing the ALE curve. To eliminate local noise in the ALE curve, a sliding window averaging method is used to smooth it. The sliding window width is set to 5 intervals, meaning that for each point... The average of the ALE values ​​from the two intervals before and after the given value (a total of five intervals) is taken as the smoothed ALE value. The smoothed ALT curve is denoted as... .

[0052] S604, calculate the first and second derivatives of the ALE curve. The first derivative reflects the trend of the ALE curve, and the second derivative reflects the turning points in the trend. The first derivative is calculated using the central difference method, and the formula is: ; In the formula: The distance between the center points of adjacent intervals. and They are respectively Smoothed ALE values ​​of the center points of adjacent intervals. Features In taking values The first derivative at point A. The second derivative is obtained by further difference of the first derivative, and the formula is: ; In the formula: and They are respectively The first derivative values ​​of the center points of adjacent intervals Features In taking values The second derivative at point .

[0053] S605, Identify geometric feature points of the ALE curve. Based on the calculation results of the first and second derivatives, identify the following three types of geometric feature points: points where the first derivative is zero, i.e., satisfying... The point where f = 0 corresponds to a local extremum of the ALE curve, representing the characteristic curve. At this value, the marginal effect on the model's predicted value is zero; the point where the second derivative is zero satisfies... The point where = 0 corresponds to an inflection point in the ALE curve's trend, representing a characteristic. The marginal impact on the model's predicted values ​​changes from increasing to decreasing or from decreasing to increasing; the peak point of the curve, i.e., the local maximum point on the ALE curve, is identified by comparing the ALE values ​​of adjacent points. If the ALE value of a point is greater than the ALE values ​​of its two adjacent points, it is determined to be a peak point. For the above geometric feature points, their precise location in the feature value space is obtained through linear interpolation.

[0054] S606, Establish a mapping rule base between physical control parameters and the geometric characteristics of ALE curves. The mapping rule base contains the following rules.

[0055] Rule A: Extraction of the critical threshold for freezing energy. For the 30-day cumulative freezing time feature x22, find the point on its ALE curve where the first derivative is zero and the ALE value changes from negative to positive. The feature value corresponding to this point is denoted as the critical threshold for freezing energy, Cfreeze. This threshold indicates that after the accumulated cold energy reaches this value, the ice thickness begins to increase significantly.

[0056] Rule B, freezing energy saturation threshold extraction: For the 30-day cumulative freezing time feature x22, find the point on its ALE curve where the second derivative, after the first derivative changes from rising to falling, is zero. The feature value corresponding to this point is denoted as the freezing energy saturation threshold Csat. This threshold indicates that after the accumulation of cold energy reaches this value, the rate of ice thickness growth begins to decline, and the promoting effect of cold energy accumulation on ice thickness growth tends to saturate.

[0057] Rule C, three-stage partition threshold extraction. Based on Cfreeze and Csat obtained from Rule A and Rule B, the slow growth start threshold Cslow is further calculated, and Cslow is 50% of Cfreeze. The three thresholds satisfy Cslow<Cfreeze<Csat, dividing the cooling capacity accumulation state into four intervals: a severe cooling capacity deficiency interval (0 to Cslow), a slow growth interval (Cslow to Cfreeze), a rapid thickening interval (Cfreeze to Csat), and a thickening saturation interval (greater than Csat).

[0058] Rule D, water level inflection point of inhibition and optimal hydrological window extraction. For the water level elevation feature x3, find the point on its ALE curve where the first derivative is zero and the ALE value reaches the maximum, and the corresponding feature value at this point is recorded as the water level inhibition inflection point Hsup. This threshold indicates that after the water level reaches this value, a further rise will exert an inhibitory effect on ice thickness growth. On the ALE curve, identify the water level interval where the ALE value is greater than 90% of the ALE value corresponding to Hsup, and this interval is recorded as the optimal hydrological window [Hlow, Hhigh].

[0059] Rule E, optimal ice-promoting air temperature and air temperature inhibition threshold extraction. For the 168-hour rolling minimum air temperature feature x13, find the point on its ALE curve where the first derivative is zero and the ALE value reaches the maximum, and the corresponding feature value at this point is recorded as the optimal ice-promoting air temperature Topt. This threshold indicates that when the rolling minimum air temperature reaches this value, the promoting effect on ice thickness growth is the strongest. On the ALE curve, find the two corresponding feature values when the ALE value drops from the maximum to 50% of the maximum, wherein the smaller value is recorded as the air temperature inhibition threshold Tsup, and the larger value is recorded as the upper air temperature threshold Tup.

[0060] Rule F, multi-scale cooling accumulation synergistic effect discrimination criterion establishment. For the 14-day cumulative freezing degree-hour feature x21 and the 30-day cumulative freezing degree-hour feature x22, calculate the integral values of their ALE curves within the respective feature value ranges, which are recorded as A21 and A22 respectively, and calculate the synergistic enhancement mode discrimination variable γsyn=A22 / A21. When γsyn>1, it indicates that the effect of long-term cooling accumulation on ice thickness growth is greater than that of short-term cooling accumulation, and a multi-scale synergistic enhancement effect exists. Calculate the single-scale dominant mode discrimination variable γdom=max(A21,A22) / min(A21,A22). When γdom>2, it indicates that the cooling accumulation of a single time scale is dominant. Calculate the thermal-hydrological hedging intensity coefficient γcancel, and its calculation formula is: ; In the formula: is the maximum absolute value of the ALE curve of the 30-day cumulative freezing degree-hour feature, is the maximum value of the absolute value of the ALE curve of the water level elevation feature. The value ranges from 0 to 1, and a smaller value indicates a greater hedging intensity between the cold accumulation effect and the water level effect.

[0061] S607: Summarize all the physical control parameters extracted in S606 to form a physical control parameter set. This set includes: thermodynamic parameters Cslow, Cfreeze, Csat, hydrological parameters Hsup, Hlow, Hhigh, air temperature parameters Topt, Tsup, Tup, and collaborative discrimination parameters γsyn, γdom, γcancel. The above parameters together constitute the input for generating the graded early warning rules in step S7 With reference to the accompanying Figure 8 , Figure 8 is a flowchart for generating graded early warning rules based on physical parameters according to an embodiment of the present invention. The step S7 of generating graded early warning rules based on physical parameters is described in detail below.

[0062] S701: Summarize the physical control parameter set extracted in step S6. This set includes thermodynamic parameters: slow growth initiation threshold Cslow, freezing energy critical threshold Cfreeze, freezing energy saturation threshold Csat; hydrological parameters: water level suppression turning point Hsup, lower limit of the optimal hydrological window Hlow, upper limit of the optimal hydrological window Hhigh; air temperature parameters: optimal ice-promoting air temperature Topt, air temperature suppression threshold Tsup, upper air temperature threshold Tup; collaborative discrimination parameters: collaborative enhancement mode discrimination variable γsyn, single-scale dominant mode discrimination variable γdom, thermal-hydrological hedging intensity coefficient γcancel. The above parameters constitute the input for generating the graded early warning rules.

[0063] S702: Generate graded monitoring and early warning rules for cold accumulation status. Based on the thermodynamic parameters, according to the 30-day cumulative freezing degree-hour at the current moment and the relationship with the thresholds, the cold accumulation status is divided into four grades, and each grade corresponds to a different monitoring frequency and early warning level. The status division conditions are: when <Cslow, it is determined as a state of severe cold deficiency, the monitoring frequency is maintained at 1 time per hour, and no early warning is issued; when Cslow≤ <Cfreeze, it is determined as a slow growth state, the monitoring frequency is increased to 1 time per 30 minutes, and a blue early warning is issued; when Cfreeze≤ <Csat, it is determined as a rapid thickening state, the monitoring frequency is increased to 1 time per 15 minutes, and a yellow early warning is issued; when ≥Csat, it is determined as a thickening saturation state, the monitoring frequency is reduced to 1 time per 2 hours, and a green early warning is issued, indicating that the cold accumulation is saturated and the ice thickness growth tends to stagnate.

[0064] S703 generates a joint early warning rule for rapid thickening. This rule requires both cold accumulation and water temperature conditions to be met simultaneously. The cold accumulation condition is: the cumulative freezing hours over 30 days at the current time are greater than or equal to the critical freezing energy threshold, i.e., CFDH. 30 ≥Cfreeze, and the change in the cumulative freezing time over the past 24 hours (30 days) is greater than zero, i.e. >0, where The calculation formula is: ; In the formula: The cumulative freeze time over the past 30 days at the current moment. This represents the cumulative freezing time over the past 30 days (24 hours prior). The water temperature condition is: the current 168-hour rolling minimum air temperature must be between the air temperature suppression threshold and the optimal ice-inducing air temperature, i.e., Tsup ≤ air_temp_min_168h ≤ Topt. When both conditions are met simultaneously, a joint early warning for rapid ice thickening is triggered, an orange alert is issued, and it is recommended to increase the ice thickness monitoring frequency to once every 10 minutes.

[0065] S704 generates high water level instability early warning and thermal-hydrological offsetting early warning rules. The high water level instability early warning condition is: the current water level elevation is greater than or equal to the water level suppression inflection point, i.e., H ≥ Hsup. When the condition is met, a high water level instability early warning is triggered, issuing a yellow warning, indicating that the high water level may lead to ice instability or rupture risk. The thermal-hydrological offsetting early warning condition is: the current 30-day cumulative freezing hours are greater than or equal to the freezing energy critical threshold, i.e. The conditions are: ≥Cfreeze, and the current water level elevation is greater than or equal to the upper limit of the optimal hydrological window (H≥Hhigh), and the thermal-hydrological offset intensity coefficient is greater than 0.3 (γcancel>0.3). When all three conditions are met simultaneously, a thermal-hydrological offset warning is triggered, and an orange warning is issued, indicating that the cold accumulation effect and the water level effect are offsetting each other, increasing the uncertainty of ice thickness growth prediction, and requiring comprehensive judgment based on field observation data.

[0066] S705 generates a multi-scale cold accumulation synergistic enhancement early warning rule. This rule is based on synergistic discriminant parameters. The conditions for multi-scale cold accumulation synergistic enhancement are: the discriminant variable of the synergistic enhancement mode is greater than 1, i.e., γsyn>1, and the discriminant variable of the single-scale dominant mode is between 0.5 and 2, i.e., 0.5≤γdom≤2. When both conditions are met simultaneously, a multi-scale cold accumulation synergistic enhancement effect is determined to exist, triggering a multi-scale cold accumulation synergistic enhancement early warning, issuing a red alert, and implementing the following two measures: First, the monitoring frequency is increased to once every 3 hours; second, the early warning lead time is extended by 12 hours, that is, the current warning status is extended forward by 12 hours.

[0067] S706 establishes a localization and migration mechanism for early warning thresholds. All thresholds used in S702 to S705 are automatically extracted from the data by the algorithm in step S6, without relying on manual presets. When this method is migrated to a new river section, steps S1 to S6 are repeated. The model is retrained using observation data from the new river section, and physical control parameters are extracted to obtain localized early warning thresholds suitable for the new river section. This mechanism ensures that the early warning rules can adapt to the differences in hydrothermal conditions in different river sections. See attached document Figure 9 , Figure 9 This is a flowchart of the model validation consistency verification process according to an embodiment of the present invention. The following provides a detailed description of step S8, the model validation consistency verification step.

[0068] S801, calculate the residuals of the simplified forecast model on the test set. The residuals are defined as the difference between the measured ice thickness and the model's predicted value, calculated using the following formula: ; In the formula: For the first Measured ice thickness values ​​for each test sample. For the first The simplified forecast model prediction values ​​for each test sample. For the first The residuals of each test sample, expressed in meters (m). The sign of the residual indicates the direction of the deviation between the model's predicted and measured values, while the absolute value of the residual indicates the magnitude of the deviation.

[0069] S802, perform distribution analysis on the residual sequence. Plot a residual histogram, dividing the residual range into 20 equal intervals, and calculate the proportion of samples in each interval to the total sample size, forming a residual distribution histogram. Plot a normal probability plot of the residuals (QQ plot), comparing the quantiles of the residuals with the theoretical quantiles of the standard normal distribution to graphically test whether the residuals follow a normal distribution. Calculate the mean of the residual sequence. and standard deviation The formula for calculating the mean is: ; In the formula: The total number of test samples, The mean of the residuals is the value; the closer the mean is to zero, the less systematic the model's predictions are. The standard deviation is calculated using the following formula: ; In the formula: This represents the standard deviation of the residuals; the smaller the value, the smaller the fluctuation in the model's predictions.

[0070] S803, perform an autocorrelation test on the residual series. The Durbin-Watson statistic is used to test for the existence of first-order autocorrelation in the residual series. The calculation formula is: ; In the formula: For the first The residuals of each test sample For the first The residuals of each test sample This is the Durbin-Watson statistic, with values ​​ranging from 0 to 4. When ∈[1.5,2.5], it is determined that the residual sequence has no significant first-order autocorrelation, and the model residuals are independent; when When <1.5, the residuals are considered to have positive autocorrelation; when When the value is greater than 2.5, the residuals are considered to have negative autocorrelation. If the residuals have significant autocorrelation, it indicates that the model may have missed time series structure information, and it is necessary to return to step S5 to re-select features or adjust the model structure.

[0071] S804, Statistical analysis of ice thickness error is performed based on the measured ice thickness range. The measured ice thickness values ​​in the test set are... The range of values ​​is divided into multiple intervals. The interval division standard is as follows: with a step size of 0.1m, the range is divided from 0m to the maximum ice thickness value, and each interval contains an ice thickness range of 0.1m. For each ice thickness interval, the number of samples falling into that interval is counted, and the root mean square error of the samples within that interval is calculated using the following formula: ; In the formula: For ice thickness interval index, For the first A set of indices for samples within a given ice thickness range. For the first Number of samples within each ice thickness range For the first Measured ice thickness for each sample For the first Predicted ice thickness for each sample For the first The root mean square error (RMSE) within each ice thickness interval is expressed in meters (m). The mean absolute error (MAE) of the samples within each interval is also calculated using the following formula: ; In the formula: For the first The average absolute error within a given ice thickness interval is expressed in meters (m). or is significantly higher than other intervals, it indicates that the model has a systematic deviation within this ice thickness range, and it is necessary to analyze the cause of the deviation and carry out targeted optimization.

[0072] S805: Perform model verification consistency check. The model verification consistency check includes two aspects. First, verify whether the physical control parameters extracted in step S6 satisfy the inherent constraint relationships of physical laws. The inequality constraints that need to be checked include: the slow growth start threshold is less than the critical freezing energy threshold, that is, Cslow<Cfreeze; the critical freezing energy threshold is less than the freezing energy saturation threshold, that is, Cfreeze<Csat; the lower limit of the optimal hydrological window is less than the upper limit of the optimal hydrological window, that is, Hlow<Hhigh; the air temperature suppression threshold is less than the optimal ice-promoting air temperature, that is, Tsup<Topt; the optimal ice-promoting air temperature is less than the upper air temperature threshold, that is, Topt<Tup. If any of the above constraint conditions is not satisfied, it is determined that the model verification consistency check fails, and it is necessary to return to step S6 to check whether there is an error in the identification process of the geometric feature points of the ALE curve. Second, verify whether the correspondence between the geometric features of the ALE curve and the physical process is reasonable. Specifically: the critical freezing energy threshold Cfreeze corresponds to the point where the ALE curve changes from negative to positive, its physical meaning is that after the cold accumulation reaches this value, the ice thickness begins to increase significantly, and this point should be located in the medium-low level interval of cold accumulation; the freezing energy saturation threshold Csat corresponds to the point where the first derivative of the ALE curve changes from rising to falling, its physical meaning is that after the cold accumulation reaches this value, the ice thickness growth tends to be saturated, and this point should be located in the high level interval of cold accumulation; the water level suppression turning point Hsup corresponds to the point where the ALE curve reaches the maximum value, its physical meaning is that after the water level exceeds this value, the promoting effect on ice thickness growth turns into a suppressing effect, and this point should be located in the medium-high interval of the water level. If the positions of the above geometric feature points obviously conflict with the physical expectations, it is determined that the model verification consistency check fails, and it is necessary to return to step S6 to adjust the number N of ALE curve interval divisions or the width of the smoothing window.

[0073] S806: Summarize the check results from S801 to S805 to generate a model verification report. If the residual analysis, autocorrelation test, interval error statistics and model verification consistency check all pass, it is determined that the streamlined prediction model meets the application requirements, and the model and early warning rules can be deployed to the actual monitoring system. If any check fails, return to the corresponding step for adjustment according to the check result, and re-execute the model training and verification process until all checks pass.

[0074] Although the embodiments of the present invention have been shown and described, it will be understood by those of ordinary skill in the art that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principle and spirit of the present invention, and the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics, characterized in that, Includes the following steps: S1, Multi-parameter continuous monitoring data acquisition: Collect hourly air temperature, hourly ice surface temperature, hourly water level elevation and hourly ice thickness data of the target river section; S2, Construction of a multi-scale hydrothermal feature system: Based on the hourly air temperature, hourly ice surface temperature and hourly water level elevation, a multi-scale hydrothermal feature matrix is ​​constructed. The feature matrix includes real-time observation features, rolling statistical features of air temperature, rolling average features of ice surface temperature, cumulative freezing time features across multiple time scales, and temperature change rate features. S3, Random Forest Full Feature Model Training: Using the feature matrix as input and hourly ice thickness as the target variable, train the random forest regression model to obtain the full feature ice thickness prediction model; S4, Global importance ranking and factor decoupling based on SHAP: Calculate the SHAP value of all features in the full feature model, sort them from largest to smallest according to the mean absolute value of SHAP, and divide the features into factors, auxiliary factors and negligible factors; S5, Feature selection and simplified prediction model construction: Based on the ranking results, the optimal feature subset is determined by recursive feature elimination cross-validation and forward feature selection. The random forest regression model is then retrained with the optimal feature subset to obtain the simplified prediction model. S6, Physical control parameter extraction based on ALE curve geometric features: A one-dimensional ALE curve is generated for the factor, and physical control parameters are extracted by identifying the points where the first derivative of the ALE curve is zero, the points where the second derivative is zero, and the peak points of the curve. S7, Generation of hierarchical early warning rules based on physical parameters: Based on the physical control parameters, establish hierarchical monitoring and early warning rules for cold accumulation status, joint early warning rules for rapid thickening, early warning rules for high water level instability and thermo-hydrological offsetting, and multi-scale cold accumulation synergistic enhancement early warning rules. S8, Model validation consistency check: Perform residual distribution analysis, residual autocorrelation test, ice thickness interval error statistics and model validation consistency check on the simplified forecast model.

2. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S2, the real-time observation features include real-time air temperature, real-time ice surface temperature, and real-time water level elevation; the rolling statistical features of air temperature include rolling averages and rolling minimums with window lengths of 6 hours, 12 hours, 24 hours, 72 hours, and 168 hours; the rolling average features of ice surface temperature include rolling averages with window lengths of 6 hours, 12 hours, 24 hours, and 72 hours; the cumulative freezing time features across multiple time scales include cumulative freezing time with time scales of 1 day, 3 days, 7 days, 14 days, and 30 days; and the temperature change rate features include the 1-hour change rate of air temperature, the 24-hour change rate of air temperature, the 1-hour change rate of ice surface temperature, and the 24-hour change rate of ice surface temperature.

3. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S4, the factor is the feature in the first 5 positions of the sorting table whose average absolute SHAP value is more than twice the value of the 6th position of the sorting table; the auxiliary factor is the feature in the 6th to 15th positions of the sorting table whose average absolute SHAP value is greater than the value of the last position of the sorting table; the negligible factor is the feature in the 16th to 26th positions of the sorting table whose average absolute SHAP value is less than 10% of the value of the 1st position of the sorting table.

4. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S5, recursive feature elimination cross-validation is used to determine the optimal number of features, forward feature selection is used to determine a simplified feature subset, and the intersection of the feature subset corresponding to the optimal number of features and the simplified feature subset is taken as the final feature subset. The random forest regression model is retrained with the final feature subset to obtain a simplified prediction model.

5. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S6, the physical control parameters include: the critical threshold for freezing energy, the threshold for freezing energy saturation, and the threshold for slow growth initiation calculated from the 30-day cumulative freezing time feature; the water level suppression inflection point, the lower limit of the optimal hydrological window, and the upper limit of the optimal hydrological window extracted from the water level elevation feature; the optimal ice-promoting temperature, the temperature suppression threshold, and the upper limit of the temperature extracted from the 168-hour rolling minimum temperature feature; and the discriminant variables for the synergistic enhancement mode, the discriminant variables for the single-scale dominant mode, and the thermal-hydrological offsetting intensity coefficient calculated based on the ALE curves of each factor.

6. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S7, the cold accumulation state classification monitoring and early warning rule is as follows: based on the comparison of the current 30-day cumulative freezing time with the slow growth initiation threshold, the freezing energy critical threshold, and the freezing energy saturation threshold, the cold accumulation state is divided into a severely insufficient cold energy state, a slow growth state, a rapid thickening state, and a thickening saturation state, and each corresponds to a different monitoring frequency and early warning color.

7. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S7, the rapid thickening joint early warning rule is as follows: when the cumulative freezing time over 30 days at the current moment is greater than or equal to the freezing energy critical threshold, and the change in the cumulative freezing time over 30 days in the past 24 hours is greater than zero, and the minimum rolling temperature over 168 hours at the current moment is between the temperature suppression threshold and the optimal ice-promoting temperature, the rapid thickening joint early warning is triggered.

8. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S7, the rules for high water level instability warning and thermal-hydrological offset warning include: high water level instability warning, which is triggered when the current water level elevation is greater than or equal to the water level suppression inflection point; thermal-hydrological offset warning, which is triggered when the current 30-day cumulative freezing time is greater than or equal to the freezing energy critical threshold, the current water level elevation is greater than or equal to the upper limit of the optimal hydrological window, and the thermal-hydrological offset intensity coefficient is greater than 0.

3.

9. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S7, the multi-scale cold accumulation synergistic enhancement early warning rule is as follows: when the discriminant variable of the synergistic enhancement mode is greater than 1 and the discriminant variable of the single-scale dominant mode is between 0.5 and 2, it is determined that there is a multi-scale cold accumulation synergistic enhancement effect, triggering a red warning, and increasing the monitoring frequency to once / 3 hours, and extending the warning lead time by 12 hours.

10. The method for interpretable prediction and early warning of river ice thickness based on multi-scale hydrothermal characteristics according to claim 1, characterized in that, In step S8, the model validation consistency check includes verifying whether the following inequality constraints hold: the slow growth initiation threshold is less than the freezing energy critical threshold, the freezing energy critical threshold is less than the freezing energy saturation threshold, the lower limit of the optimal hydrological window is less than the upper limit of the optimal hydrological window, the temperature suppression threshold is less than the optimal ice-promoting temperature, and the optimal ice-promoting temperature is less than the upper limit of the temperature threshold.

Images