Visual post-hoc interpretation method for integrated hydrometeorological machine learning forecasting model

Abstract

The present invention discloses a visualized post-explanation method for an integrated hydro-meteorological machine learning forecasting model, the method comprising the following steps: establishing a single-site historical database and a forecasting factor database, and performing feature selection; establishing a hydro-meteorological factor forecasting black box model, and statistically calculating forecast errors; the post-explanation method explains the machine learning forecasting model from multiple perspectives, evaluates the reliability of the model results based on the feature importance of MSE and the Shapley method, performs interpretability analysis of feature effects based on partial dependency graphs, individual conditional expectations, and cumulative local effects, and is used to diagnose the model. For the forecast of a single-day hydro-meteorological element, find the features that are abnormal in the single-day forecast through individual expectation conditions, and judge whether they are feature outliers based on historical data. The present invention evaluates the feature effects of the hydro-meteorological machine learning model based on the post-explanatory method, and improves the reliability and robustness of the model.

Description

Technical Field

[0001] The present invention relates to the technical field of hydrological and meteorological forecasting, and in particular to a visual post-interpretation method of an integrated hydrological and meteorological machine learning forecasting model. Background Art

[0002] The concept of hydrometeorological forecasting is to provide an estimate of future basin conditions based on recent and predicted meteorological conditions and hydrological observations. General forecasts include precipitation, evapotranspiration, runoff, soil moisture, etc. These factors interact and aggregate to form the natural cycle of water. Potential application areas of hydrometeorological forecasting include: forecasting of extreme events, such as flood and drought warnings; water scheduling, hydropower operations, long-term water resources planning, precision agricultural irrigation and watershed management; climate change impact assessment, etc.

[0003] Due to the highly nonlinear nature of the atmosphere, the prediction of hydrological and meteorological factors is of course difficult. In order to achieve accurate hydrological and meteorological forecasts, people have done a lot of work in both model improvement and data-driven. Physics-based models are based on the hydrological process mechanism of the water cycle and are transparent models with high interpretability, such as the Weather Research and Forecasting (WRF) model. However, the relatively low accuracy of physics-based models is the biggest obstacle to hydrological and meteorological forecasts.

[0004] With the booming development of artificial intelligence in recent decades, machine learning models can effectively improve forecast accuracy. Some scholars have used deep learning convolutional neural networks and long short-term memory to improve local precipitation forecasts to the 3-hour level. At the same time, deep learning also has a good response to extreme precipitation. Both machine learning and deep learning can be applied to evapotranspiration forecasts to improve the accuracy of short-term forecasts. However, machine learning is often considered a black box. Although these artificial intelligence methods provide more accurate predictions of hydrological and meteorological variables, their lack of explanation due to the black box characteristics constitutes a practical challenge in building data-driven models. Explainable Artificial Intelligence (XAI) aims at the communication between artificial intelligence models and human understanding, and attempts to build trust in opaque models. The lack of explainability hinders the further application of black box models due to the need to explain the decisions made by the model.

[0005] In the field of practical hydrometeorological forecasting, lack of explanations can cause the following obstacles:

[0006] (1) Since the black box model lacks prior knowledge of hydrological and meteorological processes, the results are questionable.

[0007] (2) Conventional black-box models only provide hydrological and meteorological forecast results, but do not provide the supporting conditions for the forecast results.

[0008] (3) The inherent bias of the training data may cause defects in the black box model, resulting in large deviations in the model prediction results.

[0009] (4) Hydrological and meteorological factors can cause different models to perform differently at different times and locations.

[0010] In order to overcome the interpretability challenge of black box models, many techniques have been developed, including intrinsically interpretable models and ex post interpretation methods. After training with a large amount of data, black box models can be considered to be more optimized forecasting models. Ex post interpretation methods are expected to overcome the interpretability challenge of black box models. These methods can be divided into two parts: importance evaluation and feature effect visualization. Methods for evaluating feature importance include permutation feature importance algorithm, impure importance of tree model, feature method based on mean square error, Shapley method, sequential backward selection or sequential forward selection, etc. In addition, many tools provide visualization of feature effects, such as partial dependence plots, individual conditional expectations, and cumulative local effects. Some scholars have applied sequential backward selection, sequential forward selection, and partial dependence plots to the forecast of meteorological factors at different spatiotemporal scales, achieving the purpose of interpretability of black box models. However, none of the methods introduced above can fully explain the prediction model. Summary of the invention

[0011] In order to overcome the defects and shortcomings of the prior art, the present invention provides a visual post-interpretation method for an integrated hydrological and meteorological machine learning forecasting model, which evaluates the characteristic effects of the hydrological and meteorological machine learning model based on a post-interpretable method, thereby improving the reliability and robustness of the model.

[0012] In order to achieve the above object, the present invention adopts the following technical solutions:

[0013] The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model includes the following steps:

[0014] Establish a single-site historical database and a predictor database, and perform feature selection;

[0015] Establish a black box model for hydrological and meteorological factor forecasting and calculate forecast errors;

[0016] Post hoc interpretation methods explain machine learning forecasting models from multiple perspectives, including feature importance based on MSE, Shapley's method, partial dependence plots, individual conditional expectations, and cumulative local effects.

[0017] The feature importance based on MSE is explained as follows: multiple decision trees are formed based on the training data, out-of-bag data are selected to test the decision trees, and the prediction mean square error (MSE) of the prediction model is obtained. The feature values ​​of a single feature of the out-of-bag data are randomly arranged to obtain a new prediction error MSE, and the difference in the prediction error MSE before and after the arrangement is calculated. The difference in the prediction error MSE of each decision tree is averaged, and the features are ranked in importance from the prediction model accuracy according to the average value;

[0018] The explanation is based on Shapley's method, which is as follows: the contribution of each feature to the predicted value is converted into the proportion of the predicted value, and the contribution to the predicted value is weighted summed over all possible feature combinations to obtain the importance of each feature;

[0019] The explanation is based on the partial dependence diagram, specifically: select the focus feature in the time series training data, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the average value of the forecast value obtained by each replacement, which is the partial dependence value, which is used to reflect the average change of the local hydrological and meteorological forecast object when it is affected by a certain factor;

[0020] The explanation is based on individual conditional expectation, specifically: select the focus feature in the training data at a single time point, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the forecast value obtained by each replacement, which is the partial dependence value at a single time point, used to reflect the impact of the feature on the result in the forecast of a single instance;

[0021] The explanation is based on the cumulative local effect. Specifically, the attention feature in the training data at a single time point is selected, and the training data is divided into two subsets based on the attention feature to generate a distribution grid of the attention feature. In a single grid, the data value of the attention feature distribution grid is used to replace the attention feature one by one, and the forecast value obtained by each replacement is calculated. The differences in the forecast values ​​in the grid are averaged and then accumulated on the grid to reflect the feature effect.

[0022] The reliability of the model results is evaluated based on the feature importance of MSE and the Shapley method. The interpretability analysis of the feature effects is performed based on the partial dependence diagram, individual conditional expectation and cumulative local effect. This is used to diagnose the model. For the forecast of single-day hydrological and meteorological elements, the abnormal features in the single-day forecast are found through individual expectation conditions, and whether they are abnormal values ​​of the features is determined based on historical data.

[0023] As a preferred technical solution, a single-site historical database and a forecast factor database are established, and the data in the database include time series, meteorological forcing, radiation flux, energy flux and hysteresis;

[0024] Feature selection methods include empirical selection and sequential backward selection.

[0025] As a preferred technical solution, the hydrological and meteorological factor forecasting black box model adopts any one of a random forest model, a naive Bayesian regression model, a support vector machine model or a deep learning model.

[0026] As a preferred technical solution, the interpretation based on the Shapley method includes the following specific steps:

[0027] Get the training set, extract the features of interest, and calculate the SV value of the features of interest, expressed as:

[0028]

[0029] Among them, p represents the number of features, x k Characteristics that indicate concern, represents the SV value of the kth feature of interest, S represents the set of all feature combinations, S i represents the i-th feature combination in the S set, Represents the trained forecast model.

[0030] As a preferred technical solution, the interpretation based on individual conditional expectations includes the following specific steps:

[0031] Select the feature of interest x in the training data k , based on the focus feature x k The training data is divided into two subsets, namely S = {x k} and C={x1,x2,…,x k-1 ,x k+1 …,x p};

[0032] We will focus on feature x k The continuous set of possible values ​​is x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m};

[0033] Compute the set of possible values ​​x s,j The partial dependence value of We will focus on feature x k By x s,j replace, Fixed, obtain n forecast values ​​and average them, and then calculate the average value of the forecast values ​​obtained each time, which is the partial dependence value The specific calculation formula is expressed as:

[0034]

[0035] Among them, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of Except x k All features except , n is the number of instances in the time series, i is the i-th data in the time series, It is a forecast model;

[0036] The partially dependent values ​​will be calculated By continuous x s Perform a visual presentation to show the focus feature x k Changes in the forecast value when the change occurs;

[0037] The feature is replaced iteratively and repeatedly to obtain the change of the predicted value under the change of the feature, so as to obtain whether the influence of the feature on the predicted value is linear or nonlinear.

[0038] As a preferred technical solution, the interpretation based on individual conditional expectations includes the following specific steps:

[0039] Select a data set P at a single time point = {x1, x2, ..., x k ,…,x p}, p is the number of features, x k is the feature of interest, x k P is divided into two subsets, namely S = {x k} and C={x1,x2,…,x k-1 ,x k+1 …,x p}; Generate x k The set of possible values ​​x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m}, m is x k The number of possible values;

[0040] Use x s The possible values ​​replace the feature x one by one k , calculate the forecast value obtained for each substitution, which is the partial dependence value at a single time point

[0041]

[0042] Among them, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of x c Except x k All features except is the forecast model, and obtains the continuous set x s Corresponding Value, in x s The element is the horizontal axis, The value is the vertical axis.

[0043] As a preferred technical solution, the interpretation based on cumulative local effects includes the following specific steps:

[0044] Extract the two relevant features of interest, make their possible values ​​continuous and gridded to form a two-dimensional grid;

[0045] The vertex values ​​of the grid are used to replace the feature values ​​of the data points in the grid, while other features remain unchanged, and are input into the black box of the machine learning model for prediction;

[0046] The four forecast values ​​obtained are subtracted to form second-order differences, each second-order difference is accumulated and centered to obtain the cumulative local effect, which is then plotted into a cumulative local effect diagram.

[0047] As a preferred technical solution, one-dimensional ALE or two-dimensional ALE is used to reflect the impact of feature changes on the forecast results.

[0048] A visual post-explanation system integrating hydrological and meteorological machine learning forecasting models, including: database construction module, feature selection module, model construction module, forecast error statistics module and explanation module;

[0049] The database construction module is used to establish a single-site historical database and a prediction factor database;

[0050] The feature selection module is used to perform feature selection;

[0051] The model building module is used to establish a black box model for hydrological and meteorological factor forecasting;

[0052] The forecast error statistics module is used to count the forecast errors;

[0053] The explanation module is used to explain the machine learning prediction model from multiple perspectives using post-hoc explanation methods, including feature importance based on MSE, Shapley method, partial dependence plot, individual conditional expectation and cumulative local effect;

[0054] The feature importance based on MSE is explained as follows: multiple decision trees are formed based on the training data, out-of-bag data are selected to test the decision trees, the prediction error MSE of the prediction model is obtained, the feature values ​​of a single feature of the out-of-bag data are randomly arranged, a new prediction error MSE is obtained, and the difference in the prediction error MSE before and after the arrangement is calculated, the prediction error MSE difference of each decision tree is averaged, and the features are ranked in importance from the prediction model accuracy according to the average value;

[0055] The explanation is based on Shapley's method, which is as follows: the contribution of each feature to the predicted value is converted into the proportion of the predicted value, and the contribution to the predicted value is weighted summed over all possible feature combinations to obtain the importance of each feature;

[0056] The explanation is based on the partial dependence diagram, specifically: select the focus feature in the time series training data, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the average value of the forecast value obtained by each replacement, which is the partial dependence value, which is used to reflect the average change of the local hydrological and meteorological forecast object when it is affected by a certain factor;

[0057] The explanation is based on individual conditional expectation, specifically: based on the feature of interest, the training data is divided into two subsets, a set of possible values ​​of the feature of interest is generated, the data values ​​in the set of possible values ​​of the feature of interest are used to replace the feature of interest one by one, and the forecast value obtained by each replacement is calculated, which is the partial dependence value at a single time point, which is used to reflect the influence of the feature on the result in the forecast of a single instance;

[0058] The explanation is based on the cumulative local effect. Specifically, the attention feature in the training data at a single time point is selected, and the training data is divided into two subsets based on the attention feature to generate a distribution grid of the attention feature. In a single grid, the data value of the attention feature distribution grid is used to replace the attention feature one by one, and the forecast value obtained by each replacement is calculated. The differences in the forecast values ​​in the grid are averaged and then accumulated on the grid, thereby reflecting the impact of feature changes on the forecast results.

[0059] The reliability of the model results is evaluated based on the feature importance of MSE and the Shapley method. The interpretability analysis of the feature effects is performed based on the partial dependence diagram, individual conditional expectation and cumulative local effect. This is used to diagnose the model. For the forecast of single-day hydrological and meteorological elements, the abnormal features in the single-day forecast are found through individual expectation conditions, and whether they are abnormal values ​​of the features is determined based on historical data.

[0060] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0061] (1) To address the “black box” problem of machine learning forecasting in hydrology and meteorology, the present invention adopts a technical solution that integrates several visualization post-interpretation methods into a complete post-interpretation system, which solves the problem of unclear correlation between the internal input and output of the machine learning model, achieves the technical effect of enhancing the interpretability of the machine learning model, and further promotes the application of machine learning models in the field of hydrology and meteorology.

[0062] (2) The post-interpretable system can reflect the relationship between input features and output forecasts from multiple perspectives, provide evidence from multiple aspects, and provide support conditions for forecast results, including: 1) Feature importance ranking: MSE-based feature importance and Shapley method provide global and local (single forecast) model result credibility assessment for hydrological and meteorological element forecasts; 2) Feature effect: One-dimensional and two-dimensional partial dependence diagrams and cumulative local effects can respectively reflect the average change of forecast values ​​when a single or two features change, thereby reflecting the feature effect from a global perspective. Individual expectation conditions reflect the feature effect from a local perspective. Therefore, based on the feature effect, the inherent bias of the data or the model's own data recognition bias can be evaluated based on the local hydrological and meteorological physical process, and the model can be further debugged or feature selection can be re-performed.

[0063] (3) The present invention provides a set of open-source, integrated post-interpretability tools, which provide a reliable basis for studying the internal structure of the forecast model and reflecting the inherent bias of the data, and can provide visual method guidance for model debugging, model optimization and reflecting data characteristics. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 A schematic flow chart of a post-interpretation method for an integrated hydrological and meteorological machine learning forecasting model of the present invention;

[0065] Figure 2 A schematic diagram of a post-interpretation method for an integrated soil moisture machine learning forecasting model of the present invention;

[0066] Figure 3 The optimal number of features based on sequence backward selection in the present invention;

[0067] Figure 4 (a) shows the distribution diagram of the predicted values ​​and observed values ​​in the test set;

[0068] Figure 4 (b) is a schematic diagram showing the time series comparison of the predicted values ​​and the observed values ​​in the test set;

[0069] Figure 5 It is a schematic diagram of feature ranking based on feature importance of MSE in the present invention;

[0070] Figure 6This is a schematic diagram of the change of Shapley value of the characteristic calculated based on the Shapley method of the present invention over time;

[0071] Figure 7 (a) shows the distribution diagram of Shapley value of soil moisture characteristics with 1-5 days lag (SM1, SM2, SM3, SM4, SM5) along with characteristic value;

[0072] Figure 7 (b) shows the distribution diagram of Shapley value of rainfall characteristics on the day (P0) and rainfall lagged 1-2 days (P1, P2) along with the characteristic value;

[0073] Figure 7 (c) shows the distribution diagram of Shapley values ​​of eight characteristics including day order (DAY), day order of the year (DOY), air temperature (TA), incident shortwave radiation (PISW), soil temperature (TS), latent heat flux (LE), gross primary production (GPP) and ecosystem respiration efficiency (RECO) along with the characteristic value;

[0074] Figure 8 (a) is a local cumulative effect diagram of the present invention in which the input feature is the order of days (DAY);

[0075] Figure 8 (b) is a local cumulative effect diagram of the present invention where the input feature is the order of days of the year (DOY);

[0076] Figure 8 (c) is a local cumulative effect diagram of the present invention where the input feature is temperature (TA);

[0077] Figure 8 (d) is a local cumulative effect diagram of the present invention where the input characteristic is incident shortwave radiation (PISW);

[0078] Figure 8 (e) is a local cumulative effect diagram of the present invention where the input feature is soil temperature (TS);

[0079] Figure 8 (f) is a local cumulative effect diagram of the present invention where the input characteristic is latent heat flux (LE);

[0080] Figure 8 (g) is a local cumulative effect diagram of the present invention where the input characteristic is gross primary production (GPP);

[0081] Figure 8 (h) is a local cumulative effect diagram of the input feature of the present invention, which is the ecosystem respiration efficiency (RECO);

[0082] Figure 8(i) is a local cumulative effect diagram of the present invention where the input feature is soil moisture lagged by 1 day (SM1);

[0083] Figure 8 (j) is a local cumulative effect diagram of the present invention where the input feature is soil moisture lagged by 2 days (SM2);

[0084] Figure 8 (k) is a local cumulative effect diagram of the present invention where the input feature is soil moisture lagged by 3 days (SM3);

[0085] Figure 8 (1) is a local cumulative effect diagram of the present invention where the input feature is soil moisture lagged by 4 days (SM4);

[0086] Figure 8 (m) is a local cumulative effect diagram of the present invention where the input feature is soil moisture lagged for 5 days (SM5);

[0087] Figure 8 (n) is a local cumulative effect diagram of the present invention whose input feature is the rainfall on that day (P0);

[0088] Figure 8 (o) is a local cumulative effect diagram of the present invention with the input feature of rainfall lag of 1 day (P1);

[0089] Figure 8 (p) is a local cumulative effect diagram of the present invention with the input feature of rainfall lag 2 days (P2);

[0090] Fig. 9 (a) is an individual expected condition diagram of the present invention whose input feature is the order of days (DAY);

[0091] Fig. 9 (b) is an individual expected condition diagram of the present invention whose input feature is the order of days of the year (DOY);

[0092] Fig. 9 (c) is an individual expected condition diagram of the present invention whose input feature is temperature (TA);

[0093] Fig. 9 (d) is an individual expected condition diagram of the present invention where the input characteristic is incident shortwave radiation (PISW);

[0094] Fig. 9 (e) is an individual expected condition diagram of the present invention whose input feature is soil temperature (TS);

[0095] Fig. 9 (f) is an individual expected condition diagram of the present invention whose input characteristic is latent heat flux (LE);

[0096] Fig. 9(g) is an individual expected condition diagram of the present invention whose input feature is gross primary production (GPP);

[0097] Fig. 9 (h) is an individual expected condition diagram of the present invention whose input feature is ecosystem respiration efficiency (RECO);

[0098] Fig. 9 (i) is the individual expected condition diagram of the present invention whose input feature is soil moisture lagged by 1 day (SM1);

[0099] Fig. 9 (j) is the individual expected condition diagram of the present invention whose input feature is soil moisture lagged by 2 days (SM2);

[0100] Fig. 9 (k) is the individual expected condition diagram of the present invention whose input feature is soil moisture lagged by 3 days (SM3);

[0101] Fig. 9 (1) is the individual expected condition diagram of the present invention whose input feature is soil moisture lagged by 4 days (SM4);

[0102] Fig. 9 (m) is the individual expected condition diagram of the present invention whose input feature is soil moisture lagged by 5 days (SM5);

[0103] Fig. 9 (n) is the individual expected condition diagram of the present invention whose input feature is the rainfall on that day (P0);

[0104] Fig. 9 (o) is the individual expected condition diagram of the present invention with the input feature of rainfall lag of 1 day (P1);

[0105] Fig. 9 (p) is the individual expected condition diagram of the present invention with the input feature of rainfall lag 2 days (P2);

[0106] Figure 10(a) is a schematic diagram of the two-dimensional cumulative local effect under the interaction of two sets of characteristics: soil temperature (TS) and shortwave incidence (PISW);

[0107] Figure 10(b) is a schematic diagram of the two-dimensional cumulative local effect under the interaction of two sets of characteristics: soil temperature (TS) and latent heat (LE). DETAILED DESCRIPTION

[0108] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0109] Example

[0110] like Figure 1 As shown, the visual post-interpretation method of the integrated hydrological and meteorological machine learning forecasting model provided in this embodiment mainly includes the following steps:

[0111] S1: Establish a single-site historical database and a prediction factor database, and perform feature selection;

[0112] S2: Establish a black box model for hydrological and meteorological factor forecasting and calculate the forecast error (R 2 , MAE and RMSE);

[0113] S3: Post hoc interpretation methods explain machine learning forecasting models from multiple perspectives, including feature importance based on MSE, Shapley's method, partial dependence plots, individual conditional expectations, and cumulative local effects;

[0114] S4: Develop an integrated system based on post-hoc interpretation methods to explain machine learning models for hydrometeorological forecasts.

[0115] Hydrometeorological forecast is to forecast the state and possible impact of hydrometeorological elements such as precipitation, evaporation, soil moisture, and runoff in a water body, a watershed, or a station / grid in the future hydrological cycle. The forecast of hydrometeorological elements at a certain station (local area) has certain commonalities, which are reflected in 1) time correlation, runoff and soil moisture have strong lag; 2) being greatly affected by meteorological forcing, changes in rainfall, sunshine hours, and average temperature will have a great impact on hydrometeorological forecasts; 3) the embodiment of water and heat exchange on the land surface, the changes between hydrometeorological elements are basically reflected in the interaction of water and heat energy flux between the land surface and the atmosphere. Among them, since the soil surface has the function of storing water in the process of land-atmosphere interaction, soil moisture is crucial to the study of hydrometeorology, and the forecast of soil moisture can reflect the changes of other hydrological elements to a certain extent. Previous studies have shown that in the simulation of land surface models, soil moisture has a great influence on precipitation, runoff and evaporation.

[0116] The hydro-meteorological factor forecast provided in this embodiment takes soil moisture as the object of hydrological forecast. First, soil moisture is an indispensable variable in the global water cycle, which determines the management of water and soil resources, surface energy and the natural cycle of water; second, soil moisture is very sensitive to many factors, such as soil properties (sand and clay content), vegetation, and especially meteorological factors, including precipitation, air temperature, soil temperature, etc. In particular, the process of soil moisture formation is a highly nonlinear system, which makes soil moisture difficult to predict.

[0117] Therefore, this embodiment selects soil moisture as the hydrological and meteorological forecast object, but is not limited to the application of soil moisture forecast.

[0118] like Figure 2 As shown, the visual post-interpretation method of the integrated hydrological and meteorological machine learning forecasting model disclosed in this embodiment specifically includes the following steps:

[0119] S1: Establish a single-site historical database and a prediction factor database, and perform feature selection: The data are all from the observation data of the FLUXNET station, mainly including time series, meteorological forcing, radiation flux, energy flux and lag. The feature selection methods include empirical selection and sequence backward selection.

[0120] In this embodiment, the single site provided is from the Haibei National Field Observation and Research Station in Qinghai Province, China, located in the northeast of the Qinghai-Tibet Plateau (37°37′N, 101°19′E), codenamed Ha2, with an average altitude of 3200 meters, and the ecosystem surface classification is permanent wetland, the vegetation is shrub forest, and the soil is clay loam, with strong permeability, strong water retention capacity, and low thermal conductivity. The region is extremely sensitive to the effects of climate change, is a sensitive area of ​​global climate change, and is the "starting area" of climate change in China. Studying the interpretability of soil moisture machine learning forecasts is of great significance to the study of climate change in the region.

[0121] In this implementation, the specific steps of establishing a single-site historical database and a prediction factor database and performing feature selection include:

[0122] S11: Get data.

[0123] The data in this example comes from the global long-term flux observation network FLUXNET (https: / / fluxnet.org / ), and the temporal resolution of the sample data is daily. The single-site historical database and forecast factor database include but are not limited to the following related variables:

[0124] Time series: Day of sequence (DAY), Day of year (DOY).

[0125] Meteorological forcing: Precipitation (P0), Temperature of air (TA), Wind direction (WD), Wind speed (WS), Atmospheric pressure (PA), Vapor Pressure Deficit (VFD), Relative humidity (RH), Friction velocity (USTAR), Soil temperature (TS);

[0126] Radiative flux: Shortwave radiation incoming (SW_IN), Shortwave radiation outgoing (SW_OUT), Longwave radiation outgoing (LW_OUT), Net radiation (NETRAD);

[0127] Energy flux: soil heat flux (G), latent heat flux (LE), sensible heat flux (H), net ecosystem exchange efficiency (NEE), ecosystem respiration efficiency (RECO), gross primary production (GPP), photosynthetic photon flux density (PPFD);

[0128] Lag: Precipitation lagged 1-14 days (P1-P14), soil moisture lagged 1-14 days (SM1-14).

[0129] S12: Feature selection.

[0130] Since FLUXNET provides multiple variables related to the same physical quantity at a single site, including physical quantities from different data sources, physical quantities generated by different data processing and filling methods, original observations of physical quantities, and different percentiles, standard deviations, uncertainties, etc. of physical quantities, multiple related variables of the same physical quantity not only increase the computational intensity for machine learning models, but also lead to poor model training due to data uncertainty. Therefore, the following two feature selection schemes are made:

[0131] Empirical selection: First, based on experience, remove features that are obviously linearly correlated in the feature set, such as percentiles, etc.; second, delete data indicator features that are obviously irrelevant to the forecast, such as uncertainty, standard deviation, etc. For example, FLUXNET provides 5 variables related to the average daily temperature, including TA_F_MDS (average daily temperature of the method combining MDS+ERA-interim), TA_F_MDS_QC (data quality index of the average daily temperature of the method combining MDS+ERA-interim and the MDS calibration method), TA_ERA (ERA-interim average daily temperature), TA_F, TA_F_QC (data quality index of the average daily temperature of ERA-interim). In this example, based on experience, only TA_F_MDS is retained as the feature of the average daily temperature. In this embodiment, the original input variables are 320, and 115 variables are retained after empirical selection.

[0132] Sequential backward selection (SBS): From the full set of features X d ={x1,x2,…,x d}, each time feature x is removed from feature set X, so that the evaluation function value J(x) reaches the maximum after feature x is removed. Starting from the entire feature set, each time feature x that minimizes the evaluation function J(x) is selected and removed. The specific steps are as follows:

[0133] (1) Initialize k=d, where d is the number of features of the original input, which is 115 in this embodiment.

[0134] (2) The evaluation function is determined as the difference in mean squared error (MSE) before and after deleting the feature, specifically:

[0135]

[0136] (3) From X d Remove x from the equation, and x satisfies J(x) to reach the maximum, specifically:

[0137] x=argmax J(x)

[0138] (4) Update the feature dataset, i.e. X k-1 =X k -x. Repeat (2)-(3) until the MSE of the retained feature data set meets the prediction requirements.

[0139] like Figure 3 As shown in Figure 2, 115 features are retained after empirical selection and input into SBS for feature selection. 2To achieve stability, the top 16 best features are selected as the input variables of the final prediction model in order to improve the accuracy of model training and automatically eliminate redundant variables that reduce accuracy from the perspective of contribution to model accuracy.

[0140] In this embodiment, by applying the above two feature selection methods, redundant features are eliminated from 320 original features, and finally 16 features that contribute the most to soil moisture forecast are retained, as shown in Table 1.

[0141] Table 1 16 input features of the random forest forecast model and their descriptions

[0142]

[0143]

[0144] S2: Establish a black box model for hydrological and meteorological factor forecasting and calculate the forecast error (R 2 , MAE and RMSE): S21: Establish a black box model for hydrological and meteorological factors forecasting.

[0145] The optional black box models include random forest, naive Bayesian regression, support vector machine, deep learning, etc. The first 70% of the time series data in the data set is used as the training set, and the last 30% is used as the test set to calculate the forecast error;

[0146] The machine learning model established in this embodiment is a random forest model, with the first 70% of the time series as a training set and the last 30% as a test set. Random forest is implemented as follows: first, subsets are divided from the training sample, and then the subsets are trained using a CART regression tree. The CART regression tree uses the minimum mean square error method to divide the feature data set. The principle of the CART regression tree is shown in formula (1). Generate multiple CART regression trees, and the final prediction result of each CART regression tree is the mean of the leaf nodes reached by the sample point. The final prediction result of the random forest is the mean of the prediction results of all CART regression trees, as shown in the following formula (1):

[0147]

[0148] Among them, A is the partition feature, s is the partition point of feature A, D1 and D2 are the two data sets divided by s, c1 and c2 are the mean of the sample outputs of the D1 and D2 data sets, respectively, and x i are the input features in datasets D1 and D2, y i is the output of the CART regression tree.

[0149] S22: Statistical forecast error

[0150] The statistical forecast error is R 2 (R-square,R2 ) Mean Absolute Error (MAE), Mean Squared Error (MSE), and Root Mean Squared Error (RMSE), see equations (2)-(5).

[0151]

[0152]

[0153]

[0154]

[0155] In the formula, n represents the number of soil moisture data to be predicted, represents the observed value of soil moisture at the i-th level, y i Represents the forecast value of the i-th hydrological and meteorological element.

[0156] The black box model provided in this embodiment is a random forest model. Generally, random forest is considered to be a more practical method for describing the complex interactions between soil moisture and different land surface variables. Compared with gradient descent tree and naive Bayesian regression, random forest has better performance in predicting soil moisture from FLUXNET sites in the test.

[0157] The complete data set used in this example is 1096 daily observation data from 2003 to 2005. Here, the first 70% of the time series data set (from January 1, 2003 to February 5, 2005) is used as the training set. At the same time, the last 30% of the data (from February 6, 2005 to December 31, 2005) is used as the test set. Considering the calculation and accuracy, the number of decision trees in the random forest is set to 100, and the number of nodes in each tree is set to 5. Figure 4 (a) and Figure 4 As shown in (b), the forecast results are obtained through the random forest model forecast, and the distribution of forecast values ​​and observed values ​​in the test set is obtained, as well as the time series comparison of forecast values ​​and observed values ​​in the test set.

[0158] In step S2, the statistical prediction error provided by this embodiment is R 2 , MAE, MSE and RMSE. In this example, the R 2 , MAE and RMSE are 0.975, 1.078, 0.190 and 0.036 respectively.

[0159] S3: The post hoc interpretation method explains the machine learning prediction model from multiple perspectives and explains the black box model through two approaches: one is to quantify the importance of features, including feature importance based on MSE and Shapley method; the other is to evaluate the effect of features, that is, to predict the output change when a single or multiple features change, thereby reflecting the impact of features on output changes. Methods include partial dependence diagrams, individual conditional expectations, and cumulative local effects;

[0160] In this embodiment, the specific steps of the post-explanation method for explaining the machine learning prediction model from multiple perspectives include:

[0161] S31: MSE-based feature importance (MFI)

[0162] This method is based on a tree model (such as random forest) and calculates importance by arranging out-of-bag data. Random forest uses a random and replacement extraction mechanism when constructing each tree. Therefore, for each tree, some training instances do not participate in the generation of the tree. These instances are called out-of-bag data samples of the tree. First, random forest forms several decision trees based on the training data. Taking the kth tree as an example, the out-of-bag data of the kth tree is used to test it to obtain the accuracy of the prediction model (this example is based on MSE). Then, the feature values ​​of a single feature of the out-of-bag data are randomly arranged to obtain a new MSE, and the difference between the MSE before and after the arrangement is calculated. Finally, the MSE difference of each tree is averaged, and this is used as a standard to rank the importance of the features based on the accuracy of the prediction model.

[0163] The implementation steps are as follows:

[0164] S311: Assume that the feature importance of the k-th variable is evaluated. For the trained random forest, there are m trees in total. For the j-th tree, the out-of-bag feature dataset D = {x1, x2, ..., x k ,…,x m} as input for prediction and get the prediction result Y j,k , see formula (6):

[0165] Y j,k =f j (D) (6)

[0166] S312: Compare with the soil moisture observation value and calculate the out-of-bag feature data set D = {x1, x2, ..., x k ,…,x m} prediction error MSE j,k , see formula (7):

[0167]

[0168] In the formula, is the observed soil moisture value.

[0169] S313: Randomly perturb the kth variable to form a new out-of-bag feature data set D′={x1,x2,…,x k ′,…,x m}, then use the new out-of-bag feature data set for verification, repeat steps (S311-S312), and calculate the prediction error to get MSE j,k ′.

[0170] S314: For the kth variable, calculate its importance (MFI) is the MSE after perturbation j,k ′ and the unperturbed MSE j,k The average value of the difference is shown in formula (8).

[0171]

[0172] In this embodiment, the number of random forest trees is set to 100, and the training data is input into the model to form 100 decision trees. The training instances of each tree are randomly extracted with replacement. Taking the MFI of SM1 as an example, for the kth tree, the data that does not participate in the training of the kth tree is the out-of-bag data, and the out-of-bag data is D = {SM1, SM2, ..., LE} as the input of the kth tree, and the MSE is calculated. The out-of-bag data SM1 is randomly arranged as SM1' to form new out-of-bag data D' = {SM1', SM2, ..., LE}, which is input into the kth tree to obtain the MSE' of the new out-of-bag data. Calculate the ΔMSE before and after the arrangement, and take the average of the ΔMSE of the 100 trees as the MFI value of the SM1 feature. Repeat the above process for each feature to obtain the MFI value of each feature and sort them.

[0173] like Figure 5 As shown in the figure, the results of feature importance based on MSE are obtained. The factors that have the greatest impact on the soil moisture forecast accuracy of the Ha2 station are P1, SM1, SM2, P2, SM3 and DAY. This is basically consistent with the research on this site, that is, recent lagged precipitation and soil moisture are the main control factors of soil moisture. In this area, soil moisture does not change much due to the slow evaporation or infiltration of water. Therefore, lagged soil moisture has an important impact on soil moisture. And lagged precipitation also has a similar direct impact on soil moisture. This is also in line with the characteristic that soil moisture in alpine shrub areas is more sensitive to single rainfall events. These corresponding parts, combined with the prior knowledge of hydrological and meteorological processes, verify that the trained random forest model has the ability to capture important features.

[0174] MFI provides a visualization tool for ranking feature importance, which reflects the relationship between the input and output of the black box model from the perspective of model accuracy, and helps to improve the credibility of the results of the hydrological forecast black box model. On the contrary, if the important features calculated by MFI are significantly different from the main influencing factors of soil moisture in the region, even if the forecast accuracy of the test set is high, the model still lacks credibility and its generalization ability is greatly lacking, making it difficult to use in actual hydrological forecasting.

[0175] S32: Shapley method (Shapley values, SV)

[0176] Shapley is an explanatory model inspired by game theory, where all features are considered contributors. The SV value of each feature is the weighted sum of the contribution to the predicted value over all possible feature combinations. The SV method can be applied to any subset. A training set of time series P = {x1, x2, ..., x k ,…,x p}, p is the number of features, x k is the feature of interest, calculate x k The SV of the feature is formula (9):

[0177]

[0178] In the formula, is the SV value of the kth feature of interest, S is the set of all feature combinations, S i is the i-th feature combination in the S set, is a trained forecast model.

[0179] Steps to implement the Shapley method:

[0180] S321: Obtain feature set P = {x1, x2, ..., x k ,…,x p}, extract the feature x of interest k , remove x k The feature set is T = {x1, x2, ..., x p};

[0181] S322: Feature set T forms different combinations S i (such as S1 = {x1}), all combinations S i They form a set S;

[0182] S323: Calculate feature x through formula (9) k of That is the SV value of the kth feature.

[0183] The SV method quantifies the contribution of each feature to the forecast value as a proportion of the forecast value, thereby obtaining the importance of each feature. In this embodiment, two types of SV local interpretation distributions are provided, including time series and eigenvalue distributions. The time series SV interpretation intuitively reflects the change in the contribution of each feature over time. If a feature with a positive effect (such as precipitation) has a negative effect in a certain instance in the time series SV interpretation, even if the forecast is accurate, it can be considered that the feature of the instance has an abnormal value, resulting in the model's forecast in this instance being unreliable. The eigenvalue distribution SV interpretation reflects the impact of the high and low eigenvalues ​​on the change in the size of the feature contribution. Based on the special land-atmosphere exchange characteristics of the site, if a feature appears to be contrary to the hydrological physical prior knowledge or observational experience, then in the forecast model, it may be insensitive to changes in certain features, which also causes the forecast value to be unreliable, and the model is improved from within, thereby improving the accuracy and physical consistency of the forecast.

[0184] S33: Partial dependence plot (PDP)

[0185] When the partial dependence diagram focuses on the change of a single feature, the average value of the forecast value change in the entire data set is used to characterize the average change of soil moisture in the forecast time series when a local factor changes. It can reflect the average change of local soil moisture when it is affected by a certain factor.

[0186] The partial dependence plot shows how features affect the predicted value by showing the marginal effect of the feature on the predicted value, making the model internally transparent. The marginal effect refers to the change in effect obtained by changing a certain feature when other features remain unchanged.

[0187] In the partial dependence plot, the marginal effects are expressed as partial dependence values. The form of expression:

[0188] A training set of time series P = {x1, x2, ..., x k ,…,x p}, p is the number of features, x k is the set of k-th features, x k P is divided into two subsets, namely S = {x k} and C={x1,x2,…,x k-1 ,x k+1 …,x p Generate x k The set of possible values ​​x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m}, m is x kThe number of possible values.

[0189] Use x s The possible values ​​replace the feature x one by one k , calculate the average of the forecast values ​​obtained for each substitution, which is the partial dependence value See formula (10):

[0190]

[0191] In the formula, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of Except x k All features except (set C), n is the number of instances in the time series, i is the i-th data in the time series, is the forecast model. Get the continuous set x s Corresponding Value, in x s The element is the horizontal axis, The value is the vertical axis, and the partial dependence diagram is obtained, which visualizes the impact of feature changes on the predicted value. The implementation steps of the partial dependence diagram are:

[0192] S331: Assume that each feature is independent of each other, and the input data set is P = {x1, x2, ..., x k ,…,x p}, feature x of interest k P is divided into two subsets, namely S = {x k} and C={x1,x2,…,x k-1 ,x k+1 …,x p};

[0193] S332: x k The continuous set of possible values ​​is x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m};

[0194] S333: Calculate the set x of possible values s,j The partial dependence value of The feature x k All by x s,j replace, Fixed, obtain n forecast values ​​and average them, and then calculate the average value of the forecast values ​​obtained each time, which is the partial dependence value As formula (10);

[0195] S334: The calculated value in S333 The value is continuous x s Visual presentation can show x k How the predicted value changes when the characteristics change;

[0196] S335: Replace the feature and repeat steps S332-S334 to obtain the change of the predicted value under the change of the feature, thereby obtaining whether the influence of the feature on the predicted value is linear or nonlinear.

[0197] For example, PDP obtains the possible values ​​of DOY as [1-365], and substitutes the possible values ​​(such as DOY = 1) into all DOY in the training set, and other features remain unchanged to obtain new forecast values. Then the PDP value when DOY = 1 is the average value of the new forecast values. Perform the above operations for DOY = 1, 2, ..., 365 respectively, and obtain the PDP graph of DOY features in the black box model. This graph reflects the changing trend of the average forecast assuming that DOY changes, eliminating the influence of other factors, and highlighting the influence of DOY (time period change) on the average forecast.

[0198] like Figure 6 As shown, the SV explanation of the time series is obtained, and the baseline in the figure is 0, indicating the average predicted soil moisture. A positive SV of a feature means that the feature has a positive impact on the soil moisture forecast at this time, and vice versa. In this embodiment, the lagged soil moisture and soil temperature have a major impact on the predicted soil moisture (more than 90% of the value). It is revealed that the time memory of soil moisture plays an important role in soil moisture forecasting. This phenomenon is consistent with the characteristics of higher soil moisture retention in alpine and shrub areas. Therefore, the SV of the time series can reveal the importance of dynamic characteristics. At the same time, the results are similar to the characteristics of alpine and shrub areas, which further improves the reliability of the random forest forecast model. In addition, PISW, P1, TA and DOY also contribute to SM prediction.

[0199] like Figure 7 (a)- Figure 7 As shown in (c), we get the SV interpretation of the eigenvalue distribution. Figure 7 It can be clearly seen from (a) that lagged soil moisture is the main driving factor of soil moisture. It can be seen that the longer the lag time, the greater the impact on soil moisture, while the increase in the number of lag days will weaken this effect.

[0200] exist Figure 7In (b), its distribution pattern can be seen more clearly. The precipitation values ​​of almost all instances are very low. P1 of one-day lagged precipitation is the highest, followed by P0 of the same-day precipitation. This is consistent with the phenomenon that the developing plant canopy intercepts the precipitation of the day, resulting in the contribution of the precipitation of the day P0 being lower than the precipitation of one-day lagged precipitation P1. However, low-value precipitation has a negative SV phenomenon. This is contrary to the positive effect of precipitation on soil moisture. The eigenvalue distribution SV interpretation discovered this, revealing that the random forest model is not sensitive to low-value precipitation data, leading to incorrect judgments.

[0201] like Figure 7 Other characteristic effects of SV are obtained as shown in (c). The trends of soil temperature (TS) and air temperature (TA) are similar, indicating that soil moisture generally increases with increasing temperature. This can be explained by the melting of snow in high mountainous areas. When TS increases, the melting of the shallow active soil layer causes a rapid increase in soil moisture. DOY shows mixed changes in SV and eigenvalues, which is consistent with the annual cycle of soil moisture. For the DAY variable, soil moisture has an overall increasing trend, especially in 2005, which is consistent with the Figure 4 (b) consistent.

[0202] The above SV interpretation results are basically consistent with hydrological and meteorological knowledge, which improves the trust in the prediction results of the black box model. However, it also reveals some defects of the black box model itself, such as insensitivity to less low-value precipitation data, providing a basis for further improving the black box model.

[0203] S34: Individual conditional expectation (ICE)

[0204] Individual conditional expectation is a method to characterize the process of the influence of features on the forecast value through marginal effects, reflecting the influence of features on the results in a single instance forecast. In the forecast interpretability of soil moisture, individual conditional expectation has two functions: 1) how a single feature affects the forecast value in extreme cases (such as special disasters, etc.), judging the importance of features, and further reflecting the feature effect; 2) through statistical analysis of time series, obtaining the changes in the effects of each feature on the forecast value in different seasons / years.

[0205] The partial dependence of the feature of each instance on the predicted value in the individual conditional expectation

[0206] The data set P at a single time point = {x1, x2, ..., x k ,…,x p}, p is the number of features, x k is the feature of interest, x k P is divided into two subsets, namely S = {xk} and C={x1,x2,…,x k-1 ,x k+1 …,x p Generate x k The set of possible values ​​x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m}, m is x k The number of possible values.

[0207] Use x s The possible values ​​replace the feature x one by one k , calculate the forecast value obtained for each substitution, which is the partial dependence value at a single time point

[0208]

[0209] In the formula, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of x c Except x k All features except (set C), is the forecast model. Get the continuous set x s Corresponding Value, in x s The element is the horizontal axis, The value is the vertical axis.

[0210] Taking rainfall characteristics as an example, in the historical data of this site, the daily rainfall fluctuates between 0 and 30 mm. Assuming that the number of continuous values ​​of ICE is 10, then the possible value of rainfall is x. s =[0,3,6,9…30]. In the single-day soil moisture forecast, x s Substitute them into the trained machine learning model one by one, keep other variables unchanged, and predict new results And x s and The ICE diagram is formed accordingly. In the process of predicting soil moisture for a single day, the ICE method changes the value of a single feature, thereby determining the change of the forecast value under the change of the feature, reflecting the impact of feature changes on soil moisture forecast.

[0211] Individual expectations:

[0212] like Fig. 9 (a)- Fig. 9As shown in (p), the ICE graph of the selected 16 features is obtained. Each curve in the figure represents the ICE of an instance. Due to the use of instance-based observation methods, the observation results show heterogeneity, making the changes in the observation results easy to understand. Fig. 9 Taking (d) as an example, the ICE curve is clearly divided into two branches according to different trends. The first branch is located in the upper part of the entire ICE graph. The other part is below. The results show that the well-trained random forest model is able to capture the relationship between PISW and soil moisture. This relationship is mainly reflected in two aspects and is associated with the two branches of the ICE graph. The first branch represents the "negative feedback" of PISW on soil moisture. Among them, as PISW increases, the surface radiates more energy, which suppresses evapotranspiration and inhibits the reduction of soil moisture. Therefore, it can be observed that the ICE maintains the level of soil moisture. The other branch of ICE can be explained by the "positive feedback" of PISW on soil moisture. When soil moisture is low, hydrothermal effect is activated. These explanations can reveal the changing rules discovered by the black box model and help users trust the model.

[0213] The ICE plot also shows that the random forest model can capture some special precipitation phenomena (such as Fig. 9 of (n), Fig. 9 (o) and Fig. 9 When soil moisture is less than 12%, the precipitation characteristics (P0, P1, and P2) have little effect on the predicted soil moisture. This is because, in this embodiment, the site is located on a frozen plateau, and there are very few instances of heavy rainfall, which leads to inherent bias in the data itself. Therefore, the uneven distribution of data destroys the model, making it impossible to grasp the true relationship.

[0214] In this embodiment, several examples are provided to illustrate the functionality of the ICE graph. PISW's ICE vividly demonstrates the heterogeneity of instances, while precipitation's ICE reveals the limitations of the black box model in real space. In summary, ICE is a useful visualization tool that can display feature effects from a per-instance view.

[0215] S35: Accumulated local effect (ALE)

[0216] The cumulative local effect obtains the forecast difference through the conditional distribution of the feature, and accumulates the forecast difference to reflect the feature effect. First, the features of the time series and analysis are extracted, and their possible values ​​are continuous and gridded to form a two-dimensional grid; secondly, the vertex values ​​of the grid are replaced by the values ​​of the features of the data points in the grid (a single case is a point), and the other features are fixed and forecasted; finally, the four forecast values ​​are subtracted to form a second-order difference, and each second-order difference is accumulated and centered to obtain the cumulative local effect, and then drawn into a cumulative local effect diagram.

[0217] This embodiment uses one-dimensional ALE and two-dimensional ALE methods to reflect the characteristic effect. The one-dimensional ALE and two-dimensional ALE methods are similar, and the specific steps of two-dimensional ALE are as follows:

[0218] S351: extracting the feature group of interest, making its possible values ​​continuous and gridded to form a two-dimensional grid;

[0219] S352: Replace the feature values ​​of the data points (a single point in a single case) in the grid with the vertex values ​​of the grid respectively, and keep other features unchanged, and make a prediction;

[0220] S353: Subtract the four soil moisture forecast values ​​to form a second-order difference, accumulate and center each second-order difference to obtain the cumulative local effect And plotted as ALE. The second-order difference of a single data point in a single grid is calculated as:

[0221]

[0222] Where f is the soil moisture forecast model, j is the first feature of interest, l is the second feature of interest, K is the total number of grids, k is the kth interval of the first feature of interest, m is the mth interval of the second feature of interest, i is the i-th data point in the grid, and z is k,j is the eigenvalue corresponding to the k interval in the first feature of interest, z m,l is the value of the feature corresponding to the mth interval divided by the second feature l, x i,\{j,l} is all the features of the i-th data point except the first feature j and the second feature l. The second-order difference accumulation of all points in all grids is:

[0223]

[0224] In the formula, is the cumulative second-order difference of all points in the uncentered data set, k j (x) is the total number of intervals into which the first feature j of interest is divided, k l (l) is the total number of intervals into which the second feature l of interest is divided, n{j,l} (k,m) is the total number of data points in the grid formed by the second feature l interval of interest in the kth time interval, N {j,l} (k,m) is the area of ​​the interval. All points in each grid are centered:

[0225]

[0226] In the formula, is the accumulated second-order difference after centralization. Centralization of all grids:

[0227]

[0228] In the formula, It is the globally centralized second-order difference accumulation.

[0229] Taking the interaction of the two features PISW and TS as an example, assuming that the value range of PISW is [250~500], the value range of TS is [-10~15], and the number of grid divisions K is 10, then the possible values ​​of PISW are [250,275,……500], and the possible values ​​of TS are [-10,-7.5……15]. Then, the possible values ​​of the two features form a grid. Substitute the four vertex values ​​of each grid node into the trained model, and keep other features unchanged to obtain 4 new soil forecast values. By making differences, accumulation and decentralization, the changes in soil moisture forecast under the interactive changes of PISW and TS are obtained. Therefore, the changes in soil moisture forecast can be visualized under the possible interactive changes of PISW and TS. Judging the performance of machine learning models under the possible interactive conditions based on whether they conform to the actual physical process is a method to visualize the input-output relationship of black box models.

[0230] Partial dependence plot in one dimension and cumulative local effects in one dimension:

[0231] Both the one-dimensional PDP and one-dimensional ALE show the changes in the output results when a single feature changes. The premise for the one-dimensional PDP is that there is no interaction between the features, which is obviously not true in this example. Therefore, this example is explained with a one-dimensional ALE diagram. Figure 8 (a)- Figure 8 As shown in (p), the 1D ALE graph of the RF model for all 16 features is obtained. The model established by the 1D ALE has the ability to capture the changes in soil moisture with the changes in features. A typical example is the soil temperature of the 1D ALE (see Figure 8 (e)). Figure 8(e) shows how soil temperature changes soil moisture, which can be divided into three sub-periods. The melting temperature starts from -10℃ to -2℃. As the soil temperature increases, the soil moisture content shows a slight upward trend. When the soil temperature is close to 0℃, the melting rate is very high, and most of the melted snow leaches into the soil, resulting in a significant increase in soil moisture. After that, as the soil temperature rises, the proportion of evapotranspiration increases, the higher the soil temperature, and the lower the soil moisture. The random forest model can reflect this physical phenomenon through data, which increases the credibility of the model. One-dimensional ALE can reveal the dynamic relationship between features and forecast results.

[0232] 2D ALE:

[0233] like Figure 10(a)-Figure 10(b) As shown in Figure 10, a two-dimensional ALE example is used for explanation. Figure 10(a) reflects the predicted change of SM under the interaction of TS and PISW. When TS is <0℃ and PISW is less than 300W / m 2 When PISW increases, soil moisture becomes relatively low. This is consistent with the natural phenomenon that during the frozen period, the shrub forest is covered by heavy snow, and there is less liquid water in the soil, leading to drought. As PISW increases, the random forest model predicts drought.

[0234] Another typical example is shown in Figure 10(b). The distribution of example points can be seen. When TS is above 0℃, soil moisture decreases as LE and TS increase. This can be explained by the joint effect of LE and TS, both of which have a negative impact on soil moisture content. This combination enhances this effect on soil moisture. Among them, there is an obvious blue part in Figure 10(b), indicating that at 0℃, as LE increases, the soil becomes moist. In fact, this rarely happens in nature. Therefore, two-dimensional ALE can not only intuitively show the combined effect of the two features on the predicted SM, but also reveal some unrealistic phenomena captured by the model.

[0235] 2D ALE is a viable visualization tool to explore the interactive effects of coupling features on the predicted SM, but it changes significantly with the number of bins. Its practical application depends on the coupling features of interest.

[0236] S4: Based on the post-explanation method, an integrated system is developed in the python environment to interpret the machine learning model for hydro-meteorological forecasting: First, the model results are evaluated for credibility using the feature importance ranking method, including the MFI for global interpretation and the Shapley method for local interpretation. If the top 5 features in the MFI are proven to be important factors for the site, have physical consistency, and most of the explanations in Shapley's local interpretation are consistent with the actual hydro-physical process of the site, the model can be considered to have a certain degree of credibility. Otherwise, the machine learning model needs to be reselected. Secondly, the interpretability analysis of the feature effect is performed to diagnose the model, including the partial dependence diagram and cumulative local effect of the global interpretation, and the individual conditional expectation of the local interpretation. If the feature effect shown by the global partial dependence diagram or the cumulative local effect is inconsistent with the actual hydro-physical process, it is considered that the machine learning model has defects in the identification of the feature, and the model needs to be repaired on the feature, including adding weights, modifying the decision tree structure, etc., so that the feature effect has hydro-physical consistency, and the model is retrained. For the forecast of single-day hydrological and meteorological elements, the individual expectation conditions of local interpretation can be used to find the abnormal features in the single-day forecast, and determine whether they are abnormal values ​​of the features based on historical data. If they are not abnormal values, the internal structure of the machine learning model is modified by increasing the instance weight and other methods (determined according to the choice of the machine learning model). Finally, the above methods are integrated into a complete set of model-independent post-explanation tools, including four modules: data reading, model training, accuracy evaluation and interpretability tools.

[0237] This embodiment uses soil moisture forecast as an example of hydrological and meteorological forecasting to illustrate that the proposed integrated post-interpretation tool can overcome the application areas of black box models from the following aspects:

[0238] (1) Improve the evaluation system of machine learning models. Usually, only the prediction accuracy of machine learning models is evaluated. This tool evaluates the interpretability of machine learning models, based on local hydrological physics knowledge, to assist hydrological forecasters in further selecting and interpreting models.

[0239] (2) Evaluate the credibility of model results based on feature importance. Taking the local hydrological dominant process as the standard is conducive to reflecting the credibility of the model and the reliability of the results.

[0240] (3) Characteristic effects provide a visual basis for diagnosing models. Characteristic effects can reflect defects in the model or data itself, such as inherent bias in the data and the model’s recognition bias of characteristic data. The visualization of characteristic effects can reflect the overall average effect of the model from a global perspective (PDP, ALE), and can also reflect the phenomenon of the characteristics in a single hydrological forecast instance that are inconsistent with hydrological physics due to data or model bias from a local perspective (ICE). This provides a basis for further debugging and repairing the model.

[0241] In summary, with the application of explainable AI tools in hydrometeorological forecasting, the application of AI technology will be promoted, and this tool can contribute to the establishment of a smarter and more accurate hydrometeorological forecasting system.

[0242] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

Claims

1. A visual post-interpretation method for an integrated hydro-meteorological machine learning forecast model, characterized in that: The steps include: Establish a single-site historical database and a predictor database, and perform feature selection; Establish a black box model for hydrological and meteorological factor forecasting and calculate forecast errors; Post hoc interpretation methods explain machine learning forecasting models from multiple perspectives, including feature importance based on MSE, Shapley's method, partial dependence plots, individual conditional expectations, and cumulative local effects. The feature importance based on MSE is explained as follows: multiple decision trees are formed based on the training data, out-of-bag data are selected to test the decision trees, the prediction error MSE of the prediction model is obtained, the feature values ​​of a single feature of the out-of-bag data are randomly arranged, a new prediction error MSE is obtained, and the difference in the prediction error MSE before and after the arrangement is calculated, the prediction error MSE difference of each decision tree is averaged, and the features are ranked in importance from the prediction model accuracy according to the average value; The explanation is based on Shapley's method, which is as follows: the contribution of each feature to the predicted value is converted into the proportion of the predicted value, and the contribution to the predicted value is weighted summed over all possible feature combinations to obtain the importance of each feature; The explanation is based on the partial dependence diagram, specifically: select the focus feature in the time series training data, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the average value of the forecast value obtained by each replacement, which is the partial dependence value, which is used to reflect the average change of the local hydrological and meteorological forecast object when it is affected by a certain factor; The explanation is based on individual conditional expectation, specifically: select the focus feature in the training data at a single time point, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the forecast value obtained by each replacement, which is the partial dependence value at a single time point, used to reflect the impact of the feature on the result in the forecast of a single instance; The explanation is based on the cumulative local effect. Specifically, the attention feature in the training data at a single time point is selected, and the training data is divided into two subsets based on the attention feature to generate a distribution grid of the attention feature. In a single grid, the data value of the attention feature distribution grid is used to replace the attention feature one by one, and the forecast value obtained by each replacement is calculated. The differences in the forecast values ​​in the grid are averaged and then accumulated on the grid to reflect the feature effect. The reliability of the model results is evaluated based on the feature importance of MSE and the Shapley method. The interpretability analysis of the feature effects is performed based on the partial dependence diagram, individual conditional expectation and cumulative local effect. This is used to diagnose the model. For the forecast of single-day hydrological and meteorological elements, the abnormal features in the single-day forecast are found through individual expectation conditions, and whether they are abnormal values ​​of the features is determined based on historical data.

2. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: The single-site historical database and forecast factor database are established, and the data in the database include time series, meteorological forcing, radiation flux, energy flux and hysteresis; Feature selection methods include empirical selection and sequential backward selection.

3. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: The hydrological and meteorological factor forecasting black box model adopts any one of a random forest model, a naive Bayesian regression model, a support vector machine model or a deep learning model.

4. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: The explanation is based on the Shapley method, and the specific steps include: Get the training set, extract the features of interest, and calculate the SV value of the features of interest, expressed as: Among them, p represents the number of features, x k Characteristics that indicate concern, represents the SV value of the kth feature of interest, S represents the set of all feature combinations, S i represents the i-th feature combination in the S set, Represents the trained forecast model.

5. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: The explanation based on individual conditional expectations includes the following specific steps: Select the feature of interest x in the training data k , based on the focus feature x k The training data is divided into two subsets, namely S = {x k } and C={x1,x2,…,x k-1 ,x k+1 …,x p }; We will focus on feature x k The continuous set of possible values ​​is x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m }; Compute the set of possible values ​​x s,j The partial dependence value of We will focus on feature x k By x s,j replace, Fixed, obtain n forecast values ​​and average them, and then calculate the average value of the forecast values ​​obtained each time, which is the partial dependence value The specific calculation formula is expressed as: Among them, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of Except x k All features except , n is the number of instances in the time series, i is the i-th data in the time series, It is a forecast model; The partially dependent values ​​will be calculated By continuous x s Perform a visual presentation to show the focus feature x k Changes in the forecast value when the change occurs; The feature is replaced iteratively and repeatedly to obtain the change of the predicted value under the change of the feature, so as to obtain whether the influence of the feature on the predicted value is linear or nonlinear.

6. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1, characterized in that: The explanation based on individual conditional expectations includes the following specific steps: Select a data set P at a single time point = {x1, x2, ..., x k ,…,x p }, p is the number of features, x k is the feature of interest, x k P is divided into two subsets, namely S = {x k } and C={x1,x2,…,x k-1 ,x k+1 …,x p }; Generate x k The set of possible values ​​x s ={x s,1 ,x s,2 ,…,x s,j ,…,x s,m }, m is x k The number of possible values; Use x s The possible values ​​replace the feature x one by one k , calculate the forecast value obtained for each substitution, which is the partial dependence value at a single time point Among them, x s is feature x k The set of possible values ​​of , j is the continuous set x s The jth value of x c Except x k All features except is the forecast model, and obtains the continuous set x s Corresponding Value, in x s The element is the horizontal axis, The value is the vertical axis.

7. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: The explanation based on cumulative local effects includes the following specific steps: Extract the two relevant features of interest, make their possible values ​​continuous and gridded to form a two-dimensional grid; The vertex values ​​of the grid are used to replace the feature values ​​of the data points in the grid, while other features remain unchanged, and are input into the black box of the machine learning model for prediction; The four forecast values ​​obtained are subtracted to form second-order differences, each second-order difference is accumulated and centered to obtain the cumulative local effect, which is then plotted into a cumulative local effect diagram.

8. The visual post-interpretation method of the integrated hydrological and meteorological machine learning forecast model according to claim 1 is characterized in that: One-dimensional ALE or two-dimensional ALE is used to reflect the impact of feature changes on forecast results.

9. A visual post-explanation system integrating hydrological and meteorological machine learning forecast models, characterized by: include: Database construction module, feature selection module, model construction module, forecast error statistics module and interpretation module; The database construction module is used to establish a single-site historical database and a prediction factor database; The feature selection module is used to perform feature selection; The model building module is used to establish a black box model for hydrological and meteorological factor forecasting; The forecast error statistics module is used to count the forecast errors; The explanation module is used to explain the machine learning prediction model from multiple perspectives using post-hoc explanation methods, including feature importance based on MSE, Shapley method, partial dependence plot, individual conditional expectation and cumulative local effect; The feature importance based on MSE is explained as follows: multiple decision trees are formed based on the training data, out-of-bag data are selected to test the decision trees, the prediction error MSE of the prediction model is obtained, the feature values ​​of a single feature of the out-of-bag data are randomly arranged, a new prediction error MSE is obtained, and the difference in the prediction error MSE before and after the arrangement is calculated, the prediction error MSE difference of each decision tree is averaged, and the features are ranked in importance from the prediction model accuracy according to the average value; The explanation is based on Shapley's method, which is as follows: the contribution of each feature to the predicted value is converted into the proportion of the predicted value, and the contribution to the predicted value is weighted summed over all possible feature combinations to obtain the importance of each feature; The explanation is based on the partial dependence diagram, specifically: select the focus feature in the time series training data, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the average value of the forecast value obtained by each replacement, which is the partial dependence value, which is used to reflect the average change of the local hydrological and meteorological forecast object when it is affected by a certain factor; The explanation is based on individual conditional expectation, specifically: select the focus feature in the training data at a single time point, divide the training data into two subsets based on the focus feature, generate a set of possible values ​​of the focus feature, use the data values ​​in the set of possible values ​​of the focus feature to replace the focus feature one by one, and calculate the forecast value obtained by each replacement, which is the partial dependence value at a single time point, used to reflect the impact of the feature on the result in the forecast of a single instance; The explanation is based on the cumulative local effect. Specifically, the training data is divided into two subsets based on the feature of interest, and a distribution grid of the feature of interest is generated. In a single grid, the data value of the distribution grid of the feature of interest is used to replace the feature of interest one by one. The forecast value obtained by each replacement is calculated, and the difference of the forecast values ​​in the grid is averaged and then accumulated on the grid, thereby reflecting the impact of the feature change on the forecast result. The reliability of the model results is evaluated based on the feature importance of MSE and the Shapley method. The interpretability analysis of the feature effects is performed based on the partial dependence diagram, individual conditional expectation and cumulative local effect. This is used to diagnose the model. For the forecast of single-day hydrological and meteorological elements, the abnormal features in the single-day forecast are found through individual expectation conditions, and whether they are abnormal values ​​of the features is determined based on historical data.

Images