A runoff prediction method based on Bayesian average combined with vine copula

Abstract

The present invention discloses a runoff prediction method based on Bayesian average combined with vine copula, and constructs a new model for runoff prediction based on Bayesian model average (BMA) combined with vine copula (ViC) function, which can effectively combine multiple prediction variables in high-dimensional situations, and can give a clear expression, which can significantly improve the prediction accuracy of monthly runoff and realize reliable prediction of monthly runoff under different forecast periods. The invention has high applicability and extensibility, broadens the theoretical framework of runoff prediction, improves the current runoff prediction accuracy, and can provide technical support for the comprehensive construction of water resources management, risk assessment and early warning systems in the basin.

Description

Technical Field

[0001] The invention relates to the technical field of hydrological forecasting, and in particular to a runoff prediction method based on Bayesian averaging combined with vine copula. Background Art

[0002] Changes in river runoff are directly related to water resource security, flood control and drought relief, agricultural production activities, and ecological and environmental health. In recent years, affected by global warming and human activities, the frequency and intensity of extreme events related to runoff have increased, such as drought events induced by reduced runoff and flood events caused by rapid increases in runoff. These extreme events have had a catastrophic impact on industrial and agricultural production, people's lives and property, and species diversity. Therefore, reliable runoff forecasts are of great significance for avoiding these potential risks, and the longer the forecast period, the more conducive it is to take measures to mitigate and avoid related potential risks as soon as possible.

[0003] Existing runoff statistical prediction models generally have the problem of poor prediction effect. For example, time series statistical models can only capture the linear relationship between their own variables and cannot well reflect the nonlinear characteristics between variables; machine learning models are mostly gray box or black box models, which are prone to overfitting problems and cannot give clear analytical expressions.

[0004] The introduction and application of Copula function can better solve the defects of the above statistical prediction models. Among them, the Vine Copula function can connect the marginal distribution functions of multiple random variables in the form of Copula pairs to construct a multidimensional joint distribution function, which can characterize the positive and negative dependencies and tail dependencies between different variables. It has been widely used in drought prediction and flood risk analysis research. In high-dimensional situations, the form of the Vine Copula function is closely related to the order of variables, that is, the tree structures corresponding to different variable orders are different, which affects the reliability of the prediction. Bayesian Model Average (BMA) is a collective prediction model that can effectively combine the prediction capabilities of multiple models by assigning different weights to different tree structure prediction models (ensemble members), thereby improving the prediction accuracy.

[0005] As an important part of the hydrological cycle, runoff is affected by many factors, including precipitation, temperature, evaporation, and previous water inflow conditions. At present, the existing models do not consider the influencing factors comprehensively when predicting runoff, which leads to poor prediction accuracy of the models when the forecast period is long, and cannot meet the needs of hydrological forecasting, water resources management, and policy making. Summary of the invention

[0006] The present invention is a runoff prediction method based on Bayesian average combined with vine Copula proposed based on Mach-Zehnder optical path structure. Copula is defined as: a multivariate distribution function formed by connecting one-dimensional marginal distributions in the interval [0,1].

[0007] To achieve the above object, the present invention provides the following technical solution: a runoff prediction method based on Bayesian average combined with vine copula, comprising the following steps:

[0008] Step 1: Analyze the previous precipitation P under different forecast periods t–l (X1), temperature T t–l (X2), Evaporation E t–l (X3) and runoff Q t–l (X4) and later runoff Q t (X5), where t represents the research period; l represents the forecast period in months; X i , i=1,…, 5 represents the corresponding random variables;

[0009] Step 2, according to the minimum chi-square test principle, based on maximum likelihood estimation, select the optimal marginal distribution function for variables X1 to X5 from normal distribution, gamma distribution, Weibull distribution and lognormal distribution respectively and obtain the corresponding cumulative probability;

[0010] Step 4: After fixing the prediction variable X5, for an n-dimensional vine copula function, the structure of the model is closely related to the order of variables. There are n-1 model structures in total. With the order of variables X1, X2, X3, X4 and X5, abbreviated as "12345", it is deduced that the conditional distribution function F(x5|x1,x2,x3,x4) of X5 under the conditions of given X1, X2, X3 and X4 satisfies the following formula (5):

[0011]

[0012] In formula (5), u k represents the random variable x k The cumulative probability of ij is the bivariate Copula parameter corresponding to the jth edge of the ith tree;

[0013] Step 5, based on Akaike AIC information criterion, the optimal Copula function is selected from the bivariate Copula functions and the corresponding parameter θ is obtained;

[0014] Step 6, through recursive calls, use the τ quantile curve τ∈[0,1] to invert formula (5) with respect to x5 to obtain the following formula (6):

[0015]

[0016] In the formula, h –1 and G –1 denote the inverse function of h-function and variable x5 obeying G distribution respectively; Step 7, according to Steps 4 to 6, the function expressions of the model structures corresponding to the other 23 variable orders ("12435", "13245", "13425", "14235", "14325", "21345", "21435", "23145", "23415", "24135", "24315", "31245", "31425", "32145", "32415", "34125", "34215", "41235", "41325", "42135", "42315", "43125", "43215") in the 5-dimensional ViC model can be derived; that is, there are (5–1)! = 24 ViC models in total, and! denotes factorial;

[0017] Step 8: predict the runoff x5 at a specific time in the forecast period t For example, using Monte Carlo simulation to generate 500 uniformly distributed random numbers τ in the interval 0-1, 500 corresponding runoff prediction values ​​x5 can be obtained by using formula (6) and the corresponding ViC model under the other 23 variable orders. t And take the mean of each model structure, so as to obtain the runoff prediction value x5 given by the ViC model under 24 different variable orders t ;

[0018] Step 9: Assume that the 24 prediction models contained in the Vine Copula function in the 5-dimensional case are represented by ViC1~ViC 24 Indicates that ViC=[ViC1,…,ViC 24 ], given the training set D, based on BMA combined with the Vine Copula function, that is, the BMAViC model, the expression of the predicted variable x5 is:

[0019]

[0020] In formula (7), p(ViC m |D) represents the candidate prediction model ViC m The posterior probability can be regarded as p(ViC m |D) statistical weight on the training set D and satisfies Indicates that given a candidate model ViC m The posterior distribution of the predicted value x5 under the condition of and training set D; the expectation maximization algorithm is used to solve the weights of the BMAViC model;

[0021] Step 10: Use the coefficient of certainty R 2 , Nash efficiency coefficient NSE and root mean square error RMSE to evaluate the prediction performance of the model.

[0022] Compared with the prior art, the present invention has the following beneficial effects:

[0023] The present invention constructs a new model for runoff prediction based on Bayesian model averaging (BMA) combined with vine copula (ViC) function, which can effectively combine multiple prediction variables in high-dimensional situations and give clear expressions, which can significantly improve the prediction accuracy of monthly runoff and achieve reliable prediction of monthly runoff under different forecast periods. The invention has high applicability and extensibility, broadens the theoretical framework of runoff prediction, improves the current runoff prediction accuracy, and can provide technical support for the comprehensive construction of water resources management, risk assessment and early warning systems in the basin. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 This is the correlation coefficient diagram between the runoff of the two hydrological stations in the study area and the prediction factors in the 1-3 month forecast period.

[0025] Figure 2 A graphical illustration of the Vine Copula (ViC) model with the variable order of "12345".

[0026] Figure 3 This is a comparison chart of runoff prediction between the BMAViC model and the random forest (RF) model within the 1- to 3-month forecast period within the rate period.

[0027] Figure 4 This is a comparison chart of runoff prediction among the BMAViC model, random forest (RF) model and long short-term memory neural network (LSTM) model within a 1-3 month forecast period during the validation period. DETAILED DESCRIPTION

[0028] The following will be combined with the drawings of the present invention in the embodiments of the present invention to clearly and completely describe the technical solutions of the present invention in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0029] The embodiment of the present invention is described by taking the monthly runoff prediction of two hydrological stations (Tangnaihai and Hongqi) in the upper reaches of the Yellow River as an example.

[0030] Step 1: Collect monthly data on runoff, precipitation, evaporation and temperature from Tangnaihai and Hongqi hydrological stations in the upper reaches of the Yellow River Basin from 1961 to 2016. Note that for hydrological stations without precipitation, evaporation and temperature, grid data containing precipitation, evaporation and temperature (such as the 0.25°×0.25° resolution CN05.1 grid data provided by the National Climate Center) can be considered as a substitute. The specific alternative is to take the area average of monthly precipitation, evaporation and temperature from the grid points contained in the catchment area controlled by a hydrological station to obtain the corresponding monthly value sequence.

[0031] Step 2: Analyze the previous precipitation (P t–l ; X1), temperature (T t–l ; X2), evaporation (E t–l ; X3) and runoff (Q t–l ; X4) and later runoff (Q t ; X5) Figure 1 ), where X i (i=1,…,5) is represented as the corresponding random variable.

[0032] Comparison Figure 1 , the prediction factors and prediction variables of each hydrological station are positively correlated, and the correlation gradually decreases with the increase of the forecast period. Except for the autocorrelation of runoff under the 3-month forecast period, which did not pass the significance test of p<0.05, the other prediction factors and prediction variables showed significant positive correlation (p<0.05) under the 1-3 month forecast period. This shows that the precipitation, temperature, evaporation and runoff selected by the present invention all significantly affect the runoff changes in the later period, that is, there is a close causal relationship between these hydrological and meteorological variables.

[0033] Step 3: According to the minimum chi-square test principle, the precipitation, temperature, evaporation, and runoff of the two hydrological stations in the upper reaches of the Yellow River obey the gamma distribution, normal distribution, Weibull distribution, and log-normal distribution, respectively. Based on the selected marginal distribution function, the cumulative probabilities of variables X1~X5 can be obtained accordingly, which are recorded as u1~u5 respectively.

[0034] Step 4: For a 5-dimensional Vine Copula (ViC) model, after fixing the predictor variable X5, take the variable sequence X1, X2, X3, X4 and X5 as an example, which is abbreviated as "12345" ( Figure 2 ), it is deduced that the conditional distribution function F(x5|x1,x2, x3,x4) of X5 under the given conditions of X1, X2, X3 and X4 satisfies:

[0035]

[0036] In the formula, uk represents the random variable x k The cumulative probability of ij is the bivariate Copula parameter corresponding to the jth edge of the ith tree (see Figure 2 ).

[0037] To obtain formula (5), the present invention first obtains the expression of the distribution function F(x1, ..., x5) of the 5-dimensional continuous random variable X = [X1, ..., X5] according to the Skalar criterion:

[0038] F(x1,…,x5)=C{F1(x1),…,F5(x5)}=C(u1,…,u5) (1)

[0039] Where C represents the 5-dimensional Copula function; u i =F i (x i ) is x i The cumulative probability of .

[0040] Then, the C-vine Copula function is used to combine variables X1 to X5, and the probability density function f(x1,…,x5) of F(x1,…,x5) can be expressed as:

[0041]

[0042] In the formula, f k and c ij|k Represents the variable x k And the probability density function of the two-variable Copula function. For each conditional probability density function f(x|u) in formula (2), it can be written as:

[0043]

[0044] In the formula, u –k It represents the vector after removing the kth variable from the vector u. Introducing the h-function, the expression of F(x|u) is:

[0045]

[0046] Where θ represents the two-variable Copula function Parameters.

[0047] Finally, for a 5-dimensional ViC model, given X1, X2, X3, and X4, when the fixed variable order is X1, X2, X3, X4, and X5, it is abbreviated as "12345" ( Figure 2 ), then according to formula (4), it can be deduced that the conditional distribution function F(x5|x1,x2,x3,x4) of X5 satisfies formula (5);

[0048] Step 5: From the 39 bivariate Copula functions to be selected, including Gaussian Copula, Student t Copula, Archimedean Copula, and Tawn Copula, a total of 39 bivariate Copula functions are selected, and the AIC information criterion is used to calculate the probability of each edge (i.e. Figure 2 The lines with arrows represent the variables that are combined in pairs) select the optimal Copula function type and obtain the corresponding parameters θ.

[0049] Comparison Figure 2 , circles represent variables, each line with an arrow represents an edge (i.e., variables that are connected in pairs), and each subgraph represents a subtree structure diagram of the ViC model under a specific variable order. Similarly, those skilled in the art can obtain tree structure diagrams of the ViC model under other variable orders without creative work.

[0050] Step 6, use the τ quantile curve (τ∈[0,1]) to invert formula (5) with respect to x5 and obtain through recursive calls:

[0051]

[0052] In the formula, h –1 and G –1 They represent the inverse function of the h-function and the variable x5 obeying the G distribution respectively.

[0053] Step 7. Similarly, a person skilled in the art can derive the functional expressions of the model structures corresponding to the other 23 variable orders ("12435", "13245", "13425", "14235", "14325", "21345", "21435", "23145", "23415", "24135", "24315", "31245", "31425", "32145", "32415", "34125", "34215", "41235", "41325", "42135", "42315", "43125", "43215") in the 5-dimensional ViC model according to steps 4 to 6 (i.e., there are a total of (5–1)!= 24 ViC models).

[0054] Step 8: predict the runoff x5 at a specific time in a forecast period t For example, we first use Monte Carlo simulation to generate 500 uniformly distributed random numbers τ in the interval [0,1]. Then, we can obtain 500 corresponding runoff prediction values ​​x5 by using formula (6) and the corresponding ViC model under the other 23 variable orders. tAnd take the mean, finally the present invention can obtain the runoff prediction value x5 given by the ViC model under 24 different variable orders t Note that the x5 obtained by the ViC model under different variable orders in the present invention t Predictions may vary.

[0055] Step 9, in order to reduce the uncertainty of the ViC model prediction corresponding to the single optimal variable order, the present invention performs multi-model ensemble prediction based on Bayesian model averaging (BMA). Assume that the 24 prediction models contained in the Vine Copula function in the 5-dimensional case are represented by ViC1~ViC 24 Indicates that ViC=[ViC1,…,ViC 24 ], given the training set D, the expression for the predictor variable x5 based on BMA combined with the Vine Copula function (i.e., the BMAViC model) is:

[0056]

[0057] Where, p(ViC m |D) represents the candidate prediction model ViC m The posterior probability of can be regarded as its statistical weight on the training set D and satisfies Indicates that given a candidate model ViC m The posterior distribution of the predicted value x5 under the condition of the training set D. The present invention adopts the expectation maximization algorithm to solve the weight of the BMAViC model.

[0058] Step 10, using the coefficient of certainty (R 2 ), Nash efficiency coefficient (NSE) and root mean square error (RMSE) were used to evaluate the prediction performance of the model in the proving period (1963-2006) and the validation period (2007-2016).

[0059] Comparison Figure 3 , during the rate period, the runoff prediction values ​​of the BMAViC model at each hydrological station in the 1-3 month forecast period are in good agreement with their corresponding observation values, and the runoff extremes of each station are well captured, indicating that the BMAViC model can effectively combine the advantages of different ViC models to improve the model's monthly runoff prediction ability. Compared with the random forest (RF) model, the runoff variation process of each station is also well predicted at the 1-3 month forecast period, but the prediction effect of the runoff peak (underestimation) and the runoff in the dry season (overestimation) is poor. The performance evaluation indicators of the BMAViC model at the 1-3 month forecast period are better than those of the RF model, and the good prediction performance can be maintained as the forecast period is extended, indicating that the runoff prediction effect of the BMAViC model is significantly better than that of the RF model.

[0060] Comparison Figure 4 , it can be seen that during the validation period, the linear fitting curves of the RF model or LSTM model deviate significantly from the 45° diagonal line and tilt to the right as the forecast period increases ( Figure 4 The black dotted line or gray dashed line fitting curve in the middle), especially for the 3-month forecast period. The scatter points between the runoff prediction values ​​of the BMAViC model and the observed values ​​are evenly distributed near the 45° diagonal line ( Figure 4 The black “+” mark and the black dotted line fitting curve in the middle indicate that the BMAViC model has strong runoff prediction capabilities at each hydrological station within a 1-3 month forecast period, and is significantly better than the RF model and the LSTM model ( Figure 4 top left text).

[0061] The above implementation modes are only used to illustrate the present invention rather than to limit the present invention. Ordinary technicians in the relevant technical field may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also belong to the scope of the present invention. The patent protection scope of the present invention is defined by the claims.

[0062] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A runoff prediction method based on Bayesian averaging combined with vine copula, characterized in that: The steps include: Step 1: Analyze the previous precipitation P under different forecast periods t–l (X1), temperature T t–l (X2), Evaporation E t–l (X3) and runoff Q t–l (X4) and later runoff Q t (X5), where t represents the research period; l represents the forecast period in months; X i , i=1,…,5 represents the corresponding random variable; Step 2, according to the minimum chi-square test principle, based on maximum likelihood estimation, select the optimal marginal distribution function for variables X1 to X5 from normal distribution, gamma distribution, Weibull distribution and lognormal distribution respectively and obtain the corresponding cumulative probability; Step 3, after fixing the prediction variable X5, for an n-dimensional vine copula function, the structure of the model is closely related to the order of variables. There are n-1 model structures in total. With the order of variables X1, X2, X3, X4 and X5, abbreviated as 12345, it is derived that the conditional distribution function F(x5|x1,x2,x3,x4) of X5 under the conditions of given X1, X2, X3 and X4 satisfies the following formula (5): In formula (5), u k represents the random variable x k The cumulative probability of ij is the bivariate Copula parameter corresponding to the jth edge of the ith tree; Step 4, based on Akaike AIC information criterion, the optimal Copula function is selected from the bivariate Copula functions and the corresponding parameter θ is obtained; Step 5, through recursive calls, use the τ quantile curve τ∈[0,1] to invert formula (5) with respect to x5 to obtain the following formula (6): In the formula, h –1 and G –1 They represent the inverse function of h-function and variable x5 obeying G distribution respectively; Step 6, according to steps 3 to 5, derive the function expressions of the model structures corresponding to the other 23 variable orders 12435, 13245, 13425, 14235, 14325, 21345, 21435, 23145, 23415, 24135, 24315, 31245, 31425, 32145, 32415, 34125, 34215, 41235, 41325, 42135, 42315, 43125, 43215 in the 5-dimensional model; That is, there are (5–1)! = 24 models, where! represents factorial; Step 7: predict the runoff x5 at a specific time in the forecast period t For example, Monte Carlo simulation is used to generate 500 uniformly distributed random numbers τ in the interval 0-1. 500 corresponding runoff predictions x5 are obtained by formula (6) and the corresponding models under the other 23 variable orders. t And take the mean of each model structure to obtain the runoff prediction x5 given by the model under 24 different variable orders t ; Step 8: Assume that the 24 models contained in the Vine Copula function in the 5-dimensional case are represented by ViC1~ViC 24 Indicates that ViC=[ViC1,…,ViC 24 ], given the training set D, based on BMA combined with the Vine Copula function, that is, the model, the expression of the predicted variable x5 is: In formula (7), p(ViC m |D) represents the model ViC m The posterior probability is regarded as p(ViC m |D) statistical weight on the training set D and satisfies p m (x5|ViC m ,D) indicates that in the model ViC m The posterior distribution of the predicted value x5 under the condition of and training set D; the expectation maximization algorithm is used to solve the weight of the model; Step 9: Use the coefficient of certainty R 2 , Nash efficiency coefficient NSE and root mean square error RMSE to evaluate the prediction performance of the model.

Images